我正在模拟流动反应器中的质量传递。流动反应器是一个管子，它允许我使用圆柱对称性来求解控制质量传递的对流-扩散-反应 (CDR) 方程。

我正在求解的方程是：

$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)$

这里，r 和 z 是各自方向上的坐标，u 是我要求解的属性（质量浓度），v 是速度场，S 是给定源，D 是扩散系数。

边界条件为：

$\frac{\partial c}{\partial r} = 0,$在和处（其中 R 是 r 方向上矩形的长度）。$r=0$$r=R$

系统示意图如下图：

我的方法是使用直线法，即在 z 方向上的每一步都在 r 方向上求解系统。

我将上面的方程离散化如下：

$\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$

对于这两个边界条件，我使用了中心差分近似，即在 r = 0 和 r = R：

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

我使用的是带有两个鬼点的网格，如下图所示：

现在，代码可以正常工作了，如果我没有将源代码放在任何边界上，请参见下面的第一个图。然而，在实际实验中，源位于流动反应器的中心，即源位于边界 r = 0 上。当我将源放在那里时，我得到了相当大的振荡（见下面的第二张图）。

我认为我要么实现了错误的边界条件，要么为这种计算使用了错误的网格。

我对消除这些振荡的想法是： 1. 使用交错网格，即实际上不要将源放在边界上 2. 使用混合边界条件

你会推荐什么？我还缺少其他东西吗？

MatLab 代码也附在下面。

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)