Numerical Methods in Diffusion

1. Numerical Methods in Diffusion J.C. LaCombe University of Nevada, Reno Reno, NV, USA lacomj@unr.edu These lecture notes complement an online learning module and diffusion simulation software that can be found at: http://unr.edu/homepage/lacomj/Diffusion/index.htm Portions of this lecture were adapted from Elements of Heat and Mass Transfer, 3rd ed., F.P. Incropera, and D.P. De Witt John Wiley & Sons, NY, 1990

2. Introduction • In many real-world problems, the details of the system may not correspond to a known solution to the diffusion equation. In these cases, it is usually possible to produce a solution where the equation is solved through iterative techniques. This is generally a lot of work to do by hand. However, with the aid of a computer, this becomes possible. The techniques used to reach such solutions are known as numerical methods. • Situations that often require numerical methods include • Multi-dimensional problems (not simply 1-D) • Complex Problems • Complex Geometry/Shape • Complex boundary conditions • Complex initial conditions Numerical Methods

3. Discretization of the Problem The numerical methods we will use in this course solve complex diffusion problems by breaking up the system into manageable parts and solving the diffusion equations for each part simultaneously with all the other parts. To do this, we discretize the problem mathematically by dividing up space into little elements or nodes and treating time as moving forward in small steps. Elements and Nodes Component divided into elements Element Node

4. Discretization of the Problem Each element is identified using subscripts, and has a finite dimensions. Discretization of the Problem Nodem,n

5. The Finite-Difference Approach Recall: Fick’s 1st law simply tells us how solute will flow if there is a concentration gradient. Recall: Fick’s 2nd law is simply a combination of Conservation of Mass with Fick’s 1st law. Fick’s Laws in 2-D Fick’s 2nd Law in 2-D Cartesian Coordinates. (1) Before we solve this, we need to re-write the equation into a discretized form. The approach presented here is known as a finite-difference solution approach.

6. The Diffusion Equation (1) Discretizing the Spatial Derivatives Consider first, the spatial 2nd derivatives on the RHS of Equation (1). The 2nd derivatives can be thought of more simply as the slope of the 1st derivatives. This is written (approximated) in the x-direction as Gradient at RHS of CV Gradient at LHS of CV (2) Slope of the 1st derivative Note that the 1st derivatives here are simply the concentration gradients. We can use the central difference approximation to determine these…

7. (2) The Diffusion Equation Discretizing the Spatial Derivatives m-1,n+1 m,n+1 m+1,n+1 Eq. (2) further simplifies if we can determine the concentration gradient at the midpoint between nodes. The gradient can be estimated using the concentration values at the neighboring nodes and the distance between the nodes. Thus, m-1,n m,n m+1,n Evaluate gradient here (3) C(x) (4) m-1 m m+1

8. The Diffusion Equation Equations (3) and (4) can now be substituted into (2) to produce the discretized form of the 2nd spatial derivative in the x direction. Discretizing the Spatial Derivatives (5a) The 1st and 2nd derivatives in the y (and z) directions can also be evaluated in a similar manner… (5b) These make up the RHS of Fick’s 2nd Law (Eq. 1)

9. The Diffusion Equation Now that we have discretized Fick’s 2nd law in space (1), we must discretize it in time as well. To do this, we will introduce a new variable, p, that is an integer that represents the time step. The duration of each step is Dt. Thus, the total time, t, is written… Discretizing the Time Derivative (6) The finite-difference approximation to the time derivative (the LHS of Fick’s 2nd law) is then expressed as… (7) This is the LHS of Fick’s 2nd Law (Eq. 1) The superscript, p, denotes the time dependence. The time derivative is expressed in terms as the difference in concentrations between the new time (step p+1) and the previous time (step p).

10. Fick’s 2nd Law in Discretized Form We present here a solution approach known as the explicit method. In this finite-difference scheme, the concentration at any node m,n at time t+Dt is calculated from knowledge of the concentration at the same and neighboring nodes for the preceding time t. We now can combine (5a,b) and (7) to produce the discretized form of the diffusion equation, (1). The 2-D Diffusion Equation Fick’s 2nd Law in discretized form... (8) Note: Other approaches, such as the implicit method, are more efficient with a computer, but require more complex algorithms. Nonetheless, the fundamental principles are the same as we are applying here. We will not be covering these other methods in this course.

11. The Fourier Number Equation (8) is the general form of our solution. We can simplify the notation a bit if we use square elements, so that Dx = Dy. Additionally, we can form the following group of parameters, which is commonly known as the dimensionless Fourier Number, Fo. The Diffusion Equation (9) Now, we can re-arrange (8) to solve for the concentration in node m,n at the new time step, p+1. This equation applies to any element/node on the interior of a component. The expression simplifies to… 2-D Interior Node (10) The NEW composition in an element is calculated using the PREVIOUS compositions in the element and its neighbors.

