Unstable solution with pdepe

Hi,

i tried to implement one part of a sedimentation model in Matlab. Because my mathematical background is not so strong ( i am working on it) i tried to use pdepe.

What you will see is a 3D plot of a sedimentation process. A growing sediment layer and a clear liquor interface (cli) (i think you would call that a shock and a fan ?).

The code is instable at the clear liquor interface and when the sediment and the cli meet.

Any suggestions ? Limiter functions ? (aside of the probably bad programming style)

My code:

function Direct %#ok<*NASGU> global u0 omega sigmaN k drho c g wun r0 rb xmax uc %initial Parameters from S.Berres Chem Eng J 111 (2004) 91-103 % Attempt to solve the second order parabolic equation from the phenomenological % sedimentation model. If any of the authors may read this code, yes i am % considering to use your method. % % % % % c=5; uc=0.08; sigmaN=5.7; k=9; xmax=0.66; % u0=0.1; % should be >uc drho=1660; % density difference in kg/m^3 wun=0.001; % sedimentation velocity of a single floc r0=0.06; % radius at the liquor height rb=0.3; % radius at the bottom of the tube g=9.81; % earth acceleration m/s^2 omega=sqrt(100/rb); % % m=0; Nx=200; conc0=0.07; Nt=50; tend=10; xmesh=linspace(r0,rb,Nx); tmesh=linspace(0,tend,Nt); % sol=pdepe(m,@pdefun,@icfun,@bcfun,xmesh,tmesh); Concentration=sol(:,:,1); figure surf(xmesh,tmesh,Concentration) xlabel('radius') ylabel('time') % function [a,f,s]=pdefun(xmesh,t,u,ux) %global xc omega sigma0 k drho c g wun r0 rb % if (u > uc & u < xmax) fbk=wun.*u.*(1-u).^c; sigmaE=sigmaN.*((u/uc)^k-1); sigmaEu=sigmaN.*k.*(u/uc).^(k-1)./uc; kla=-fbk.*sigmaEu./(drho.*g.*u); % f=kla.*ux+omega.^2.*xmesh./g*fbk; f = -kla.*ux-omega^2.*xmesh.*fbk; s=0; a=1; end

    function uinit = icfun(r)
        uinit = u0;%(xmesh==r);
        %u0=u0;
    end

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

      %global xc omega sigma0 k drho c g wun r0 rb
        pl=0;
        ql=1;
        pr=0;
        qr=1;
    end
end

[Youssef Khmou]

hi, just edit your code to appear clearly by moving the line "function Direct" forward with space tab .