How to understand the factors (c,f,s) in pdepe and how to formulate boundary conditions?

Hello!

Maybe some general words to set the ground. The system I am trying the simulate is a 1D porous material that is passed through by some air flow.

I took over the program some time ago and was so far rather using it more than "really understanding" it.

The system of partial differential equations looks like this:

1. A_G du(1)/dt = -m_GS_dot(u(1),u(2),u(3),u(4)) - m_G_dot * du(1)/dx
2. A_S du(2)/dt = m_GS_dot(u(1),u(2),u(3),u(4))
3. A_G * C_G(u(1)) * du(3)/dt = q_GS_dot(u(3),u(4)) - m_G_dot * C_G(u(1)) * du(3)/dx
4. A_S * C_S(u(2)) * du(4)/dt = - q_GS_dot(u(3),u(4)) + m_GS_dot(u(1),u(2),u(3),u(4)) * h(u(2),u(4))

(N.B. I couldn't find the appropriate symbol for \partial but the derivatives expressed by usual upright 'd's are meant to be partial derivatives.)

The program is running and everything seems ok. I just don't understand a couple of decisions that were made when programming the code.

Question 1:

In the actual definition of the pde where the matrices c, f, and s are defined the components of c are straight forward (the factors in front of the time derivatives on the lhs of the equations). I would have naively thought that the spatial derivatives present in eqns. (I) and (III) would go to f, representing the flux term where I thought that convective and diffusive contributions would be considered. However, they were put in s together with all the rest that does in fact represent a source while f was set equal to zero.

To clarify:

c1 = A_G;
c2 = A_S;
c3 = A_G * C_G(u1);
c4 = A_S * C_S(u2);

c = [c1; c2; c3; c4];

f1 = 0;
f2 = 0;
f3 = 0;
f4 = 0;

f = [f1; f2; f3; f4];

s1 = -m_GS_dot(u(1),...,u(4)) - cnst.m_G_dot_s * dudx(1);
s2 = m_GS_dot(u(1),...,u(4));
s3 = q_GS_dot(u(2),u(4)) - m_G_dot * C_G(u(1)) * dudx(3);
s4 = m_GS_dot(u(1),...,u(4)) * h(u(2),u(4)) - q_GS_dot(u(2),u(4));

s = [s1; s2; s3; s4];

So far so good. At this point it seems like no big deal whether to express the spatial derivative with s, which should be allowed since s = s(x,t,u,dudx) or with f. Still when checking other solvers the terms of a conservation equation are separated into a transient term, a convective, a diffusive and a source term. Matlab only uses the expression flux and source apart from the obvious transient term, where I would have counted our convective term as part of a flux as opposed to a source but that might only be semantics...?

Question 2:

When it comes to the formulation of the boundary conditions, the program code completely confuses me.

u(1)_in = 0.001;
u(2)_in = 0.05;
u(3)_in = 290;
u(4)_in = u(3)_in;

pl = [ul(1) - u(1)_in;...
      ul(2) - u(2)_in;...
      ul(3) - u(3)_in;...
      ul(4) - u(4)_in];
ql = [0;0;0;0];
pr = [0;0;0;0];
qr = [1;1;1;1];
% qr = [5;5;5;5];

With eqn. (1-6) from the pdepe documentation this is a little weird at least for me. While the definition of pl still seems ok, forcing identity of the solution variables and the user defined values u(i)_in on the left boundary, I have no idea why they even bothered setting ql = [0;0;0;0] with f == 0! The same applies for the right side where it gets even a little weirder since the calculation won't run if qr was set to [0;0;0;0]. You can see by the commented line that apparantly there was some playing around with the value of qr. Again, with f == 0 I don't see what difference this makes. Finally, I also don't understand the boundary condition pr = [0;0;0;0].

I hope the questions are posed clearly enough. If not, I will amend that any time.

Thanks for your help!

Cheers, Markus

[Torsten]

The main point is that - in principle - pdepe is not suited to solve your system of PDEs (1)-(4). The "pe" in pde"pe" stands for "parabolic-elliptic" ; your PDEs are hyperbolic in nature (there is no spatial second derivative term present). So you have to find a way out that for your system, no boundary condition can be imposed at x=xr. If pdepe works with your boundary condition (0+1*0=0), then be happy and don't care about the "why". Unfortunately, there is no hyperbolic solver in matlab directly that solves your problem, but if you have access to the NAG Toolbox,try

http://www.nag.co.uk/numeric/MB/manual64_24_1/html/D03/d03pff.html

Best wishes Torsten.

[Markus]

Hi Torsten!

Thank you for your answer. However, I think I don't fully understand it. I had to check again for the distinctions parabolic, elliptic, and hyperbolic and found that it is made by a relation between the coefficients in front of the highest derivatives. I am sure you know that but just to make sure that I am thinking in the right direction, let's say we have

A * d^2 u(x,t)/dx^2 + B * d^2 u(x,t)/dtdx + C * d^2 u(x,t)/dt^2 + ... = 0

then

4AC - B^2 < 0 would correspond to a hyperbolic equation, 4AC - B^2 = 0 to a parabolic equation, and 4AC - B^2 > 0 to an elliptic equation

For me this is so far mostly nomenclature that doesn't tell me much but as you said at least it would give a hint about whether the solver can be used or not.

Since I have no second derivative w.r.t. either one variable nor the mixed derivative w.r.t. both variables I would have thought that in my case A=B=C=0 and that the equation classifies as parabolic.

Where's the error?

Can you also comment on the formulation of the pde? Would it make any difference or say would it be possible to formulate the equations without setting f = [0;0;0;0]?

So far I do not have access to the NAG toolbox. I was looking into using a solver library written in Python as an alternative. But before doing so I wanted to understand how and especially why the code I have is "working".

Best regards, Markus

[Torsten]

For your question about the classification of your PDE-system:

Your equations (1) and (3) are one-way wave equations (with source terms) which are the prototype for all hyperbolic partial differential equations.

For your question about the formulation of the pde:

First think about which artificial boundary condition you want to apply at x=xr for u(1) and u(3) to make pdepe work (as I already mentionned above, from the physical point of view, no boundary condition is allowed to be set at x=xr because the solution is only influenced from the interior of the domain, not from exterior conditions). If this boundary condition can be implemented within pdepe, you get your f automatically (but I doubt that any senseful boundary condition can be deduced that can be implemented in pdepe).

If you don't want to include the convection term in MATLAB's s, from the equations, the only choice for f would be f=[u(1);0;u(3);0].

Now the boundary conditions for u(1) and u(3) read:

p+q*f=0,

and the only senseful boundary condition for u(1) and u(3) at x=xr are

du(1)/dx = du(3)/dx =0.

But how to choose pr(1), qr(1), pr(3) and qr(3) to get du(1)/dx = du(3)/dx =0 at x=xr if f(1)=u(1) and f(3)=u(3) ? No way.

Moreover, your equations (2) and (4) are usual ODEs without spatial derivatives - so no boundary conditions at x=xl and x=xr are allowed to be prescribed for them at all. This again is not possible within pdepe.

To summarize:

From the theoretical point of view: The use of pdepe for your problem is not justified. Use an adequate tool to solve your equations - e.g. d03pff from the NAG library.

Form the practical point of view: If your settings work to solve the problem, why changing the solver ?

Best wishes

Torsten.

[Markus]

Hi Torsten!

Thanks for the clarification on the terminology. I used your hints to do a more profound search on the web and found many references to advection problems with and without source terms.

I also checked out the NAG toolbox and saw that they offer a 30 day trial period. I was thinking that I might first catch up on some basic understanding and then giving the toolbox a try before asking our IT department whether we could licence the package.

Since you seem very knowledgeable on that topic and potentially also on the usage of the NAG toolbox it would be great if there was some way to reach you in a PM.

Is there?

Best regards,

Markus

[Torsten]

Sure.

My Email-address is

[e-mail]

Best wishes

Torsten.