# 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. Commented May 17, 2018 at 14:58
• Oh, and your two initial conditions on $u_y$ need a parenthesis around the replacement rule. Commented May 17, 2018 at 15:02
• And you define P with no arguments, and then later try and give it three arguments. Commented May 17, 2018 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. Commented May 17, 2018 at 15:25
• @KraZug Would you like to sum up your comments in an answer? Or was it a simple mistake?
– Kuba
Commented May 18, 2018 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).