Solving 3 nonlinear, transient, coupled PDEs [closed]

Question

I'm trying to solve the fluid equations, continuity equations for mass, momentum, and energy, for a transient shockwave in 1D planar geometry. A piston is traveling at some velocity $u_p$ where $u_P$ $>$ $c_{sound}$ and strikes quiescent air. I am attempting to solve the problem in the reference frame of the piston, with the air coming at the piston with velocity $- u_p$.

These are the non-dimensional equations I am trying to solve with pressure and internal energy have been substituted with the Ideal Gas EOS:

with $\rho$ being the density, $u$ being the particle velocity, and $T$ being the temperature. These are the dependent variables. $x$ and $t$ are the independent variables. Subscripts denote derivatives and exempting the $0$ subscript and $x_i$ constants. The prefactors with the subscript $0$, like $p_0$, etc. are constants.

The boundary conditions and initial conditions and solving over a domain of 0 to 30 for x and 0 to 20 for t:

I honestly don't know what boundary condition to apply at $x = 0$ for $\rho$. If the 2nd derivative of $\rho$ appeared I would apply a Neumann condition.

These boundary conditions and initial conditions would approximate the shock forming and then propagating within the domain that NDSolve (the command I've tried to use) would solve for. This is also the reason that $t_{max}$ is less than $x_{max}$, so the shockwave doesn't reach the right hand boundary. The problem I've run into is that I don't know what method would be the best in Mathematica to solve these equations. I've looked around on the stack exchange and extensively through the Wolfram documentation. I could have missed the answer no matter how thoroughly I've searched, and I haven't really found an answer on what methods are the best to solve this equation. I'm honestly a little lost.

Any help would be greatly appreciated.

Answer 1

Please next time show some efforts.

PDE1 = D[r[t, x], t] + u[t, x]*D[r[t, x], x] + r[t, x]*D[u[t, x], x] == 0

PDE2 = r[t, x]*D[u[t, x], t] + p0/(r0*u0^2)*D[r[t, x]*T[t, x], x] + 
   r[t, x]*u[t, x]*D[u[t, x], x] - xu/x0*D[u[t, x], x, x] == 0

PDE3 = r[t, x]*D[T[t, x], t] + r[t, x]*u[t, x]*ex + 
   p0/(r0*e0)*r[t, x]*T[t, x]*D[u[t, x], x] - 
   xu*u0^2/(x0*e0)*D[u[t, x], x]*D[u[t, x], x] - 
   xk/x0*D[T[t, x], x, x] == 0

p0 = 1; r0 = 1; u0 = 1; xu = 1; e0 = 1; x0 = 1; xk = 1; ex = 1; up = 0.5;

mol[n_Integer, o_: "Pseudospectral"] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}}
mol[tf : False | True, sf_: Automatic] := {"MethodOfLines", 
  "DifferentiateBoundaryConditions" -> {tf, "ScaleFactor" -> sf}}

ics = {u[0, x] == up, T[0, x] == 1, r[0, x] == 1};
With[{xmax = 5}, 
  bcs = {{u[t, 0] == 0, r[t, 0] == 0, 
     Derivative[0, 1][T][t, 0] == 0}, {u[t, xmax] == -up, 
     T[t, xmax] == 1, r[t, xmax] == 1}}];

sol = NDSolve[{PDE1, PDE2, PDE3, bcs, ics}, {r, u, T}, {t, 0, 1}, {x, 
   0, 5}, Method -> Union[mol[70, 4], mol[True, 100]]]

Plot3D[{r[t, x], u[t, x], T[t, x]} /. sol, {t, 0, 1}, {x, 0, 5}]

Note: Before you ask further, provide the missing boundary condition and numerical values for the parameters.