Solving 3 nonlinear, transient, coupled PDEs [closed]

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.

closed as off-topic by xzczd, Henrik Schumacher, Carl Lange, J. M. will be back soon♦Mar 9 at 9:37

This question appears to be off-topic. The users who voted to close gave these specific reasons:

• "The question is out of scope for this site. The answer to this question requires either advice from Wolfram support or the services of a professional consultant." – J. M. will be back soon
• "This question cannot be answered without additional information. Questions on problems in code must describe the specific problem and include valid code to reproduce it. Any data used for programming examples should be embedded in the question or code to generate the (fake) data must be included." – xzczd, Henrik Schumacher, Carl Lange
If this question can be reworded to fit the rules in the help center, please edit the question.

• If you share the code that you tried in copyable form, I am sure that somebody is capable and willing to help you... – Henrik Schumacher Mar 7 at 23:55
• @Liam then unfortunately nobody will be able to help you. Why can’t you copy the code?! – MarcoB Mar 8 at 2:13
• If you can't provide the specific code, I'm afraid this question is off-topic here. Anyway, since you've mentioned shock wave, have you read the following posts?: mathematica.stackexchange.com/q/11748/1871 mathematica.stackexchange.com/q/24417/1871 – xzczd Mar 8 at 6:44
• @Liam Can you express $e_x$ through $\rho _x, T_x$? – Alex Trounev Mar 8 at 14:16
• "Also, I can't post the code explicitly due to where I work.." Oh come on: You've already posted the differential equations in $\LaTeX$ form and they are certainly not super-secret. The point is that us people here use our free time and don't charge you and your employer any money. So the least effort we expect from you is to provide rudimental code to start with so that we don't have to retype everything. Really, if you cannot do that, your only alternative is to hire a professional consultant. – Henrik Schumacher Mar 9 at 7:22

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.

• I received several messages including NDSolve::ndsz: At t == 0.0004000246070339927, step size is effectively zero; singularity or stiff system suspected.` – Alex Trounev Mar 8 at 13:30
• There is a typo in PDE1. Must be PDE1 = D[r[t, x], t] + u[t, x]*D[r[t, x], x] + r[t, x]*D[u[t, x], x] == 0. And in PDE3 parameter $e_x=\partial e(T)/\partial x$. – Alex Trounev Mar 8 at 13:43
• @AlexTrounev First of all thank you for the correction, secondly I am not getting any such warnings. The only one I am getting is of boundary and initial conditions are inconsistent – zhk Mar 8 at 13:56
• @AlexTrounev The third point $e_x$, if your suggestion is true then we need to know what exactly is $e(T)$ and we will be needing another D/A equation for it? – zhk Mar 8 at 13:59