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

4 views (last 30 days)
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
  5 Comments
Markus
Markus on 27 Oct 2014
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

Sign in to comment.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!