Prev: Trajectory in VRML
Next: close('name') or close(handle)
From: Esther Rojas on 3 Aug 2010 09:47
Having to solve the following equation.
dT/dt=9.7e-7d2T/d2z-6.4e-3dT/dz-3.3e-7T
in z€[0, 6], t€[0, 725] and T€[0.72,1].
I manage to run it with pdepe, using:
Sol0=pdepe(0,@pd1D1SFpdev2,@pd1D1SFicv2,@pd1D1SFbcv2,z,t), where

function [c,f,s]=pd1D1SFpdev2(z,t,u,DuDz)
c=1;
f=9.74e-7*DuDz-6.43e-3*u;
s=-3.28E-07*u;

function [p1 q1 pr qr]=pd1D1SFbcv2(z1,u1,zr,ur,t)
p1=u1-1;q1=0;
qr=1;pr=6.43e-3*ur;

function u0 = pd1D1SFicv2(z)
if z==0, u0=1; else u0=0.72; end

The problem is that the discretization in z-coordenate I have to use to have a physical meaningfull solution is too big. If I use a refiner zmesh, the solutions are oscillating. According to the help, the pdepe is the one who chooses which time-step is the best. I've tried reducing RelTol, AbsolTol without any improvement. Switching NormControl to 'on' improves a little the solution, but not enough to avoid osillations. The same happens reduccing Initial Step and maxStep.

Any idea where to look at?. Thank you in advance,
Esther
From: Torsten Hennig on 3 Aug 2010 06:53
> Having to solve the following equation.
> dT/dt=9.7e-7d2T/d2z-6.4e-3dT/dz-3.3e-7T
> in z€[0, 6], t€[0, 725] and
> T€[0.72,1].
> I manage to run it with pdepe, using:
> Sol0=pdepe(0,@pd1D1SFpdev2,@pd1D1SFicv2,@pd1D1SFbcv2,z
> ,t), where
>
> function [c,f,s]=pd1D1SFpdev2(z,t,u,DuDz)
> c=1;
> f=9.74e-7*DuDz-6.43e-3*u;
> s=-3.28E-07*u;
>
> function [p1 q1 pr qr]=pd1D1SFbcv2(z1,u1,zr,ur,t)
> p1=u1-1;q1=0;
> qr=1;pr=6.43e-3*ur;
>
> function u0 = pd1D1SFicv2(z)
> if z==0, u0=1; else u0=0.72; end
>
> The problem is that the discretization in
> z-coordenate I have to use to have a physical
> meaningfull solution is too big. If I use a refiner
> zmesh, the solutions are oscillating. According to
> the help, the pdepe is the one who chooses which
> time-step is the best. I've tried reducing RelTol,
> AbsolTol without any improvement. Switching
> NormControl to 'on' improves a little the solution,
> but not enough to avoid osillations. The same happens
> reduccing Initial Step and maxStep.
>
> Any idea where to look at?. Thank you in advance,
> Esther

Two factors make your problem difficult to solve
for a parabolic-elliptic solver like pdepe:

1. The discontinuity of the solution at the moving
(temperature ?) front
2. The convection dominance of the problem (i.e.
the term v*dT/dz strongly dominates the term
D*d^2T/dz^2).

The only remedy is to _increase_ the number of
discretization points.
For your parameter settings, I get oscillation-free solutions for a z-grid of about 3000 uniformly
distributed points.

Best wishes
Torsten.
From: Torsten Hennig on 4 Aug 2010 00:29
> > Having to solve the following equation.
> > dT/dt=9.7e-7d2T/d2z-6.4e-3dT/dz-3.3e-7T
> > in z€[0, 6], t€[0, 725] and
> > T€[0.72,1].
> > I manage to run it with pdepe, using:
> >
> Sol0=pdepe(0,@pd1D1SFpdev2,@pd1D1SFicv2,@pd1D1SFbcv2,z
>
> > ,t), where
> >
> > function [c,f,s]=pd1D1SFpdev2(z,t,u,DuDz)
> > c=1;
> > f=9.74e-7*DuDz-6.43e-3*u;
> > s=-3.28E-07*u;
> >
> > function [p1 q1 pr qr]=pd1D1SFbcv2(z1,u1,zr,ur,t)
> > p1=u1-1;q1=0;
> > qr=1;pr=6.43e-3*ur;
> >
> > function u0 = pd1D1SFicv2(z)
> > if z==0, u0=1; else u0=0.72; end
> >
> > The problem is that the discretization in
> > z-coordenate I have to use to have a physical
> > meaningfull solution is too big. If I use a
> refiner
> > zmesh, the solutions are oscillating. According to
> > the help, the pdepe is the one who chooses which
> > time-step is the best. I've tried reducing RelTol,
> > AbsolTol without any improvement. Switching
> > NormControl to 'on' improves a little the
> solution,
> > but not enough to avoid osillations. The same
> happens
> > reduccing Initial Step and maxStep.
> >
> > Any idea where to look at?. Thank you in advance,
> > Esther
>
> Two factors make your problem difficult to solve
> for a parabolic-elliptic solver like pdepe:
>
> 1. The discontinuity of the solution at the moving
> (temperature ?) front
> 2. The convection dominance of the problem (i.e.
> the term v*dT/dz strongly dominates the term
> D*d^2T/dz^2).
>
> The only remedy is to _increase_ the number of
> discretization points.
> For your parameter settings, I get oscillation-free
> solutions for a z-grid of about 3000 uniformly
> distributed points.
>

5000 discretization points give a perfect solution.

