Circumvent NDSolve::bdord: error for 1-D Euler Equations

Question

I attempted to use NDSolve for the 1-D isentropic unsteady flow equations with low subsonic inflow velocity and prescribed inflow total enthalpy; along with a non-reflective subsonic outflow boundary condition (BC) based on the method of characteristics. The latter BC is commonly used in Computational Fluid Dynamics (CFD) and in Computational Aero-Acoustics. NDSolve rejects this BC outright in Mathematica 9.0.1.

Code:

Auxiliary Eq's and parameters.

gamma = 1.4;
sigma = 0.25;
uInfty = 0.02;
M$ref = uInfty/Sqrt[gamma - 1];
dhInfty = (-(1/2)) uInfty^2;
h[x_, t_] := 1 + dh[x, t];
c[x_, t_] := ((gamma - 1)*h[x, t])^(1/2);
k = sigma*(1 - M$ref^2);
rule = {
 residu -> D[dh[x, t] + (1/2)*u[x, t]^2, x],
 residh -> u[x, t]*D[dh[x, t], x] + (gamma - 1)*D[u[x, t], x]};

PDE's for dependent variables u and dh vs. x and t.

AllEqs = {
 ueqn -> D[u[x, t], t] + residu == 0,
 dheqn -> D[dh[x, t], t] + residh == 0} /. rule;

BC and IC.

BCICRules = {
 bcInFlow -> {(*Inflow Boundary x=0 (subsonic)*)
  u[0, t] == uInfty,
  dh[0, t] + (1/2)*u[0, t]^2 == 0},
 bcOutFlow -> {(*Outflow Boundary x=1 (subsonic, non-reflective)*)
  Derivative[0, 1][dh][1, t] - c[1, t]*Derivative[0, 1][u][1, t] + 
   k*c[1, t]*(dh[1, t] - dhInfty) == 0},
 ic -> {(*IC*)
   u[x, 0] == (1 - x)*uInfty,
   dh[x, 0] == (1 - x)*dhInfty}
};

Error:

Solve the equations numerically to steady state

NDSolve[Flatten@({ueqn, dheqn, bcInFlow, bcOutFlow, ic} /. AllEqs /. BCICRules), 
{dh[x, t], u[x, t]}, {t, 0, 3.}, {x, 0, 1.}]

Following error occurs!

NDSolve::bdord: Boundary condition 0.157956 (0.0002 +dh[1,t]) Sqrt[1+dh[1,t]]+(dh^(0,1))[1,t]-0.632456 Sqrt[1+dh[1,t]] (u^(0,1))[1,t] should have derivatives of order lower than the differential order of the partial differential equation. >>

Observation:

• The reason for failure is bcOutFlow specification. Changing it to the following by vanishing the time derivatives results into a solution without any complain.

New bcOutFlow.

bcOutFlow -> {
  Derivative[0, 0][dh][1, t] - c[1, t]*Derivative[0, 0][u][1, t] + 
  k*c[1, t]*(dh[1, t] - dhInfty) == 0}

Here is the plot.

