From: Esther Rojas on
Thank you, Torsten

Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote in message <1807341531.54202.1280919593899.JavaMail.root(a)gallium.mathforum.org>...
> > 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.

I agree with you. In my opinion, the main problem with dT/dz is due to the initial conditions I've imposed: a hard step. This step, since the convection term is so important, moves along the time, keeping the stepness it has from the very beggining. At the temperature front -as you named it before- dT/Dz does not properly exists, so the pedepe function has it hard to solve. If a softer initial step is imposed -using a tanh function, for example-, the no-physical oscillations disappear.
>
> 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)).)
>

I'll try the backward differences to check if there was any difference

> 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.
You are rigth, nothing changes (I've tryed the above formulation some time ago, finding no mayor differences)


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

Your comments are always helpful. Thank you again,

Esther
>
> Best wishes
> Torsten.
From: Torsten Hennig on
> Thank you, Torsten
>
> Torsten Hennig <Torsten.Hennig(a)umsicht.fhg.de> wrote
> in message
> <1807341531.54202.1280919593899.JavaMail.root(a)gallium.
> mathforum.org>...
> > > 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.
>
> I agree with you. In my opinion, the main problem
> with dT/dz is due to the initial conditions I've
> imposed: a hard step. This step, since the convection
> term is so important, moves along the time, keeping
> the stepness it has from the very beggining. At the
> temperature front -as you named it before- dT/Dz does
> not properly exists, so the pedepe function has it
> hard to solve. If a softer initial step is imposed
> -using a tanh function, for example-, the no-physical
> oscillations disappear.
> >
> > 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)).)
> >
>
> I'll try the backward differences to check if there
> was any difference
>

I checked it within my solver: the oscillations
disappeared.
But the above backward difference scheme is only first-order accurate and instead of oscillations it introduces
artificial diffusion to the solution - the solution
curves become smeared (similarily as if you choose
a higher value for the thermal conductivity).
To get both - an oscillation-free _and_ exact solution,
the upwind scheme should be at least second-order
accurate.
Such a scheme is described in the following article:
http://portal.acm.org/citation.cfm?id=155272

But I think the easiest solution is to choose
a very high number of discretization points in pdepe -:)

Best wishes
Torsten.

> > 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.
> You are rigth, nothing changes (I've tryed the above
> formulation some time ago, finding no mayor
> differences)
>
>
> >
> > 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.
>
> Your comments are always helpful. Thank you again,
>
> Esther
> >
> > Best wishes
> > Torsten.