The Kitchin Research Group

We want to solve for the concentration profile of component that diffuses into a 1D rod, with an impermeable barrier at the end. The PDE governing this situation is:

\(\frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2}\)

at \(t=0\), in this example we have \(C_0(x) = 0\) as an initial condition, with boundary conditions \(C(0,t)=0.1\) and \(\partial C/ \partial x(L,t)=0\).

We are going to discretize this equation in both time and space to arrive at the solution. We will let \(i\) be the index for the spatial discretization, and \(j\) be the index for the temporal discretization. The discretization looks like this.

Note that we cannot use the method of lines as we did before because we have the derivative-based boundary condition at one of the boundaries.

We approximate the time derivative as:

\(\frac{\partial C}{\partial t} \bigg| _{i,j} \approx \frac{C_{i,j+1} - C_{i,j}}{\Delta t} \)

\(\frac{\partial^2 C}{\partial x^2} \bigg| _{i,j} \approx \frac{C_{i+1,j} - 2 C_{i,j} + C_{i-1,j}}{h^2} \)

We define \(\alpha = \frac{D \Delta t}{h^2}\), and from these two approximations and the PDE, we solve for the unknown solution at a later time step as:

\(C_{i, j+1} = \alpha C_{i+1,j} + (1 - 2 \alpha) C_{i,j} + \alpha C_{i-1,j} \)

We know \(C_{i,j=0}\) from the initial conditions, so we simply need to iterate to evaluate \(C_{i,j}\), which is the solution at each time step.

See also: http://www3.nd.edu/~jjwteach/441/PdfNotes/lecture16.pdf

import numpy as np
import matplotlib.pyplot as plt

N = 20  # number of points to discretize
L = 1.0
X = np.linspace(0, L, N) # position along the rod
h = L / (N - 1)          # discretization spacing

C0t = 0.1  # concentration at x = 0
D = 0.02

tfinal = 50.0
Ntsteps = 1000
dt = tfinal / (Ntsteps - 1)
t = np.linspace(0, tfinal, Ntsteps)

alpha = D * dt / h**2
print alpha

C_xt = [] # container for all the time steps

# initial condition at t = 0
C = np.zeros(X.shape)
C[0] = C0t

C_xt += [C]

for j in range(1, Ntsteps):
    N = np.zeros(C.shape)
    N[0] =  C0t
    N[1:-1] = alpha*C[2:] + (1 - 2 * alpha) * C[1:-1] + alpha * C[0:-2]
    N[-1] = N[-2]  # derivative boundary condition flux = 0
    C[:] = N
    C_xt += [N]

    # plot selective solutions
    if j in [1,2,5,10,20,50,100,200,500]:
        plt.plot(X, N, label='t={0:1.2f}'.format(t[j]))

plt.xlabel('Position in rod')
plt.ylabel('Concentration')
plt.title('Concentration at different times')
plt.legend(loc='best')
plt.savefig('images/transient-diffusion-temporal-dependence.png')

C_xt = np.array(C_xt)
plt.figure()
plt.plot(t, C_xt[:,5], label='x={0:1.2f}'.format(X[5]))
plt.plot(t, C_xt[:,10], label='x={0:1.2f}'.format(X[10]))
plt.plot(t, C_xt[:,15], label='x={0:1.2f}'.format(X[15]))
plt.plot(t, C_xt[:,19], label='x={0:1.2f}'.format(X[19]))
plt.legend(loc='best')
plt.xlabel('Time')
plt.ylabel('Concentration')
plt.savefig('images/transient-diffusion-position-dependence.png')

plt.show()

0.361361361361

The solution is somewhat sensitive to the choices of time step and spatial discretization. If you make the time step too big, the method is not stable, and large oscillations may occur.

Copyright (C) 2013 by John Kitchin. See the License for information about copying.

org-mode source