*** Transient diffusion - partial differential equations
:PROPERTIES:
:categories: PDE
:tags: mass transfer
:date: 2013/04/02 09:27:42
:updated: 2013/04/02 09:44:05
:END:
index:pde
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.
[[./images/pde-diffusion-discretization-scheme.png]]
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
#+BEGIN_SRC python
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()
#+END_SRC
#+RESULTS:
: 0.361361361361
[[./images/transient-diffusion-temporal-dependence.png]]
[[./images/transient-diffusion-position-dependence.png]]
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.