输运方程的解析解

计算科学 pde matlab 有限差分 数值建模 平流扩散
2021-11-28 06:45:02

我在看对流扩散方程的解析解

Ct=D2Cx2vCx
初始条件为 和狄利克雷边界条件
c(x,0)=4000

C(x=0,t>0)=4100

纽曼边界条件

Cx=0 at x=L;t>0.

在此处输入图像描述

但是,当对上述解析解进行编码时,我会得到 C 的负值,这是不现实的。

function AnalyticalSoln2()
format long
R = 1;
D = 900;
v = 200;
L = 60; 
co = 4100;
ci = 4000;
X = linspace(0,60,10);
t = 0:0.001:2;
sol=[];
for pos = 1:length(X)
    x = X(pos);
A1 = 0.5*erfc((x.*R-t.*v)./(2*(t.*D*R).^0.5));
A2 = 0.5*exp(x.*v/D).*erfc((x.*R+t.*v)./(2*(D*t.*R).^0.5));
A31 = 0.5*(2+v*(2*L-x)/D + (t.*v^2)/(D*R))*exp(v*L/D);
A32 = erfc((R*(2*L-x)+t.*v)./(2*(D*t.*R).^0.5));
A41 = -((t.*v^2)./(pi*D*R)).^0.5;
A42 = exp((v*L/D)-(R./(4*t.*D)).*(2*L -x + t.*v/R).^2);
A = A1 + A2 + A31.*A32 + A41.*A42;
if t==0 
C = ci + (co - ci)*A'
else
C = ci + (co - ci)*A' - co*A';
end
sol = horzcat(sol,C);
end
sol(1:100:end,:)
end

而从 pdepe 求解器获得的数值解是非负的。

这是使用 pdepe 获得的解决方案

function DiffusionConvectionMATLAB
format short
global D
m = 0;
x = linspace(0,60,10);
t = 0.1:0.1:2;
D = 900;
sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t)
    function [g,f,s] = pdefun(x,t,c,DcDx)
    v = 200;
    g = 1;
    f = D*DcDx;
    s = -v*DcDx;
    end

    function c0 = icfun(x)
    c0 = 4000;
    end

    function [pl,ql,pr,qr] = bcfun(xl,cl,xr,cr,t)
    pl = cl-4100;
    ql = 0;
    pr = 0;
    qr = 1;
    end
end

有人可以建议我从分析表达式计算解决方案的方式是否有任何问题?

1个回答

这是实现分析解决方案的一种方法。t0 = 0但是,如果我使用(正如您在评论中提到的),我似乎没有得到正确的情节。如果您设置t0的时间大于您拥有的最后时间(t0 = 10例如),它会起作用。

我还使用了函数句柄,以便A(x,t)可以在任何值上进行x评估t

format long
R = 1;
D = 900;
v = 200;
L = 60; 
co = 4100;
ci = 4000;
x = linspace(0,60,100);
tvec = 0.0:0.01:2;
sol=[];

%Use function handles so that A(x,t) can be evaluated for any x and t
A1 = @(x,t) 0.5*erfc((x.*R-t.*v)./(2*(t.*D*R).^0.5));
A2 = @(x,t) 0.5*exp(x.*v/D).*erfc((x.*R+t.*v)./(2*(D*t.*R).^0.5));
A31 = @(x,t) 0.5*(2+v*(2*L-x)/D + (t.*v^2)/(D*R))*exp(v*L/D);
A32 = @(x,t) erfc((R*(2*L-x)+t.*v)./(2*(D*t.*R).^0.5));
A41 = @(x,t) -((t.*v^2)./(pi*D*R)).^0.5;
A42 = @(x,t) exp((v*L/D)-(R./(4*t.*D)).*(2*L -x + t.*v/R).^2);
A = @(x,t) A1(x,t) + A2(x,t) + A31(x,t).*A32(x,t) + A41(x,t).*A42(x,t);

t0 = 0; %this doesn't seem to give the right solution
%t0 = 10; %this seems to work

for t = tvec %loop over time
    if (t <= t0)
        C = ci + (co - ci)*A(x,t);
    else
        C = ci + (co - ci)*A(x,t) - co*A(x,t-t0);
    end

    sol = [sol; C];
end