# NDSolve problem with describing thermodynamical evolution of a box

I use the following code to try and model 2D convection in time,

P = ρ[x, y, t]/(μ mu) kb T[x, y, t]  ;
e = P/(γ - 1);
cp = 5/2 kb/(μ mu);
Rgas = 8.3144598;
cv = 5/2 kb/(μ mu) - Rgas;
γ = cp/cv;
g = 28.02*9.81;
μ = 0.6163328197226503;
mu = 1.66053904*10^-27;
kb = 1.38064852*10^-23;
NDSolve[{
D[ρ[x, y, t]*u[x, y, t],
t] == -D[ρ[x, y, t]*(u[x, y, t])^2 + P, x],
D[ρ[x, y, t]*v[x, y, t],
t] == -D[ρ[x, y, t]*(v[x, y, t])^2 + P, y] +
g ρ[x, y, t],
D[ρ[x, y, t], t] == -D[ρ[x, y, t]*u[x, y, t], x] -
D[ρ[x, y, t]*v[x, y, t], y],
D[e[x, y, t], t] == -D[u[x, y, t]*e[x, y, t], x] -
D[v[x, y, t]*e[x, y, t], y] -
P[x, y, t]*(D[u[x, y, t], x] - D[v[x, y, t], y]),

v[0, y, t] == v[12000000, y, t], u[0, y, t] == u[12000000, y, t],
T[0, y, t] ==
T[12000000, y, t], ρ[0, y, t] == ρ[12000000, y, t],
D[u[x, y, t], y] /. y -> 4000000 == 0,
D[u[x, y, t], y] /. y -> 0 == 0, v[x, 4000000, t] == 0,
v[x, 0, t] == 0,
D[T[x, y, 0], y] == 0.41 g ρ [x, y, 0] T[x, y, 0]/P[x, y, 0],
D[P[x, y, 0], y] == g ρ[x, y, 0],
T[x, 0, 0] == 5770, ρ[x, 0, 0] == 1.42*10^-7*1.408*10^3
}, {u, v,
T, ρ}, {{x, 0, 12000000}, {y, 0, 4000000}, {t, 0, 100}}]


When evaluating this I get the following error

    NDSolve::dsvar: 0 cannot be used as a variable.


Is that because I use u[x,y,0] and NDSolve sees the zero as a variable or are there other problems?

• Your last two initial conditions, T[x, 0, 0] == 5770, \[Rho][x, 0, 0] == 1.42*10^-7*1.408*10^3 specify two variables to be zero, not just one. – KraZug May 17 '18 at 14:58
• Oh, and your two initial conditions on $u_y$ need a parenthesis around the replacement rule. – KraZug May 17 '18 at 15:02
• And you define P with no arguments, and then later try and give it three arguments. – KraZug May 17 '18 at 15:02
• Oh, and finally the arguments at the end should not be in a list, this is what is causing the specific error message. – KraZug May 17 '18 at 15:25
• @KraZug Would you like to sum up your comments in an answer? Or was it a simple mistake? – Kuba May 18 '18 at 12:16

After fixing up the obvious typos, you get to:

cp = 5/2 kb/(μ mu);
Rgas = 8.3144598;
cv = 5/2 kb/(μ mu) - Rgas;
γ = cp/cv;
g = 28.02*9.81;
μ = 0.6163328197226503;
mu = 1.66053904*10^-27;
kb = 1.38064852*10^-23;
P[x_, y_, t_] = ρ[x, y, t]/(μ mu) kb T[x, y, t];
e[x_, y_, t_] = P[x, y, t]/(γ - 1);

eqns = {D[ρ[x, y, t]*u[x, y, t], t] == -D[ρ[x, y, t]*(u[x, y, t])^2 + P[x, y, t], x],
D[ρ[x, y, t]*v[x, y, t], t] == -D[ρ[x, y, t]*(v[x, y, t])^2 + P[x, y, t], y] +  g ρ[x, y, t],
D[ρ[x, y, t], t] == -D[ρ[x, y, t]*u[x, y, t], x] - D[ρ[x, y, t]*v[x, y, t], y],
D[e[x, y, t], t] == -D[u[x, y, t]*e[x, y, t], x] - D[v[x, y, t]*e[x, y, t], y] - P[x, y, t]*(D[u[x, y, t], x] - D[v[x, y, t], y])}

bcs = {v[0, y, t] == v[12000000, y, t], u[0, y, t] == u[12000000, y, t],
T[0, y, t] == T[12000000, y, t], ρ[0, y, t] == ρ[12000000, y, t],
(D[u[x, y, t], y] /. y -> 4000000) == 0, (D[u[x, y, t], y] /. y -> 0) == 0,
v[x, 4000000, t] == 0, v[x, 0, t] == 0}

ics= {D[T[x, y, 0], y] == 0.41 g ρ[x, y, 0] T[x, y, 0]/P[x, y, 0],
D[P[x, y, 0], y] == g ρ[x, y, 0], T[x, y, 0] == 5770,
ρ[x, y, 0] == 1.42*10^-7*1.408*10^3}

NDSolve[Join[{eqns,bcs,ics}, {u, v, T, ρ}, {x, 0, 12000000}, {y, 0, 4000000}, {t, 0, 100}]


However, this then gives an error to say that you can't include $u_y$ in the boundary conditions as the equations don't include a higher order derivative in $y$ (i.e. $u_{yy}$). If you remove these derivatives, it gives a different error. I suspect that you need to tell Mathematica explicitly which method you want it to use to discretize your equations, and it would be very sensible to nondimensionalise properly (i.e. don't have your $x$ and $y$ ranges go up to the millions).