> Best wishes
> Torsten.

Best wishes
Torsten.
From: Esther Rojas on 4 Aug 2010 05:56
Thank you Torsten for your support (I really appreciate very much you are there).
have you imposed any special requirements when using 5000 discretization points?. I've run the programme with such a fine mesh but I still keep on having some fluctuations: no higher that 1.0112 for Nz=5000 and not close to the temperature step, but they are there.
Is there any possibility to impose that dT/dz is always < or = to 0?. In that case, oscillations would not be allowed from the very beggining, since we know how the temperature profile looks like.
thank you again,

Esther


Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <1234978681.53535.1280910605986.JavaMail.root(a)gallium.mathforum.org>...
> > > Having to solve the following equation.
> > > dT/dt=9.7e-7d2T/d2z-6.4e-3dT/dz-3.3e-7T
> > > in z&#8364;[0, 6], t&#8364;[0, 725] and
> > > T&#8364;[0.72,1].
> > > I manage to run it with pdepe, using:
> > >
> > Sol0=pdepe(0,@pd1D1SFpdev2,@pd1D1SFicv2,@pd1D1SFbcv2,z
> >
> > > ,t), where
> > >
> > > function [c,f,s]=pd1D1SFpdev2(z,t,u,DuDz)
> > > c=1;
> > > f=9.74e-7*DuDz-6.43e-3*u;
> > > s=-3.28E-07*u;
> > >
> > > function [p1 q1 pr qr]=pd1D1SFbcv2(z1,u1,zr,ur,t)
> > > p1=u1-1;q1=0;
> > > qr=1;pr=6.43e-3*ur;
> > >
> > > function u0 = pd1D1SFicv2(z)
> > > if z==0, u0=1; else u0=0.72; end
> > >
> > > The problem is that the discretization in
> > > z-coordenate I have to use to have a physical
> > > meaningfull solution is too big. If I use a
> > refiner
> > > zmesh, the solutions are oscillating. According to
> > > the help, the pdepe is the one who chooses which
> > > time-step is the best. I've tried reducing RelTol,
> > > AbsolTol without any improvement. Switching
> > > NormControl to 'on' improves a little the
> > solution,
> > > but not enough to avoid osillations. The same
> > happens
> > > reduccing Initial Step and maxStep.
> > >
> > > Any idea where to look at?. Thank you in advance,
> > > Esther
> >
> > Two factors make your problem difficult to solve
> > for a parabolic-elliptic solver like pdepe:
> >
> > 1. The discontinuity of the solution at the moving
> > (temperature ?) front
> > 2. The convection dominance of the problem (i.e.
> > the term v*dT/dz strongly dominates the term
> > D*d^2T/dz^2).
> >
> > The only remedy is to _increase_ the number of
> > discretization points.
> > For your parameter settings, I get oscillation-free
> > solutions for a z-grid of about 3000 uniformly
> > distributed points.
> >
>
> 5000 discretization points give a perfect solution.
>
> > Best wishes
> > Torsten.
>
> Best wishes
> Torsten.
From: Torsten Hennig on 4 Aug 2010 02:59
> Thank you Torsten for your support (I really
> appreciate very much you are there).
> have you imposed any special requirements when using
> 5000 discretization points?. I've run the programme
> with such a fine mesh but I still keep on having some
> fluctuations: no higher that 1.0112 for Nz=5000 and
> not close to the temperature step, but they are
> there.
> Is there any possibility to impose that dT/dz is
> always < or = to 0?. In that case, oscillations would
> not be allowed from the very beggining, since we know
> how the temperature profile looks like.
> thank you again,
>
> Esther
>

I don't have MATLAB installed, but I did put the
options in my own solver similar to those of pdepe -
so slight differences in the solution behaviour
might occur.

The reason for the oscillations lies in the
discretization of the first-order derivative
of T with respect to z.

pdepe uses centred differences:
dT/dz|z=z_i ~ (T(z_(i+1))-T(z_(i-1)))/(z_(i+1)-z_(i-1))
But it's well-known that for convection-dominated
flows upwind schemes have to be used to avoid
oscillations (the simplest one for your equation
would be
dT/dz|z=z_i ~ (T(z_(i))-T(z_(i-1)))/(z_(i)-z_(i-1)).)

Unfortunately you can't choose the discretization scheme
in pdepe.

The only remedies are:
1. Increase the number of discretization points
2. Increase the diffusion coefficient artificially
(this won't change the solution very much, but will
help to cure numerical instabilities ; e.g. try
9.74e-5 d^2T/dz^2 instead of 9.74e-7 d^2T/dz^2)

Maybe you can try the setting

function [c,f,s]=pd1D1SFpdev2(z,t,u,DuDz)
c=1;
f=9.74e-7*DuDz;
s=-3.28E-07*u-6.43e-3*DuDz;

function [p1 q1 pr qr]=pd1D1SFbcv2(z1,u1,zr,ur,t)
p1=u1-1;q1=0;
qr=1;pr=0;

function u0 = pd1D1SFicv2(z)
if z==0, u0=1; else u0=0.72; end

but I think this won't change anything.

As far as I know, there is no 'simple' MATLAB
routine you can call and which has upwind schemes
implemented.

If you have access to the NAG numerical library:
there are routines especially suited for
convection-dominated flows.
If you are interested, I can search for their names.

Sorry, not much help.

Best wishes
Torsten.
First  |  Prev  |  Next  |  Last
Pages: 1 2 3
Prev: Trajectory in VRML
Next: close('name') or close(handle)