Modeling Diffusion

 

Let us now consider diffusion as the driving force of change in the concentration in our system. The force that makes the stuff move in this case is the difference between concentrations in adjacent segments. It is also good to remember that in this discrete approximation we are actually dealing with points on a continuum, in this case a line 0x. The concentrations that we are considering are located at these points. We are dealing with average concentrations for the whole segments and we are assuming that these averages are located at these nodes. Therefore if there is no outside force to move the stuff, it would be reasonable to assume that the farther away the points we consider are, the less stuff can be moved between them by the concentration gradient. .

Just like before let us define the concentration of stuff in any given segment at time t+Dt, assuming that we know the concentration there at time t. Let us write the equation for the total amount of stuff in segment at time t+Dt:

C(t+Dt, x). Dx =
      C(t, x). Dx + [C(t, x-Dx) - C(t, x)]/ Dx .D. Dt + [C(t, x+Dx) - C(t, x)]/ Dx .D. Dt

In this equation [C(t, x+Dx) - C(t, x)]/ Dx is the empirically derived equation for the diffusive flux between two adjacent segments. D is the diffusion coefficient that characterizes the environment, the media; it tells us how fast diffusion can occur in this kind of media.
After some rearranging we get:

Once again, if we let Dx ➔ 0 and Dt ➔ 0, we get the well known diffusion equation as a partial differential equation:

In discrete notation, the equation for concentration at the next time step becomes:

C(t+Dt, x) = C(t, x) + [C(t, x-Dx) - 2C(t,x) + C(t, x+Dx)].D .Dt / Dx2.        (2)

If we know the concentration at the previous time step we can calculate the concentration at the next time step. Just like in the Advection example, to calculate this equation at any (x,t) we need to define the initial condition:

C(0,x) = co(x)

As for the boundary conditions, in this case we will need two of them. We cannot use equation (2) to calculate both the value on the left hand side boundary C(t,0) and on the right hand side boundary C(t,N), where N is the number of the maximal segment that we consider. Therefore we need two boundary conditions:

C(t,0) = b1(t);       C(t,N) = b2(t).

Similarly, there may be other types of boundary conditions, such as:

C(t,0) = C(t,1),      C(t,N-1)=C(t,N).

This will be a condition of no flow across the boundaries.