Diffusion by Numerical Analysis

Ever been overwhelmed by the mathematics of diffusion? Stop looking for your error function table and grab your personal computer. All you need is a spreadsheet program, a little knowledge of numerical analysis and you can solve all sorts of complicated diffusion problems and that includes problems that do not have analytical solutions.

Most diffusion problems can be solved using a finite difference method of numerical analysis and Fick's first law of diffusion. Consider the solute profile shown in figure 1 and the change in composition at C1 produced by a flux into ( Jin) and a flux out of (Jout) a region of thickness h. The flux equations would look like the following:

J_in= [q_in/delta(t)] = - D [(C1 - Co)/h]

J_out= [q_out/delta(t)] = - D [(C2 - C1)/h]

where delta(t) is a finite time interval, D is the diffusivity of the solute, q represents the amount of solute per unit area entering or leaving the volume of interest and Ci is the solute per unit volume. The net amount of solute collected in time delta(t) will be

C₁^new - C₁ = [(q_in - q_out)/h] = [(D x delta(t))/h^2 x (C_o + C₂ - 2C₁)

This equation can be greatly simplified by setting [D x delta(t)/h^2] equal to 1/2 to give

C₁^new = (C_o + C₂/2

The only constraint is that the time interval and distance between composition measurements are fixed. This requirement is called the stability criterion for the finite difference method.

Consider a carburization problem where the surface of a 0.22 weight percent (wt.%) carbon steel is held in a carburizing atmosphere. The carbon potential of the atmosphere is 1.1 wt.% at 925°C (1700°F) and the carbon profile is required after 8 hours. To solve this problem a spreadsheet would be set up as shown in figure 2. Only a portion of the spreadsheet is shown to indicate the simplicity of the technique and how each successive column is calculated. Each time interval (gt) requires a new column to be added and the time interval is set by the choice of distance interval (h) and the diffusivity (D). For a diffusivity of 1.59x10^-5 mm^2/sec, the gt time interval would be 314 seconds.

Thus, 92 columns or iterations of Æt would be required to calculate the profile for an 8 hour carburization cycle. To simplify the spreadsheet, the depth was run to 6 mm and the composition at 6 mm was kept constant throughout the 92 iterations at a value of 0.22 wt.% C. This is equivalent to a 6 mm sheet of steel that is held between two different atmospheres: one side of the steel sheet at a carbon potential of 1.1 wt.% and the other side at a potential of 0.22 wt.%. This constrains the carburization cycle to less than 10 hours. Longer times can be run by increasing the spreadsheet to incorporate greater depths, but the last row will still have to be held constant at the base carbon composition. A plot of the resulting carbon profile for the 8 hour cycle is shown in figure 3 along with the analytical solution using a mathematical approximation for the error function.

With a little effort, the spreadsheet can be turned into a powerful tool where carbon potential, base carbon content and diffusivity can be entered as variables to test new carburization cycles. Carbon leveling operations can also be modeled by lowering the surface carbon potential after the case depth is obtained. In the problem shown here, the time to attain a uniform carbon concentration, less than 0.9 wt.%, can be determined by changing the surface carbon concentration to 0.9 in column 93 and then inspecting the subsequent columns to determine the number of intervals required to level the carbon. Many spreadsheets are equipped with graphics packages and the effect of changing the variables can be immediately observed in these graphs at fixed number of iterations or times.