# Transport Equation in a Tube: Source Term on Boundary

I'm modeling mass transport in a flow reactor. The flow reactor is a tube, which allows me to use cylindrical symmetry in solving the Convection-Diffusion-Reaction (CDR) Equation, which governs the mass transport.

The equation that I'm solving is:

$D ( \frac{\partial^2 u}{\partial r^2} + \frac{1}{r}\frac{\partial u}{\partial r} ) - v(r) \frac{\partial u}{\partial z} = S(r)$

Here, r and z are the coordinates in the respective directions, u is the property I want to solve for (the mass concentration), v is the velocity field, S is a given source and D is the diffusion coefficient.

The Boundary conditions are:

$\frac{\partial c}{\partial r} = 0,$ at $r=0$ and at $r=R$ (where R is the length of the rectangle in the r direction).

The system is sketched in the figure below:

My approach is to use the Method of Lines, i.e. to solve the system in the r direction for each step in the z direction.

I discretised the above equation as follows:

$\frac{du}{dz} = \frac{D}{v_i}(\frac{u_{i+1} - 2u_i + u{i-1}}{2h}+\frac{1}{r_i}\frac{u_{i+1}-u_{i-1}}{2h}) - S_i$

For both boundary conditions I used the central difference approximation, i.e. at r = 0 and r = R:

$\frac{u_{i+1} - u_{i-1}}{2h} = 0$

I am using a grid with two ghost points, as shown in the figure below:

Now, the code works ok, if I do not put the source on any of the boundaries, see the first plot below. However, in the actual experiment the source is in the centre of the flow reactor, i.e. the source is located on the boundary r = 0. When I put the source there, I get quite large oscillations (see the second plot below).

I think that I've either implemented the wrong boundary condition or that I'm using the wrong grid for this kind of computation.

My thoughts on removing those oscillations are: 1. use a staggered grid, i.e. do not actually put a source on the boundary 2. use a mixed boundary condition

What would you recommend? Am I missing something else?

The MatLab code is also attached below.

clear all
close all

%radius of flow reactor in mm
R = 13

%step-size in r direction
hr = 0.2
r = -hr : hr : R + hr

%parameters
D = 1.0 %diffusion
v = 1./(13.1 - r)
v = ones(1,length(r))

%build diagonals a, b, c of matrix A
%
%      |b c 0 0 0|
%      |a b c 0 0|
%  A = |0 a b c 0|
%      |0 0 a b c|
%      |0 0 0 a b|
%
% with a = alpha - beta, b = -2alpha, c = alpha + beta

alpha = ones(length(r)-2,1) ./ v(2:end-1)' .* (D/hr^2)
beta = ones(length(r)-2,1) ./ (v(2:end-1).* r(2:end-1))' .* (D/(2*hr))

diag_a = alpha - beta
diag_b = -2*alpha
diag_c = alpha + beta

A= spdiags([diag_a diag_b diag_c],[-1 0 1],length(r)-2,length(r)-2)

%implement boundary conditions du/dr = 0 at r = 0
A(1,1) = -4*alpha(1)
A(2,1) =  4*alpha(1)

%implement boundary conditions du/dr = 0 at r = R
A(end-1,end)= A(end-1,end)*2

ht = 0.00001
u = zeros(length(r)-2,1);

source = zeros(length(r)-2,1);
source(1:length(r)/2-2) = 1.

n_steps = 500

clear ustored

ustored = []

%step forward in time with explicit Euler
for i = 1:n_steps

ustored = [ustored u];
u_next = u(:,1) + ht*(A*u(:,1)+source(:,1));
u = u_next;

end

ustored = [ustored u]

time = 0:ht:(n_steps)*ht;

mesh(time,r(2:end-1),ustored)

• What happens if you use a finer grid? Or if you locate the source at the other boundary? – Jan Jan 31 '14 at 13:05
• @Jan I tried to use both a finer grid and a higher time step, a.k.a. z step, the oscillations remain. – seb Jan 31 '14 at 15:07
• What happens if you move your innermost point from r=0 to a slightly larger value? Maybe the singularity at r=0, is causing problems? (And it relation to that, how does your ghost point behave there?) – Bernhard Jan 31 '14 at 15:12
• @Bernhard thanks, I will try that out. Although I used L'Hopital's rule to take care of the 1/r * du/dr singularity at r = 0, i.e. I approximated it to be d^2/d^r2 c at this particular grid point. – seb Jan 31 '14 at 15:32