用 Matlab 求解偏微分方程

Question

`sol = pdepe(m,@ParticleDiffusionpde,@ParticleDiffusionic,@ParticleDiffusionbc,x,t);
% Extract the first solution component as u.
u = sol(:,:,:); 

function [c,f,s] = ParticleDiffusionpde(x,t,u,DuDx)
global Ds 
c = 1/Ds;
f = DuDx;
s = 0;

function u0 = ParticleDiffusionic(x)
global qo
u0 = qo;

function [pl,ql,pr,qr] = ParticleDiffusionbc(xl,ul,xr,ur,t,x)

global Ds K n
global Amo Gc kf rhop
global uavg
global dr R nr

sum = 0;
for i = 1:1:nr-1
r1 = (i-1)*dr; % radius at i
r2 = i * dr; % radius at i+1
r1 = double(r1); % convert to double precision
r2 = double(r2);
sum = sum + (dr / 2 * (r1*ul+ r2*ur));
end;
uavg = 3/R^3 * sum;

ql = 1;
pl = 0;

qr = 1;
pr = -((kf/(Ds.*rhop)).*(Amo - Gc.*uavg - ((double(ur/K)).^2).^(n/2) ));`


dq(r,t)/dt = Ds( d2q(r,t)/dr2 + (2/r)*dq(r,t)/dr )

q(r, t=0) = 0

dq(r=0, t)/dr = 0

dq(r=dp/2, t)/dr = (kf/Ds*rhop) [C(t) - Cp(at r = dp/2)]

q = solid phase concentration of trace compound in a particle with radius dp/2

C = bulk liquid concentration of trace compound

Cp = trace compound concentration at particle surface

我想用给定的初始和边界条件来解决上述 pde。尝试了 Matlab 的 pdepe，但效果不佳。也许边界条件给我带来了问题。我还使用这个等温方程来平衡：q = K*Cp^(1/n)。这是对流扩散方程，但我找不到任何可以正确解决此类方程的文章。

Answer 1

目前的实现有两个问题。

不正确的来源术语

您尝试求解的 PDE 具有以下形式

具有等价形式

最后一项是由于原始 PDE 中的因子 2 而出现的。 最后一个词需要通过源词并入pdepe。

计算q 平均值

当前实现尝试使用传递给边界条件函数的q 的左右值来计算q 的平均值。 这是不正确的。 q 的平均值需要根据数量的最新值向量来计算。

然而，我们有一个复杂的情况是接收所有网格值的唯一函数是ParticleDiffusionpde；但是，传递给该函数的网格值不能保证来自我们提供的网格。

解决方案：使用事件（如pdepe 文档中所述）。 这是一个 hack，因为事件函数旨在检测过零，但它的优点是该函数在我们提供的网格上被赋予了 q 的所有值。

所以，下面的工作示例（您会注意到我将所有参数设置为 1，因为我不知道更好）使用 events 函数来更新可以通过以下方式访问的变量 qStore边界条件函数（see here 解释），边界条件函数执行向量化梯形积分以进行平均计算。

工作示例

function [] = ParticleDiffusion()

    %   Parameters
    Ds   = 1;
    q0   = 0;
    K    = 1;
    n    = 1;
    Amo  = 1;
    Gc   = 1;
    kf   = 1;
    rhop = 1;

    %   Space
    rMesh = linspace(0,1,10);
    rMesh = rMesh(:) ;
    dr    = rMesh(2) - rMesh(1) ; 

    %   Time
    tSpan = linspace(0,1,10);

    %   Vector to store current u-value
    qStore = zeros(size(rMesh));

    options.Events = @(m,t,x,y) events(m,t,x,y);

    %   Solve
    [sol,~,~,~,~] = pdepe(1,@ParticleDiffusionpde,@ParticleDiffusionic,@ParticleDiffusionbc,rMesh,tSpan,options);



    %   Use the events function to update qStore
    function [value,isterminal,direction] = events(m,~,~,y)
        qStore     = y; % Value of q on rMesh
        value      = m; % Since m is constant, it will never be zero (no event detection)
        isterminal = 0; % Continue integration
        direction  = 0; % Detect all zero crossings (not important)
    end


    function [c,f,s] = ParticleDiffusionpde(r,~,~,DqDr)

        %   Define the capacity, flux, and source
        c = 1/Ds;
        f = DqDr;
        s = DqDr./r;

    end

    function u0 = ParticleDiffusionic(~)
        u0 = q0;
    end

    function [pl,ql,pr,qr] = ParticleDiffusionbc(~,~,R,ur,~)

        %   Calculate average value of current solution
        qL    = qStore(1:end-1);
        qR    = qStore(2: end );
        total = sum((qL.*rMesh(1:end-1) + qR.*rMesh(2:end))) * dr/2;
        qavg  = 3/R^3 * total;

        %   Left boundary
        pl = 0;
        ql = 1;

        %   Right boundary
        qr = 1;
        pr = -(kf/(Ds.*rhop)).*(Amo - Gc.*qavg - (ur/K).^n);
    end

end