Plot3D[Evaluate[{#[x, t]} /. soln0], {t, 0, 3}, {x, 0, 1}, 
Mesh -> None, ColorFunction -> "DarkRainbow", 
PlotStyle -> Opacity[.7], PlotPoints -> 40, ImageSize -> 400, 
PlotLabel -> #] & /@ {dh, u} // Row

• Above error states that BC is not allowed to contain time derivatives of the same order as the time derivatives in the PDE's (namely, first-order time derivatives). THIS CLAIM IS PATENTLY INCORRECT FOR THE COMPRESSIBLE FLOW EQUATIONS. The non-reflective boundary condition containing time derivatives is widely used, as anyone can verify from the CFD literature of the past two decades.
• I had hoped that, within Mathematica's Method of Lines, the non-reflective outflow boundary condition should just give rise to another ODE in time at the outflow-boundary grid point $x=1$. That ODE should be coupled automatically to those at the interior grid points through the stencil of the built-in default spatial differencing scheme.
• The pre-test has caused Mathematica to give the NDSolve::bdord: error and reject the problem outright. If the pre-test were bypassed, then there might be hope that NDSolve would translate the time-derivative outflow boundary condition into a time-dependent ODE for the dependent variables at the boundary grid point $x=1$. It could then go forward with the numerical integration of the ODE system.

Question

1. Does anyone out there know of a way to circumvent this difficulty with the NDSolve implementation in Mathematica 9.0.1?
2. Specifically, does anyone know if it is possible to bypass NDSolve's pre-test for derivatives in the boundary conditions?
3. Any suggestion how to accommodate reflective BC for a PDE like the above one in Mathematica?

Answer 1

I managed to solve your equation set by

1. adding a very small viscosity μ to the equation to circumvent the computational difficulty caused by the discontinuities of spatial derivative i.e. D[u[x, t], x] and D[dh[x, t], x] in the solution.

2. discretizing the PDEs to a set of ODEs to circumvent the bdord warning, via pdetoode.

The following is the code:

gamma = 14/10; sigma = 25/100; uInfty = 2/100;
M$ref = uInfty/Sqrt[gamma - 1];
dhInfty = (-(1/2)) uInfty^2;
h[x_, t_] = 1 + dh[x, t];
c[x_, t_] = ((gamma - 1) h[x, t])^(1/2);
k = sigma (1 - M$ref^2);
residu = D[dh[x, t] + (1/2) u[x, t]^2, x]; 
residh = u[x, t] D[dh[x, t], x] + (gamma - 1) D[u[x, t], x];

With[{μ = 5 10^-4}, {ueqn, dheqn} = {D[u[x, t], t] + residu == μ D[u[x, t], x, x], 
    D[dh[x, t], t] + residh == 0}];
bcInFlow = {u[0, t] == uInfty, dh[0, t] + (1/2) u[0, t]^2 == 0}; 
bcOutFlow = {Derivative[0, 1][dh][1, t] - c[1, t] Derivative[0, 1][u][1, t] + 
    k c[1, t] (dh[1, t] - dhInfty) == 0};
ic = {u[x, 0] == (1 - x) uInfty, dh[x, 0] == (1 - x) dhInfty};

domain = {L, R} = {0, 1}; tend = 3;
points = 500; xdifforder = 8;
grid = Array[# &, points, domain];
var = {u /@ grid, dh /@ grid};
(* Definition of pdetoode isn't included in this code piece,
   please find it in the link above. *)
ptoo = pdetoode[{u, dh}[x, t], t, grid, xdifforder];
odeu = Delete[#, {{1}, {-1}}] &@ptoo@ueqn;
odedh = Delete[#, {{1}(*, {-1}*)}] &@ptoo@dheqn;
odeic = ptoo@ic;
{odebcin, odebcout} = ptoo@{bcInFlow, bcOutFlow};
sollst = NDSolveValue[{odeu, odedh, odeic, odebcin, odebcout}, var, {t, 0, tend}, 
   MaxSteps -> Infinity, SolveDelayed -> True];
{solu, soldh} = rebuild[#, grid, -1] & /@ sollst;
Manipulate[GraphicsRow[Plot[#[x, t], {x, L, R}, 
     PlotRange -> #2] & @@@ {{solu, {0, 0.03}}, {soldh, {0, 0.016}}}, 
  ImageSize -> 800], {t, 0, tend}]

Remark

1. In principle SolveDelayed -> True is equivalent to Method -> {"EquationSimplification" -> "Solve"}, and the latter is recommended, but actually, when added as a counter measure to the warning ndsdtc, the latter will make NDSolve return unevaluated, at least in v9.0.1, I think it's a bug.

2. The definition of odeu and odedh can also be odeu = Delete[#, {{1}(*,{-1}*)}] &@ptoo@ueqn; odedh = Delete[#, {{1}, {-1}}] &@ptoo@dheqn;, but this definition takes more time to solve.

3. It seems that points should be large enough, and xdifforder should be a large enough even number.

4. μ should be small enough, but not too small, or NDSolve may spit out the mconly warning and fails.