Solving Diffusion Problem Crank Nicholson Scheme

1. Solving Diffusion Problem Crank Nicholson Scheme The 1D Diffusion Problem is: John Crank Phyllis Nicolson 1916 –2006 1917 –1968 Here the diffusion constant is a function of T: We first define a function that is the integral of D: Or equivalently, with constant f = 5/7.

2. Insert z back into the diffusion equation we find: The RHS is then turned into the Laplacian of function z: Now split both sides at the nth time step, we get: If we substitute z back we get an expression for the function value at n+1 time step ... ? But wait…

3. Let’s look at the stability requirement. By simple algebra we get the time step for this scheme to be stable is: This is roughly the diffusion time across one cell. A domain containing hundreds of cells would require enormous time for the diffusion to cause some noticeable effect. Thus we need some scheme that allows us to take larger time steps but retains stability. The trick here is to link the function value at n-th time step with n+1 –th time step. We write the differencing scheme as: This graph shows how this scheme works:

4. It can be proven that by using this implicit method, the scheme becomes unconditionally stable for any step size chosen. Now let’s do the back substitution. It should be: Instead of stability issues, it requires us to solve a complicated nonlinear equation system at each time step. But there’s a way to get around this problem. We can write the integral of D at the n+1 –th time step, which causes the nonlinearity as: Now we make the n+1 dependence on the RHS linear. Do the substitution, we get:

5. where We define the coefficients as above then the equation becomes: The problem is then transformed into solving a tridiagonal matrix:

6. * taken from Numerical Recipies This matrix can easily be solved by back substitution scheme, which can be found in any elementary numerical method book. Go 2D Now consider the 2D diffusion problem. The 2D Crank-Nicholson scheme is essentially the same as the 1D version, we simply use the operator splitting technique to extend the method to higher dimensions. Explicitly, the scheme looks like this:

7. Step 1. evolve half time step on x direction with y direction variance attached where

8. Step 2. evolve another half time step on y direction with x direction variance attached. where Since the sweeps on different directions are identical, it is possible to solve a multidimensional diffusion problem by a single subroutine.

9. In solving Euler equation with diffusion, we can use operator splitting: solve the usual Euler equation by splitting on different directions thru time step dt to get the density, velocity and pressure. Then from the pressure at each grid, find the temperature distribution, do a Crank-Nicholson calculation with the same time step dt (here we still need to split dt) to find a new temperature distribution. Then use this temperature distribution to find pressure distribution again. This is the actual pressure distribution at dt. Then go on to the next time step.