Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: PDE unstable, need help
Date: Wed, 6 Feb 2013 15:54:08 +0000 (UTC)
Organization: Universit&#228;t f&#252;r Bodenkultur Wien
Lines: 92
Message-ID: <ketub0$j2f$1@newscl01ah.mathworks.com>
References: <ketp47$qmn$1@newscl01ah.mathworks.com> <kettkg$gcq$1@newscl01ah.mathworks.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: www-03-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: newscl01ah.mathworks.com 1360166048 19535 172.30.248.48 (6 Feb 2013 15:54:08 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Wed, 6 Feb 2013 15:54:08 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 3192309
Xref: news.mathworks.com comp.soft-sys.matlab:788405

"Torsten" wrote in message <kettkg$gcq$1@newscl01ah.mathworks.com>...
> "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.

THANK YOU !