Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

enter image description here

  • 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 MA 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?
share|improve this question
4  
I find two problems with your post. First, your code is unreadable. Please format it in the normal Mathematica input form. Second, it isn't clear what kind of help you are expecting. –  m_goldberg May 1 '13 at 3:15
2  
You should copy the InputForm version of your code to make it better readable. The warning message you get from Mathematica is probably not about what is theoretically possible, but about BC conditions that Mathematica can handle. –  Sjoerd C. de Vries May 1 '13 at 5:22
    
From the Author: I improved code readability by converting Mathematica StandardForm commands to InputForm. I also tried to be more specific about what kind of help I hope someone can give. Would you please consider re-opening the Question? I hope my editing will answer the objections posed by commentors @m_goldberg and Sjoerd C. de Vries –  P. D. Thomas May 1 '13 at 20:05

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.