我正在求解一个在出口处具有 Neumann 边界条件且在入口处具有恒定浓度 C 的一维扩散方程。最后，我想观察浓度如何沿 x 方向随时间扩散。

$$\frac{\partial C}{\partial t} = D\frac{\partial ^2 C}{\partial x^2}$$

下面是用 MATLAB 编写的代码。

function main()
C0 = [500 480 460 440 420 400 380 360 340];
tspan = [0 10];
[t C]  = ode23s(@(t,s) fun(t,s), tspan , C0);
plot(t,C)

function dC = fun(t,C)
D  = 3000;
retain_rows = 2:8;
tri = -(D/(1/5^2))*full(gallery('tridiag',9,-1,2,-1)); % 5 is del x
tri = tri(retain_rows,:);
dC(1,1) = 0; % constant concentration at inlet
dC(retain_rows,1) = tri*C;
dC(9,1) = -D*((C(9)-C(8))/2); Neumann boundary condition at outlet
end
end

我有从 ode 求解器获得的解决方案。我想知道如何从获得的溶液中检查摩尔平衡。另外，有人可以解释如何更改上述代码以在入口处应用恒定通量吗？

有什么建议？