Unstable solution with pdepe
5 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
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
1 Commento
Youssef Khmou
il 6 Feb 2013
hi, just edit your code to appear clearly by moving the line "function Direct" forward with space tab .
Risposte (0)
Vedere anche
Categorie
Scopri di più su Numerical Integration and Differential Equations in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!