Diffusion in heterogeneous media¶
Diffusion in heterogeneous media normally implies a non-constant diffusion coefficient $\alpha = \alpha (x)$. A 1D diffusion model with such a variable diffusion coefficient reads
Written out, this becomes
where, e.g., an arithmetic mean can to be used for $\dfc_{i+\frac{1}{2}}$:
Implementation¶
Suitable code for solving the discrete equations is very similar to what we created for a constant $\dfc$. Since the Fourier number has no meaning for varying $\dfc$, we introduce a related parameter $D=\Delta t /\Delta x^2$.
def solver_theta(I, a, L, Nx, D, T, theta=0.5, u_L=1, u_R=0,
user_action=None):
x = linspace(0, L, Nx+1) # mesh points in space
dx = x[1] - x[0]
dt = D*dx**2
Nt = int(round(T/float(dt)))
t = linspace(0, T, Nt+1) # mesh points in time
u = zeros(Nx+1) # solution array at t[n+1]
u_n = zeros(Nx+1) # solution at t[n]
Dl = 0.5*D*theta
Dr = 0.5*D*(1-theta)
# Representation of sparse matrix and right-hand side
diagonal = zeros(Nx+1)
lower = zeros(Nx)
upper = zeros(Nx)
b = zeros(Nx+1)
# Precompute sparse matrix (scipy format)
diagonal[1:-1] = 1 + Dl*(a[2:] + 2*a[1:-1] + a[:-2])
lower[:-1] = -Dl*(a[1:-1] + a[:-2])
upper[1:] = -Dl*(a[2:] + a[1:-1])
# Insert boundary conditions
diagonal[0] = 1
upper[0] = 0
diagonal[Nx] = 1
lower[-1] = 0
A = scipy.sparse.diags(
diagonals=[diagonal, lower, upper],
offsets=[0, -1, 1],
shape=(Nx+1, Nx+1),
format='csr')
# Set initial condition
for i in range(0,Nx+1):
u_n[i] = I(x[i])
if user_action is not None:
user_action(u_n, x, t, 0)
# Time loop
for n in range(0, Nt):
b[1:-1] = u_n[1:-1] + Dr*(
(a[2:] + a[1:-1])*(u_n[2:] - u_n[1:-1]) -
(a[1:-1] + a[0:-2])*(u_n[1:-1] - u_n[:-2]))
# Boundary conditions
b[0] = u_L(t[n+1])
b[-1] = u_R(t[n+1])
# Solve
u[:] = scipy.sparse.linalg.spsolve(A, b)
if user_action is not None:
user_action(u, x, t, n+1)
# Switch variables before next step
u_n, u = u, u_n
The code is found in the file diffu1D_vc.py.
Stationary solution¶
As $t\rightarrow\infty$, the solution of the problem (1)-(4) will approach a stationary limit where $\partial u/\partial t=0$. The governing equation is then
with boundary conditions $u(0)=U_0$ and $u(L)=U_L$. mathcal{I}_t is possible to obtain an exact solution of (5) for any $\alpha$. Integrating twice and applying the boundary conditions to determine the integration constants gives
Piecewise constant medium¶
Consider a medium built of $M$ layers. The layer boundaries are denoted $b_0, \ldots, b_M$, where $b_0=0$ and $b_M=L$. If the layers potentially have different material properties, but these properties are constant within each layer, we can express $\alpha$ as a piecewise constant function according to
The exact solution (6) in case of such a piecewise constant $\alpha$ function is easy to derive. Assume that $x$ is in the $m$-th layer: $x\in [b_m, b_{m+1}]$. In the integral $\int_0^x (a(\xi))^{-1}d\xi$ we must integrate through the first $m-1$ layers and then add the contribution from the remaining part $x-b_m$ into the $m$-th layer:
Remark. mathcal{I}_t may sound strange to have a discontinuous $\alpha$ in a differential equation where one is to differentiate, but a discontinuous $\alpha$ is compensated by a discontinuous $u_x$ such that $\alpha u_x$ is continuous and therefore can be differentiated as $(\alpha u_x)_x$.
Implementation of diffusion in a piecewise constant medium¶
Programming with piecewise function definitions quickly becomes
cumbersome as the most naive approach is to test for which interval
$x$ lies, and then start evaluating a formula like
(8). In Python, vectorized expressions may
help to speed up the computations.
The convenience classes PiecewiseConstant and
IntegratedPiecewiseConstant in the Heaviside
module were made to simplify programming with
functions like (7) and expressions like
(8). These utilities not only represent
piecewise constant functions, but also smoothed versions of them
where the discontinuities can be smoothed out in a controlled fashion.
The PiecewiseConstant class is created by sending in the domain as a
2-tuple or 2-list and a data object describing the boundaries
$b_0,\ldots,b_M$ and the corresponding function values
$\alpha_0,\ldots,\alpha_{M-1}$. More precisely, data is a nested
list, where data[i][0] holds $b_i$ and data[i][1] holds the
corresponding value $\alpha_i$, for $i=0,\ldots,M-1$. Given $b_i$ and
$\alpha_i$ in arrays b and a, it is easy to fill out the nested
list data.
In our application, we want to represent $\alpha$ and $1/\alpha$ as piecewise constant functions, in addition to the $u(x)$ function which involves the integrals of $1/\alpha$. A class creating the functions we need and a method for evaluating $u$, can take the form
class SerialLayers:
"""
b: coordinates of boundaries of layers, b[0] is left boundary
and b[-1] is right boundary of the domain [0,L].
a: values of the functions in each layer (len(a) = len(b)-1).
U_0: u(x) value at left boundary x=0=b[0].
U_L: u(x) value at right boundary x=L=b[0].
"""
def __init__(self, a, b, U_0, U_L, eps=0):
self.a, self.b = np.asarray(a), np.asarray(b)
self.eps = eps # smoothing parameter for smoothed a
self.U_0, self.U_L = U_0, U_L
a_data = [[bi, ai] for bi, ai in zip(self.b, self.a)]
domain = [b[0], b[-1]]
self.a_func = PiecewiseConstant(domain, a_data, eps)
# inv_a = 1/a is needed in formulas
inv_a_data = [[bi, 1./ai] for bi, ai in zip(self.b, self.a)]
self.inv_a_func = \
PiecewiseConstant(domain, inv_a_data, eps)
self.integral_of_inv_a_func = \
IntegratedPiecewiseConstant(domain, inv_a_data, eps)
# Denominator in the exact formula is constant
self.inv_a_0L = self.integral_of_inv_a_func(b[-1])
def __call__(self, x):
solution = self.U_0 + (self.U_L-self.U_0)*\
self.integral_of_inv_a_func(x)/self.inv_a_0L
return solution
A visualization method is also convenient to have. Below we plot $u(x)$ along with $\alpha (x)$ (which works well as long as $\max \alpha(x)$ is of the same size as $\max u = \max(U_0,U_L)$).
%matplotlib inline
class SerialLayers:
...
def plot(self):
x, y_a = self.a_func.plot()
x = np.asarray(x); y_a = np.asarray(y_a)
y_u = self.u_exact(x)
import matplotlib.pyplot as plt
plt.figure()
plt.plot(x, y_u, 'b')
plt.hold('on') # Matlab style
plt.plot(x, y_a, 'r')
ymin = -0.1
ymax = 1.2*max(y_u.max(), y_a.max())
plt.axis([x[0], x[-1], ymin, ymax])
plt.legend(['solution $u$', 'coefficient $a$'], loc='upper left')
if self.eps > 0:
plt.title('Smoothing eps: %s' % self.eps)
plt.savefig('tmp.pdf')
plt.savefig('tmp.png')
plt.show()
Figure shows the case where
b = [0, 0.25, 0.5, 1] # material boundaries
a = [0.2, 0.4, 4] # material values
U_0 = 0.5; U_L = 5 # boundary conditions
Solution of the stationary diffusion equation corresponding to a piecewise constant diffusion coefficient.
By adding the eps parameter to the constructor of the SerialLayers
class, we can experiment with smoothed versions of $\alpha$ and see
the (small) impact on $u$. Figure
shows the result.
Solution of the stationary diffusion equation corresponding to a smoothed piecewise constant diffusion coefficient.
Axi-symmetric diffusion¶
Suppose we have a diffusion process taking place in a straight tube with radius $R$. We assume axi-symmetry such that $u$ is just a function of $r$ and $t$, with $r$ being the radial distance from the center axis of the tube to a point. With such axi-symmetry it is advantageous to introduce cylindrical coordinates $r$, $\theta$, and $z$, where $z$ is in the direction of the tube and $(r,\theta)$ are polar coordinates in a cross section. Axi-symmetry means that all quantities are independent of $\theta$. From the relations $x=\cos\theta$, $y=\sin\theta$, and $z=z$, between Cartesian and cylindrical coordinates, one can (with some effort) derive the diffusion equation in cylindrical coordinates, which with axi-symmetry takes the form
Let us assume that $u$ does not change along the tube axis so it suffices to compute variations in a cross section. Then $\partial u/\partial z = 0$ and we have a 1D diffusion equation in the radial coordinate $r$ and time $t$. In particular, we shall address the initial-boundary value problem
The condition (10) is a necessary symmetry condition at $r=0$, while (11) could be any Dirichlet or Neumann condition (or Robin condition in case of cooling or heating).
The finite difference approximation will need the discretized version of the PDE for $r=0$ (just as we use the PDE at the boundary when implementing Neumann conditions). However, discretizing the PDE at $r=0$ poses a problem because of the $1/r$ factor. We therefore need to work out the PDE for discretization at $r=0$ with care. Let us, for the case of constant $\dfc$, expand the spatial derivative term to
The last term faces a difficulty at $r=0$, since it becomes a $0/0$ expression caused by the symmetry condition at $r=0$. However, L'Hosptial's rule can be used:
The PDE at $r=0$ therefore becomes
For a variable coefficient $\dfc(r)$ the expanded spatial derivative term reads
We are interested in this expression for $r=0$. A necessary condition for $u$ to be axi-symmetric is that all input data, including $\alpha$, must also be axi-symmetric, implying that $\alpha'(0)=0$ (the second term vanishes anyway because of $r=0$). The limit of interest is
The PDE at $r=0$ now looks like
The fictitious value $u^n_{-1}$ can be eliminated using the discrete symmetry condition
which then gives the modified approximation to the term with the second-order derivative of $u$ in $r$ at $r=0$:
The discretization of the term with the second-order derivative in $r$ at any internal mesh point is straightforward:
To complete the discretization, we need a scheme in time, but that can be done as before and does not interfere with the discretization in space.
Spherically-symmetric diffusion¶
Discretization in spherical coordinates¶
Let us now pose the problem from the section Axi-symmetric diffusion in spherical coordinates, where $u$ only depends on the radial coordinate $r$ and time $t$. That is, we have spherical symmetry. For simplicity we restrict the diffusion coefficient $\dfc$ to be a constant. The PDE reads
for $r\in (0,R)$ and $t\in (0,T]$. The parameter $\gamma$ is 2 for spherically-symmetric problems and 1 for axi-symmetric problems. The boundary and initial conditions have the same mathematical form as in (9)-(12).
Since the PDE in spherical coordinates has the same form as the PDE in the section Axi-symmetric diffusion, just with the $\gamma$ parameter being different, we can use the same discretization approach. At the origin $r=0$ we get problems with the term
but L'Hosptial's rule shows that this term equals $\gamma\partial^2 u/ \partial r^2$, and the PDE at $r=0$ becomes
The associated discrete form is then
for a Crank-Nicolson scheme.
Discretization in Cartesian coordinates¶
The spherically-symmetric spatial derivative can be transformed to the Cartesian counterpart by introducing
Inserting $u=v/r$ in
yields
The two terms in the parenthesis can be combined to
The PDE for $v$ takes the form
For $\alpha$ constant we immediately realize that we can reuse a solver in Cartesian coordinates to compute $v$. With variable $\alpha$, a "reaction" term $v/r$ needs to be added to the Cartesian solver. The boundary condition $\partial u/\partial r=0$ at $r=0$, implied by symmetry, forces $v(0,t)=0$, because