Massachusetts Institute of Technology
Department of Materials Science and Engineering
77 Massachusetts Avenue, Cambridge, MA 02139-4307

3.21 Kinetic Processes in Materials --- Spring 2006

Some numerical work with diffusion boundary-value problems

Miguel Marioni, DMSE Ph.D. 2003, made this impressive Mathematica worksheet while he was TA in 3.21, adapting a MathCad worksheet I had previously developed with similar content.

Some of the coding in this notebook is pretty sophisticated; the main thing to be grasped is the behavior of the various forms of the solutions to the diffusion equation that are evaluated and plotted.  I encourage you to experiment with changing values of variables, etc., to get comfortable with how the equations depict the physical problem being represented.

Those of you who are learning Mathematica and using it for solving a variety of problems may find some useful sections of code in this notebook that can be adapted for other purposes.

A word on Mathematica

Executable inputs to Mathematica receive a label of the form In[number], where number is a counter for each expression. Each  In[number] statement will be followed by a corresponding  Out[number] right after it, unless printing of the output is suppressed by appending a colon ";" to the end of each In[number] statement.

A couple of hints

The infinity sign (∞) is typed in from the pallette, or otherwise pressing the sequence of keys:
<ESC> + "inf"  +<ESC>   (Overscript[⟶, produces]) ∞
Works the same for other symbols:
<ESC> + "a"  +<ESC>   (Overscript[⟶, produces]) α
<ESC> + "b"  +<ESC>   (Overscript[⟶, produces]) β
One particular example:
<ESC> + "ee"  +<ESC>   (Overscript[⟶, produces]) e   (which is the same as capital "E", and is a shortcut for the Exp[] operator.

Note: like capital "E", other letters are reserved. Such is the case of "I", equivalently written i (<ESC>+"ii"+<ESC>) which stands for the imaginary constant, and (to our great grief) "D", which is the derivation operator. To circumvent this inconvenience to our usual symbol for diffusion, we use the calligraphic D,  <ESC>+"scD"+<ESC> ⟶𝒟.

Function arguments are given between brackets "[" and "]" and are separated by commas ",".  In the definition, the functions take named patterns really (name followed by an underscore "_"), so for instance one would define the function f[x_] := (x^2)   A statement like  f[x_:  1.234] := (x^2)  defines the same function but using 1.234 as the default value for x.   
For more information look in Mathematica's Help Browser!

Some default packages and definitions for this notebook

<< Graphics`Legend`
<< Calculus`FourierTransform`
<< Graphics`MultipleListPlot`

The following lines set options that define styles for all plots in this notebook.

myt = AbsoluteThickness[2];
Color = {PlotStyle -> {{myt, Hue[0]}, {myt, Hue[ .2]}, {myt, 
                Hue[ .4]}, {myt, Hue[ .6]}, {myt, Hue[1]}}, 
          LegendShadow -> None, 
          LegendPosition -> {1, - .2}};
Black = {PlotStyle -> {AbsoluteDashing[{}], 
              AbsoluteDashing[{2, 2}], AbsoluteDashing[{4, 4}], 
              AbsoluteDashing[{4, 2, 1, 2}], AbsoluteDashing[{2, 2, 2, 4}]}, 
          LegendShadow -> None, 
          LegendPosition -> {1, - .2}};

\(SetOptions[Plot, 
      AxesLabel -> {"x [m]", \*""c [\!\(m\^\(-3\)\)]""}, 
      DefaultFont -> {"Times", 14}]; 
    SetOptions[MultipleListPlot, 
      AxesLabel -> {"x [m]", \*""c [\!\(m\^\(-3\)\)]""}, 
      DefaultFont -> {"Times", 14}, PlotJoined -> False, 
      SymbolShape -> PlotSymbol[Box], 
      SymbolStyle -> {Hue[0], Hue[0.2`], Hue[0.4`], Hue[0.6`], Hue[1]}, 
      PlotStyle -> {Hue[0], Hue[0.2`], Hue[0.4`], Hue[0.6`], 
          Hue[1]}];\)
\(\(SetOptions[ListPlot, PlotJoined -> False, 
        PlotStyle -> AbsolutePointSize[5], 
        DefaultFont -> {"Times", 14}, 
        AxesLabel -> {"x [m]", \*""c [(m^(-3))]""}];\)\)

Now check to see that they have been set properly…

Options[ListPlot]

equation and rule manipulation

ClearAll[S, M, Op];

    
S[a_ < b_, c_ < d_] := a + c < b + d
S[a_ ≤ b_, c_ ≤ d_] := 
      a + c ≤ b + d
S[a_ > b_, c_ > d_] := a + c > b + d
S[a_ ≥ b_, c_ ≥ d_] := 
      a + c ≥ b + d
S[a_ == b_, c_ == d_] := 
      a  + c == b + d

        
S[a_ < b_, c_] := a + c < b + c
S[a_ ≤ b_, c_] := 
      a + c ≤ b + c
S[a_ > b_, c_] := a + c > b + c
S[a_ ≥ b_, c_] := 
      a + c ≥ b + c
S[a_ == b_, c_] := 
      a  + c == b + c

        
M[a_ < b_, c_] := a c < b c
M[a_ ≤ b_, c_] := 
      a c ≤ b c
M[a_ > b_, c_] := a c > b c
M[a_ ≥ b_, c_] := 
      a c ≥ b c
M[a_ == b_, c_] := 
      a c == b c

        
Op[a_ < b_, c_] := c[a] < c[b]
Op[a_ ≤ b_, c_] := 
      c[a] ≤ c[b]
Op[a_ > b_, c_] := c[a] > c[b]
Op[a_ ≥ b_, c_] := 
      c[a] ≥ c[b]
Op[a_ == b_, c_] := c[a] ==  c[b]

L[a_ < _] := a 
L[a_ ≤ _] := a
L[a_ > _] := a
L[a_ ≥ _] := a
L[a_ == _] := a

    
R[_ < a_] := a 
R[_ ≤ a_] := a
R[_ > a_] := a
R[_ ≥ a_] := a
R[_ == a_] := a

RS[Rule[a_, b_]] := Rule[a, Simplify[b]]
RS[Rule[a_, b_], op_] := Rule[a, op[b]]

Diffusion couple: two semi-infinite solids

Error function

Mathematica has a built in definition for the error function, \(TraditionalForm\`Erf[
          z] = \(2\/\@π\) \(∫\_0\%z\( e\^\(-u\^2\)\) du\)\).  Some of it's properties are exacted below:

Erf[0]

Erf[∞]

Erf[-∞]

Erf[-z]

Diffusion couple comprised of two semi-infinite bars

Consider two semi-infinite solids joined at x = 0,  and a constant and uniform diffusivity (Subscript[𝒟, 1]):

(Subscript[𝒟, 1]) = (10^(-16));

(we have suppressed the units, (m^2)/s for convenience but we should keep in mind that they are implicit.) We know that the solution to the diffusion equation with constant diffusivity 𝒟 for the step composition profile:

c[x, t = 0]={({{0}, {(Subscript[c, o])}})     ({{if x < 0}, {if x > 0}})

is the following c1[x,t]:

c1[x_, t_, 𝒟_:  (Subscript[𝒟, 1]), 
        c_:  (Subscript[c, o])] := (c/2) (1 + 
              Erf[(x/(Sqrt[4  𝒟 t]))]) /; 
        t > 0
c1[x_, t_, 𝒟_:  (Subscript[𝒟, 1]), c_:  (Subscript[c, o])] := 
      If[x > 0, c, 0] /; t ≤ 0

where the second line is the initial condition, written explicitly so as not to generate division by zero errors when t = 0, and where in this case we define

(Subscript[c, o]) = 1

We can plot the solutions for different times:

\(\(Plot[{c1[x, 0], c1[x, 10\^1], c1[x, 10\^2], 
          c1[x, 10\^3]}, {x, \(-\ 10\^\(-6\)\), 10\^\(-6\)}, Color, 
        PlotLegend -> {"t=0", \*""t=(10^1)s"", \*""t=(10^2)s"", \*""t=(10^3)s""}];\)\)

Diffusion from a plate of finite thickness

Uniform concentration

Consider a semi-infinite slab of thickness L = (Subscript[L, o]), at uniform concentration (Subscript[c, o]) of a diffusant. The diffusant can leave the material at it's boundaries, so that at x = 0 and x = L the concentration is always 0.

With:

(Subscript[L, o]) = (10^(-6));

where again, we left out the units which are implicitly MKS, we obtain for the solution c2[x,t] (see Kinetics of Materials, Eq. (5.49):

c2[x_, t_, n_:  10, 𝒟_:  (Subscript[𝒟, 1]), 
        c_:  (Subscript[c, o]), 
        L_:  (Subscript[L, o])] := ((4  
                c)/π) (∑_(j = 0)^n) (1/(2  j + 1)) 
              Sin[((2  j + 
                        1)/L) π x] (E^(-(((((2  j + 1)/L) π))^2) 𝒟 t)) /; 
        t > 0
c2[x_, t_, 𝒟_:  (Subscript[𝒟, 1]), c_:  (Subscript[c, o]), 
        L_:  (Subscript[L, o])] := 
      If[0 < x && x < L, c, 0] /; t ≤ 0

where n = ∞, but we will generally use a sufficiently large finite number, in this case 10. (More terms may be necessary for shorter times. The effect of the number of terms is shown on the next subsection.) We can plot the solutions for different times:

\(\(Plot[{c2[x, 0], c2[x, 10\^1], c2[x, 10\^2], c2[x, 10\^3]}, {x, 0, 
          L\_o}, PlotRange -> {0, c\_o}, Color, 
        PlotLegend -> {"t=0", \*""t=(10^1)s"", \*""t=(10^2)s"", \*""t=(10^3)s""}];\)\)

The effect of the number of terms in the summation:

We now plot the profile as above for time t=10s, but taking different number of terms for the summations. More terms give a more accurate result. About ten terms are sufficient for this particular case. If the time was shorter, more terms in the summation would be necessary.

Plot[{c2[x, (10^1), 10], c2[x, (10^1), 5], c2[x, (10^1), 3]}, {x, 0, 
          (Subscript[L, o])}, PlotRange -> {0, 1.2  (Subscript[c, o])}, Color, 
        PlotLegend -> {"n=10", "n=5", "n=3"}];

Arbitrary initial concentration

Now look at the problem with arbitrary c[x, 0], in this example a linear function, and again solve for the same boundary conditions (at x = 0 and x = L the concentration is always 0.) So we just use the Fourier sine-series:

With

(Subscript[c, 1]) = 0.5;
(Subscript[c, 2]) = 1;

we write the initial condition as a linear function of x:

c3[x_, t_, n_:  10, L_:  (Subscript[L, o]), cL_:  (Subscript[c, 1]), 
        cR_:  (Subscript[c, 2]), 𝒟_:  (Subscript[𝒟, 1])] := 
      If[0 < x && x < L, cL + (cR - cL) x/L, 0] /; 
        t ≤ 0

and find its Fourier coefficients  (see Kinetics of Materials, eq. (5.45):

A[n_] := ((2/(Subscript[L, o])) (∫_0^((Subscript[L, o])))c3[x, 0] 
              Sin[((n π x)/(Subscript[L, o]))] \[DifferentialD]x)

There is a little numerical problem for (Subscript[L, o]) = (10^(-6)) m with the precision loss during the numerical integration (a computational problem -- note that other programs wouldn't even report this):

A[30]

There are several tactics to avoid these problems:  in this simple case we could integrate analytically (done below in Ab[n]). We could instead change variables u = n x/(Subscript[L, o]) which circumvents the difficulty and is more general (done below in Ac[n]) or we could turn off the warnings with
Off[FactorSquareFree::"lrgexp"];
and hope for the best (but we want to check that the result is okay, and will therefore compare the coefficients below):

Ab[n_] := ((2/(Subscript[L, o])) (∫_0^((Subscript[L, o])))c3[x, 0] 
                Sin[((nn π x)/(Subscript[L, o]))] \[DifferentialD]x) /. 
        nn -> n

Ac[n_] := 
      N[(2/n) (∫_0^n) c3[((u (Subscript[L, o]))/n), 0] 
                Sin[u π] \[DifferentialD]u] /; n > 0

A comparison of coefficients follows: (NOTE The table below likely extends off the right-hand side of your screen and you will have to scroll over to see that part if this is of interest to you.)

Off[FactorSquareFree::"lrgexp"];
Prepend[
        Table[{n, A[n], Ab[n], Ac[n]}, {n, 1, 
            30}], {"n", "A[n]", "Ab[n]", "Ac[n]"}] // 
      TableForm
On[FactorSquareFree::"lrgexp"];

We finally define the composition

c3[x_, t_, n_:  10, L_:  (Subscript[L, o]), cL_:  (Subscript[c, 1]), 
        cR_:  (Subscript[c, 2]), 𝒟_:  (Subscript[𝒟, 1])] := (∑_(k = 1)^n) Ac[k] 
            Sin[((k π x)/L)] (E^(-((((k/L) π))^2) 𝒟 t)) /; t > 0

We can plot the solutions for different times. We have the problem of underflow, however, when the argument of the exponential (T[t])

\(\(Plot[{c3[x, 0], c3[x, 10\^1, 20], c3[x, 10\^2], c3[x, 10\^3]}, {x, 0, 
          L\_o}, PlotRange -> {0, c\_o}, Color, 
        PlotLegend -> {"t=0", \*""t=(10^1)s"", \*""t=(10^2)s"", \*""t=(10^3)s""}];\)\)

Finite-difference method of solving the diffusion equation

Finite-difference methods are useful for solving all sorts of diffusion problems when closed-form analytical solutions are not available.  Here we set out the basic method for representing derivatives in terms of a discrete representation of the composition profile, that is, the value of c is known at a set of discrete points and these data are used to approximate derivatives and perform an iterative calculation to project the evolution of the concentration field.
We do this below for one-dimensional diffusion and ultimately apply the result to diffusion in a finite diffusion couple with a concentration-dependent diffusivity.

We now consider the case of a diffusion couple of pure A and pure B. Let each piece be 1 μm thick:

(Subscript[L, 1]) = (10^(-6));

and pick a time short enough so that they are effectively semi-infinite. The maximum time will be 10 s.

We will tackle this problem with finite differences.

Finite differences

We start by expanding an arbitrary (𝒞^1) function u at x + Δx in powers of Δx.

f1p = 
      u[x + Δx] == 
          Series[u[x + Δx], {Δx, 0, 2}] // 
        Normal
f1m = 
      u[x - Δx] == 
          Series[u[x - Δx], {Δx, 0, 2}] // 
        Normal

We now construct the centered differences by solving the above equations for u'[x] and u''[x].

{Solve[S[f1p, M[f1m, -1]],u'[x]],
Solve[S[f1p, f1m],(u')'[x]]}//Flatten

For the time domain, we take the forward difference, so that we can refer the concentration at time t + Δt to that of time t. For consistency we use a similar derivation to the above one, but will end up obtaining the very expression of the derivative, before taking the limit Δt → 0.

u[t + Δt] == 
        Series[u[t + Δt], {Δt, 0, 1}] // 
      Normal

Solve[%, u'[t]] // Flatten

Putting these results together:

fd={Derivative[1, 0][c][x, t]->((c[x + Δx, t] - 
                    c[x - Δx, t])/(2 Δx)),
Derivative[2, 0][c][x, t]->((c[x + Δx, t] + 
                    c[x - Δx, t] - 
                    2  c[x, t])/Δ)(x^2),
Derivative[0, 1][c][x, t]->((c[x, t + Δt] - 
                    c[x, t])/Δ)t};

Diffusion equation with finite differences

The diffusion equation can now be written in its full non-linear form in one dimension as

eq1 = (∂_t) c[x, 
              t] == (∂_x)(𝒟[
                  c[x, t]] (∂_x) c[x, t]) // Simplify

We turn to the finite difference version:

eq2 = eq1 /. fd

And solve for the concentration at point t + Δt.

eq3 = Solve[eq2, c[x, t + Δt]] // Flatten

Solution for concentration-dependent diffusivity 𝒟2 (quadratic in c)

Assume 𝒟2 varies parabolically with c and increases by an order of magnitude as c goes from 0 to 1

𝒟2[c_] := 
        (10^(-17)) +  .9 (10^(-16)) (c^2);

Plot[𝒟2[c], {c, 0, 1}, 
        AxesLabel -> {"c", "𝒟2[c]"}];

we compute and define the derivative of 𝒟2:

(∂_c)((10^(-17)) +  .9 (10^(-16)) (c^2))

d𝒟2[c_] := 1.8 (10^(-16)) c

Our initial condition c[x, t = 0] and subsequent composition profiles c[x, t] will be discretized as well, ie., written as a list of 2 n + 1 points. The spatial difference or increment will thus be L/n

We now have to note that the finite difference scheme is approximate, and that the truncation errors might destabilize our solution or even preclude convergence to the real solution. It can be shown for the uniform 𝒟-case that  stability and convergence are assured provided that \(TraditionalForm\`Δt < Δx\^2\/\(2  𝒟\)\)

Without analyzing the full problem at hand, we take \(TraditionalForm\`Δt < Δx\^2\/\(2  𝒟\_max\) = \((L\_1/n)\)\^2\/\(2  𝒟2[1]\)\). For example \(TraditionalForm\`Δt = \(1\/2\) \((L\_1/n)\)\^2\/\(2  𝒟2[1]\)\)

c[x, t + Δt] /. 
            eq3 /. {t -> j Δt, 
            x -> i Δx} /. j -> j - 1 // 
      Simplify

With c[i, j] = 
        c[i Δx,  jΔt]:

Recursive option   (slow!!!)

ClearAll[c4]
\(c4[i_, j_, n_:  20, L_:  L\_1] := 0 /; \(\(i == \(-n\)\)\(
\)
        \( (*c4[i_, j_, n_:  20, L_:  L\_1] :=  .5 /; i == 0*) \)\)\)
c4[i_, j_, n_:  20, L_:  (Subscript[L, 1])] := 1 /; i == n
c4[i_, j_, n_:  20, L_:  (Subscript[L, 1])] := 
      0.5 +  .5 Sign[i] /; j == 0
Δt4[n_:  20, 
        L_:  (Subscript[L, 1])] := (((L/n))^2)/(4  𝒟2[1])

c4[i_, j_, n_:  20, 
        L_:  (Subscript[L, 1])] := 
c4[i, j, n, 
          L] = Δt (c4[i , j - 1 ]/Δt + (((c4[i - 1, j - 1] - 2 c4[i , j - 1] + c4[i + 1, j - 1]) 𝒟2[c4[i, j - 1]])/Δ)(x^2) + (((((c4[i - 1, j - 1] - c4[i + 1, j - 1]))^2) d𝒟2[c4[i, j - 1]])/(4 Δ(x^2)))) /. {Δx -> 
              L/n, Δt -> Δt4[n, 
                L]}

Note: the above definition is convenient in that it "remembers" previously calculated values. Let's now calculate the profile for a given time:

CL4[m_, n_:  20, L_:  (Subscript[L, 1])] := 
        Table[{1.  i L/n, c4[i, m, n, L]}, {i, -n, n}];

examples:

c4[1, 1]

CL4[1]

Unfortunately, caluculation proceeds very slowly. Let's try something else. We have an array cc of the form \(TraditionalForm\`{t, {c[x\_\(-n\), t],  ... , c[0, t], \(\(...\)\(.\)\), c[x\_n, t]}}\_t\)

Storing results in array

These are the x-coords of the points considered:

np = 20;   

    xx = Table[(Subscript[L, 1])/np  i, {i, -np, np}];

These is the concentration array, which contains only the initial conditions for now.

cc = {{0, Table[c4[i, 0, np, (Subscript[L, 1])], {i, -np, np}]}};

Now c4[i, j] is really cc[[j + 1, 2, 
        i + np + 1]]:

ClearAll[CL5];
\(CL5[j_, debug_:  True] := 
      Module[
{δt, δx}, 
δx = 
          L\_1/np; 
δt = Δt4[np, 
            L\_1]; 
If\ [
          j + 1 == Length[cc] + 1, 
If[debug, 
            Print["new calc"]]; 
AppendTo[
            cc, 
{j\ δt, 
Flatten[
{0, 
Table[
δt\ \((cc[[j, 2, 1 + i + np]]\/δt + \(d𝒟2[cc[[j, 2, 1 + i + np]]]\ \((cc[[j, 2, i + np]] - cc[[j, 2, 2 + i + np]])\)\^2\)\/\(4\ δx\^2\) + \(1\/δx\^2\) \((\((cc[[j, 2, i + np]] - 
                                    2\ cc[[j, 2, 
                                        1 + i + np]] + 
                                    cc[[j, 2, 
                                      2 + i + 
                                        np]])\)\ 𝒟2[cc[[j, 2, 
                                    1 + i + 
                                      np]]])\))\), 
{i, \(-np\) + 1, np - 1}
], 
                  1}
]}
]; 
          If[debug, 
            cc[[j + 1, 
              2]]]
, 
If[
            j + 1 > Length[cc] + 1, 
If[debug, 
              Print["Needs iteration ", j - 1]], 
If[
              debug, Print["old calc"]]; 
If[debug, 
              cc[[j + 1, 
                2]]]
]
]
]\)

CL5[1]

Now we plot different profiles! First, we compute the data:

For[i = 1, i ≤ 100, i++, CL5[i, False]];

Next we un-scale the x coordinates:

prof[n_] := 
        Transpose[{xx, 
            cc[[n + 1, 
              2]]}];
label[n_] := "t = " <> 
        ToString[
          cc[[n + 1, 
            1]]] <> " s"

MultipleListPlot[prof[0], prof[10], prof[20], prof[30], 
        PlotLegend -> {label[0], label[10], label[20], label[30], }, 
        PlotJoined -> True];

Note: in this approximation we set the concentration constant on the edges. This is not really correct and only valid for short times. The real boundary would be that the flux on x = -L and x = L is positive and negative respectively. The difference for long times is that the profile will become uniform, instead of the sloped line that we will eventually get.

Teaser

Select all the graphs and press <CTRL>+"y"

For[i = 100, i ≤ 200, i++, 
        CL5[i, False]];

mv = Table[ListPlot[prof[i]], {i, 0, 200}];