Solving a system of partial differential and algebraic equations? setting the step-size to get reasonable results

Question

I'm trying to use Mathematica to get some early approximate solutions to a system of algebraic and partial differential equations. It is actually 1D model of an ideal gas in a tube. I'll divide my question into two parts so it will be more understandable:

• Part one:

I want to solve Navier–Stokes equations for a compressible inviscid isothermal fluid in a 1D tube:

1. $\frac{\partial P}{\partial t}+\frac{\partial \left( \nu P \right)}{\partial u}=0$

2. $\frac{\partial \left(P \nu\right) }{\partial t}+\frac{\partial \left( P \nu^2 \right) }{\partial u}+\mathring{R} T_a\frac{\partial P}{\partial u}=0$

Initial condition:

1. $\nu\left( u,0 \right)=0$
2. $P\left( u,0 \right)=1$

Dirichlet boundary condition (step signal)

$\left\{\begin{matrix} if \, t<1 \, then & \nu\left( 0,t \right)=0 \\ else & P\left( 0,t \right)=2 \end{matrix}\right.$

and $P\left( 1,t \right)=1$

I know I can use NDSolveValue and NDSolve to have some numerical solutions. for example I used the command below to solve a simplified version of the above system of PDEs:

pval = NDSolveValue[{D[p[x, t]*v[x, t], x] + D[p[x, t], t] == 0, 
D[p[x, t], x] + D[p[x, t]*v[x, t]^2, x] + D[p[x, t]*v[x, t], t] == 
0, p[0, t] == UnitStep[t - 1] + 1, p[1, t] == 1, p[x, 0] == 1, 
v[x, 0] == 0}, p, {x, 0, 1}, {t, 0, 3}, 
Method -> {"MethodOfLines", "TemporalVariable" -> t, 
"SpatialDiscretization" -> {"TensorProductGrid", 
  "MaxPoints" -> 4000, "MinPoints" -> 4000}}, MaxStepSize -> 0.001]
Dimensions[sol1["Grid"]]

but the result I get is pretty awful:

ContourPlot[pval[x, t], {x, 0, 1}, {t, 0, 2.2}]

and it does not go furthure than t=2.232...

My first guess I don't have the right temporal resolution, but I don't know how to set it like what I have found for spatial resolution.

But I would like to ask my question in a general form:

1. How should I solve this system of algebraic and PDEs. DSolve, NDSolveValue or NDSolve?

2. how to get reasonable solutions by setting the spatial and temporal step size?

   • Part two:

I would like to solve the Navier–Stokes equations for a compressible viscose adiabatic fluid in a 1D tube:

1. Ideal gas: $P=\mathring{R} T \rho$
2. conservation of mass: $\frac{\partial \rho}{\partial t}+\frac{\partial \left( \nu \rho \right)}{\partial u}=0$
3. conservation of linear momentum: $\frac{\partial \left(\rho \nu\right) }{\partial t}+\frac{\partial \left( \rho \nu^2 \right) }{\partial u}+\frac{\partial P}{\partial u}-\frac{4\mu\nu}{D}=0$
4. conservation of heat for a control volume, assuming the Newtonian viscose fluid, the conduction/radiation is negligible (adiabatic for the whole tube): $\mu\nu^2=\frac{cD}{4}\left( \frac{\partial\left( T\rho\nu \right)}{\partial u}+\frac{\partial \left(\rho T \right)}{\partial t} \right)$

these equations in Mathematica format:

p[x, t] == T[x, t]*q[x, t]
D[q[x, t]*v[x, t], x] + D[q[x, t], t] == 0
D[p[x, t], x] + D[q[x, t]*v[x, t]^2, x] + D[q[x, t]*v[x, t], t] -v[x,t]==0
v[x, t]^2 - D[T[x, t]*q[x, t]*v[x, t], x] == 0

Where $\mathring{R}$, $\nu$, $D$ and $c$ are constants. $P$, $\rho$, $\nu$ and $T$ are functions of $x$ and $t$

I would like to ask my question in a general form:

1. How should I solve this system of algebraic and PDEs. DSolve, NDSolveValue or NDSolve?
2. how to get reasonable solutions by setting the spatial and temporal step size?

P.S. I'm not able to get the right tags for my post!

Answer 1

Solution of "Simplified Version" with Original BC

The following modifications to the "simplified version" of this hydrodynamics problem with the boundary condition, p[0, t] == 0.1*t + 1 (now replaced by the OP), allows the computation to be carried out to t == 10, but not much beyond.

{pval, vval} = NDSolveValue[{D[p[x, t]*v[x, t], x] + D[p[x, t], t] == 0 , 
   D[p[x, t], x] + D[p[x, t]*v[x, t]^2, x] + D[p[x, t]*v[x, t], t] == 0, 
   p[0, t] == 0.1*t + 1, p[1, t] == 1, p[x, 0] == 1, v[x, 0] == 0}, 
   {p, v}, {x, 0, 1}, {t, 0, 10}, 
   Method -> {"MethodOfLines", "TemporalVariable" -> t, "SpatialDiscretization" -> 
       {"TensorProductGrid", "MaxPoints" -> 1000, "MinPoints" -> 1000}}, 
   MaxStepSize -> {Automatic, 5 10^-4}]

NDSolveValue::eerr: Warning: scaled local spatial error estimate of 168.55610168929198at t = 10. in the direction of independent variable x is much greater than the prescribed error tolerance. Grid spacing with 1000 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.

In fact, reducing the spatial and temporal step sizes does not help much, and I suspect that some form of artificial viscosity may be needed. Here are plots of the results. Spurious fine-scale behavior is evident in the first plot.

ContourPlot[pval[x, t], {x, 0, 1}, {t, 0, 10}, PlotRange -> All]

ContourPlot[vval[x, t], {x, 0, 1}, {t, 0, 10}, PlotRange -> All]

Note that a similar problem is solved in 11748, although I am a bit suspicious of the two approaches employed.

Solution of "Simplified Version" with Present BC

The desired boundary condition formally is If[t < 1, v[0, t] == 0, p[0, t] == 2]. However, I do not believe that NDSolve can handle this boundary condition for nonlinear PDEs. There is, however, a straightforward equivalent approach. The initial and boundary conditions for t < 1 yield the solution {pval[x, t] == 1, vval[x, t] ==0}, both analytically and numerically. So, beginning the calculation with p[0, t] == 2 might seem like the natural course of action. Unfortunately, NDSolve (correctly) complains about inconsistency between the initial and boundary conditions for p at {x == 0, t == 0} and replaces the boundary condition by p[0, t] == 2 - Exp[1 - t]], after which it runs more or less as desired. Specifically, with the boundary conditions just described explicitly included, the code and results are

{pval, vval} = NDSolveValue[{D[p[x, t]*v[x, t], x] + D[p[x, t], t] == 0, 
    D[p[x, t], x] + D[p[x, t]*v[x, t]^2, x] + D[p[x, t]*v[x, t], t] == 0, 
    p[0, t] == 2 - Exp[1 - t], p[1, t] == 1, p[x, 1] == 1, v[x, 1] == 0}, 
    {p, v}, {x, 0, 1}, {t, 1, 10}, Method -> {"MethodOfLines", "TemporalVariable" -> t, 
    "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 1000, 
    "MinPoints" -> 1000}}, MaxStepSize -> {Automatic, 5 10^-4}];

NDSolveValue::eerr: Warning ... (as above)

ContourPlot[Piecewise[{{pval[x, t], t > 1}}, 1], {x, 0, 1}, {t, 0, 10}]

ContourPlot[Piecewise[{{vval[x, t], t > 1}}, 0], {x, 0, 1}, {t, 0, 10}]

In addition to the obvious noise, which as before is due to numerical instability, a sheath has formed very near x == 0 for t > 5 (approximately). The cause of this sheath is unclear to me. Note that increasing the spatial and temporal resolution has little effect.