Math Forum Discussions

Torsten

Posts: 1,129
Registered: 11/8/10

Re: PDE unstable, need help
Posted: Feb 6, 2013 10:42 AM

Plain Text

"Moritz" wrote in message <ketp47$qmn$1@newscl01ah.mathworks.com>...
> 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

pdepe is not suited to solve your problem.
This software is:

http://www.google.de/url?sa=t&rct=j&q=&esrc=s&frm=1&source=web&cd=1&ved=0CDIQFjAA&url=http%3A%2F%2Fwww.nag.co.uk%2Fnumeric%2FMB%2Fmanual_21_1%2Fpdf%2FD03%2Fd03ps.pdf&ei=-HgSUYzFC4rEsgaopIHABg&usg=AFQjCNEENu2rvqxsmFpasSbcLofNKyQkFA&sig2=9JoMeLt7yylOsPFHffNBTw

Best wishes
Torsten.