# 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.

• If you share the code that you tried in copyable form, I am sure that somebody is capable and willing to help you... Commented Mar 7, 2019 at 23:55
• @Liam then unfortunately nobody will be able to help you. Why can’t you copy the code?! Commented Mar 8, 2019 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 Commented Mar 8, 2019 at 6:44
• @Liam Can you express $e_x$ through $\rho _x, T_x$? Commented Mar 8, 2019 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. Commented Mar 9, 2019 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. Commented Mar 8, 2019 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$. Commented Mar 8, 2019 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
Commented Mar 8, 2019 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
Commented Mar 8, 2019 at 13:59