用 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

用 Matlab 求解偏微分方程

Solving PDE with Matlab

matlab

pde

源词不正确

计算`q`平均值

工作示例

用 Matlab 求解偏微分方程

Solving PDE with Matlab

matlab

pde

源词不正确

计算q平均值

工作示例

计算`q`平均值