12. Explicit Method & Solution Stability For the case of 1-D transport, Equation (8) would instead develop into the form of The 1-D Diffusion Equation 1-D Interior Node (11) The accuracy of finite-difference solutions may be improved by decreasing the values of Dx and Dt (I.e., finer discretization). On the other hand, making these values larger will allow the calculation to proceed more quickly. One additional limitation of the explicit method is that it is not always a stable solution. If the values of Dx and Dt are not small enough, it can cause the solution to oscillate (even when this is physically impossible).

13. Stability Criteria To prevent such erroneous results when the solution is “unstable”, the values of Dx and Dt must meet certain criteria (details omitted). For interior nodes, these are, Critical Values of Fo Recalling, 1-D Stability Criteria 2-D Stability Criteria So, once you pick a value of either Dx or Dt , the other value must be chosen so that the stability criteria is met. Simply re-arrange the equation for Fo to calculate the acceptable value.

14. m-1,n+1 m,n+1 External Surface(zero-flux plane) m-1,n m,n m-1,n-1 m,n-1 m+1 Imaginary node m-1 m Zero-Flux Elements and Surfaces The equations presented so far (10, 11) are for interior nodes. I.e., each element’s surroundings are geometrically the same in all directions. We can develop similar equations for different element types, but we need to be clever, or it gets messy. A surface node (with no flux flowing through the surface), can be modeled using the same equation as an interior node. All we need to do is include an imaginary node just outside the surface and set its composition to the same as the node just inside the surface. This has the effect of producing a zero net gradient through the surface. I.e., if the surface node is then the no-flux condition is modeled by adding an imaginary node at m+1 and setting to achieve a state of no-flux at node m (no gradient means no net flux). Other Element Configurations

15. Zero-Flux Elements and Surfaces So, we can model a surface node by modifying the equation for an interior node. Recalling Equation (10) for 2-D, Other Element Configurations 2-D Interior Node (10) We then incorporate the imaginary node… And are left with the equation for a surface node… 2-D Surface Node (12) In 1-D, this would work out to… 1-D Surface Node (13)

16. Other Solution Approaches The method used on the previous slides to discretize the problem is not the only way to produce equations such as Equations (10)-(13). Another method can be used to provide even greater flexibility with boundary conditions. It is simply based on conservation of mass (Recall that Fick’s 2nd law is also essentially this as well). Let us consider the element surrounding each node to be subject to conservation of mass. This would be written as… Developing Expressions for OtherElement Types Solid State Diffusive Flux Solute “Generated” Stored Mass + = In practice, this can be something like solute entering an element at external surface

17. Other Solution Approaches Writing this for a generic interior node, we account for all possible influences. As before, minor changes can be made for an external node. Note that here, flux into the node is considered “positive”. Other Ways to Discretize the Problem (14) Where,

18. Diffusion with “Mass Generation” When massaged, Equation (14) evolves into the same form as the earlier equations (10)-(13), except now, we have added in the solute generation term. Other Ways to Discretize the Problem (15) And in 1-D, this is (16) Thus, there are a variety of approaches to produce the discretized diffusion equations for a variety of different element types.

19. 1 2 3 4 5 6 Dx110-3 cm Example An earlier topic presented the analytical solution to the case of a binary diffusion couple. Let’s analyze this using a finite-difference model. Assume D = 110-9 cm2/s, and the initial compositions are Cl= 0.75, and CR = 0.25. 1-D Thick Diffusion Couple First, the stability criteria for this 1-D arrangement is that Fo ½. Thus the maximum time step for a stable solution of this problem is 500 seconds.

20. Example (1-D Interior Node) Equation (11) is the suitable solution form: 1-D Thick Diffusion Couple To handle the “infinite” ends, we treat them as having no-flux conditions (I.e., the concentration gradient is zero at the ends) using Equation (13). This will be ok, provided that the concentration field never reaches the end during our simulation. The equations are written for each of the 6 elements… These 6 equations must be solved at each time step, p. There will be one equation for each node. Models with lots of elements involve solving lots of equations.

21. Example This is expressed more concisely in matrix form. 1-D Thick Diffusion Couple The new concentrations at each node are calculated by solving this matrix at each time step. At each step, you use the resulting concentrations from the previous step, Cp+1, as the new values of Cp. Likewise, to get it all started, you just use the initial concentrations at each node. You can use whatever methods or software you want to solve the matrix. Even a spreadsheet will work…

22. Example MS Excel Worksheet… 1-D Thick Diffusion Couple Double-click above (ppt only) to open the actual spreadsheet!

23. m,n+1 m,n+1 m,n m,n m-1,n m+1,n m-1,n m+1,n m,n-1 m,n-1 m,n+1 m,n+1 m,n m,n m-1,n m+1,n m-1,n m+1,n m,n-1 m,n-1 2-D Equation Summary Plane Surface Interior Node 2-D Explicit Finite Difference Equations (Dx=Dy) Interior Corner Exterior Corner