Thread Subject:

PDE unstable, need help

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

"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.

"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 !