Question 7.PS.1: Unsteady Diffusion Through a Membrane Consider the membrane ...
Unsteady Diffusion Through a Membrane
Consider the membrane in Fig. 7.4 separating two well-mixed fluids, each containing solute A at a fixed concentration. On the left side of the membrane the solute has a concentration c_i, while on the right side the concentration is c_0. Within the membrane itself, the initial concentration is also c_0. We want to study how the concentration changes as time progresses, and solute diffuses from the left side of the membrane to the right. The diffusion equation for a binary, dilute system with a constant diffusion coefficient is
\frac{∂c}{∂t} = D \frac{∂^2c}{∂x^2} (7.3.1)
At this point you might notice that the formulation of this problem is identical to determining the temperature profile due to conduction through a solid wall, or to determining the velocity profile for Couette flow between two flat plates. Although the physical mechanisms governing these three problems are different, the mathematics describing the processes are the same.
It is convenient before going further to put this problem in dimensionless form. If we choose θ = (c − c_o)/(c_i − c_o), ξ = x/L, and τ = tD/L^2, then we get
\frac{∂θ}{∂τ} = \frac{∂^2θ}{∂ξ^ 2} (7.3.2)
with
θ(ξ , 0) = 0, θ(0, τ ) = 1, θ(1, τ ) = 0. (7.3.3)
We are going to consider two cases. In the first case, both reservoirs are of infinite size. As a result, transport of solute from one reservoir to the other does not change the boundary concentration. This is a solvable problem via separation of variables. We just need to be a bit careful about making the boundary conditions homogeneous. To solve the problem, we define another function
y = θ − (1 − ξ ) (7.3.4)
which is just subtracting off the steady state solution for θ. The transformed variable also obeys the diffusion equation
\frac{∂y}{∂τ} = \frac{∂^2y}{∂ξ^ 2} (7.3.5)
with the homogeneous boundary conditions y(0, τ ) = 0 and y(1, τ ) = 0. The initial condition is now different,
y(ξ , 0) = −1 + ξ (7.3.6)
The differential equation for y is solved by separation of variables to give
y =\sum\limits_{n=1}^{\infty }{c_n sin(nπξ )e^{−(nπ)^2τ }} (7.3.7)
which you can readily verify satisfies the homogeneous Dirichlet boundary conditions at ξ = 0 and ξ = 1. From the initial condition, we have
ξ − 1 = \sum\limits_{n=1}^{\infty }{c_n sin(nπξ )} (7.3.8)
We can compute the coefficients
c_n =\frac{ 2}{nπ} (7.3.9)
following the approach we used in Section 1.2.6, integrating by parts to handle the ξ term on the left-hand side. The solution for the concentration is then
θ = 1 − ξ −\sum\limits_{n=1}^{\infty }{\frac{2}{nπ} sin(nπξ )e^{−(nπ)^2τ}} (7.3.10)
We will use this solution to test out the accuracy of the numerical solution. In the second case, we will consider a finite size reservoir on the right boundary.
In this case, the transport from the left boundary will lead to an increase in the concentration on the right boundary. The steady state for this finite reservoir problem is θ = 1, since eventually the reservoir on the right will be filled up from diffusion from the left one.
Infinite Reservoir
Let’s begin with the infinite reservoir case, using a second-order centered finite difference approximation for the spatial derivative, and forward Euler integration for the time derivative. We will use n equally spaced nodes. For an interior node i, we get
\frac{dθ_i}{dτ} =\frac{ θ^{(k)}_{i+1} − 2θ^{(k)}_i + θ^{(k)}_{i−1}}{Δx^2} (7.3.11)
Solving the system of ODEs using explicit Euler leads to
θ^{(k+1)}_i = θ^{(k)}_i +\frac{h}{Δx^2} θ^{(k)}_{i+1} − 2θ^{(k)}_i + θ^{(k)}_{i−1} (7.3.12)
The left and right nodes are unchanged in time, so their explicit Euler form becomes
θ^{(k+1)}_1 = θ^{(k)}_1 (7.3.13)
and
θ^{(k+1)}_n = θ^{(k)}_n (7.3.14)
This is the problem we solved in Example 7.2, just with different boundary conditions and initial conditions. As a result, we don’t need to write any new MATLAB programs yet.
Figure 7.5 compares the result from explicit Euler to the exact solution in Eq. (7.5).For the finite differences solution, we used 51 nodes and a time step size of 10^{−4}. For very short times, the finite difference solution seems to be more accurate than evaluating the Fourier series because we need many terms in the Fourier series to accurately capture the boundary layer near the left reservoir. (There are more sophisticated analytical techniques, known as perturbation methods, that circumvent this shortcoming with Fourier series solutions.) The numerical solution does not oscillate because it is not actually a sum of sine waves. Moreover, there is an interesting issue related to the rate of propagation of information in finite differences. Recall that the concentrations at each node are connected to the previous node, so it takes a finite number of time steps before the concentration at a node can propagate to a distant node. While it seems like the numerical integration is superior to the exact solution, remember that neither approach is perfect.
The numerical solution has multiple sources of error, including the truncation errors of finite differences and explicit Euler. If we look at the two solutions for a longer time point, we see that the agreement is quite good. The Fourier series becomes accurate because we need only a few modes to capture most of the features of the solution; the exponential terms in Eq. (7.3.10) decay quickly with n when τ is large. So long as we use a sufficient number of nodes and a small enough time step, the truncation errors in the numerical solution are also small with respect to the resolution of the plot.
While Eq. (7.3.10) is an exact result, it requires an infinite number of terms until it gives a result with infinite precision. Obviously, when we produced Fig. 7.5, we had to truncate the sum at a finite number of terms. To see the relative difference between the two solutions, Fig. 7.6 presents the quantity ||θ_{analytical} − θ_{numerical}|| as a function of the K terms that we included to compute the sum. The result is interesting. The plot of the error versus K steadily decreases for small K, but levels out for large K. This is because, when we reach large enough K, the exponential term in Eq. (7.3.10) becomes very small. As a result, when we add another term, we are below machine precision and nothing happens. The total error is then given by the error in the numerical solution. Note that if you use a smaller grid and time step size, then the plateau is reached at a lower value of error and a higher value of K.
Finite Reservoir
Now that we understand the infinite reservoir case, let’s consider the case where the reservoir on the right has a finite size L_T . Returning briefly to dimensional variables, the mass flux at the right boundary of the membrane is
j = −D\frac{∂c}{∂x}|_{x=L} (7.3.15)
Over a time step of duration Δt, we then expect the amount of mass inside the tank to increase by an amount
Δm = jAΔt (7.3.16)
where A is the area of the membrane. To convert this into concentration in the tank, c_T , we should divide by the volume V = AL_T of the tank,
Δc_T = j \frac{Δt}{L_T} (7.3.17)
If we assume that the tank is well mixed and that there is no resistance to mass transfer to the tank (i.e., the partition coefficient is unity), then the concentration of the tank is equal to the concentration c(1, t). With the finite reservoir model, all that has changed in our system of ODEs is the concentration at node n. Putting everything together, using an explicit Euler approximation for all of the changes in time, and a backward difference equation for the flux, we have
c^{(k+1)}_ n − c^{(k)}_n = −D(\frac{c^{(k)}_{n} − c^{(k)}_{n−1}}{Δx}) \frac{Δt}{L_T} (7.3.18)
Converting to dimensionless notation gives
θ^{(k+1)}_n = θ^{(k)}_n − \epsilon \frac{h}{Δx}\left(θ^{(k)}_n − θ^{(k)}_{n−1}\right) (7.3.19)
where \epsilon = L/L_T is the size of the membrane relative to the size of the tank on the right. This is generally a small number since we think of the membrane as being thin.
To solve the finite tank problem, we just need to make small changes to our original programs. First, we need to alter the explicit Euler program to take the value of \epsilon:
1 function [xout,yout] = ...
ivp_explicit_euler_sys(x0,y0,h,n_out,i_out,eps)
2 % x0 = initial value for x
3 % y0 = initial value for y
4 % h = step size
5 % i_out = frequency to output x and y to [xout, yout]
6 % n_out = number of output points
7 % requires the function getf(x,y)
8 xout = zeros(n_out+1,1); xout(1) = x0;
9 yout = zeros(n_out+1,length(y0)); yout(1,:) = y0;
10 x = x0; y = y0;
11 for j = 2:n_out + 1
12 for i = 1:i_out
13 y=y+h*getf(y,eps);
14 x = x + h;
15 end
16 xout(j,1) = x;
17 yout(j,:) = y;
18 end
We then need to modify the function getf to account for the finite reservoir at node n:
1 function f = getf(y,eps)
2 n = length(y);
3 dx = 1/(n-1);
4 f = zeros(n,1);
5 f(1) = 0;
6 for i = 2:n-1
7 f(i) = y(i+1) - 2*y(i) + y(i-1);
8 end
9 f = f/dx^2;
10 f(n) = -eps/dx*(y(n)-y(n-1));
The changes are in lines 9 and 10; instead of setting f(n) = 0, we now use Eq. (7.3.19). For simplicity, we moved the division by Δx^2 before we set the value of node n so that line 10 in our program corresponds exactly to the notation in Eq. (7.3.19).
Figure 7.7 shows how the finite size of the tank affects the solution. As we see in Fig. 7.7a, the concentration profile inside the membrane looks the same as that for an infinite reservoir at short times. This makes sense, since for short times there is no accumulation in the right tank. However, as we see in Fig. 7.7b, we start to see substantial accumulation in the right reservoir as material finally reaches the tank. This raises the concentration of the right boundary of the membrane, reducing the driving force for mass transfer. As a result, the mass transfer slows down as a function of time.
