0

I need to solve a system of DE's in 2 dimensions over time, The problem is that mathematica complains when I have a boundary condition of the same order as the differential order of the DE's, so I tried using the finiteelement method since it worked for other people who had similar problems, but it just returns itself as if nothing ever happened

P[x_, y_, t_] = e[x, y, t]/(γ - 1) ; 
e[x_, y_, t_] = (γ - 1) ρ[x, y, t]/(μ mu ) kb T[x, y, t];
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;
sol1 = NDSolve[{
D[ρ[x, y, t]*u[x, y, t], 
 t] == -D[ρ[x, y, t]*u[x, y, t]*u[x, y, t] + P[x, y, t], x] -
  D[ρ[x, y, t]*u[x, y, t]*
    v[x, y, t], 
  y],
D[ρ[x, y, t]*v[x, y, t], 
 t] == -D[ρ[x, y, t]*v[x, y, t]*
     u[x, y, t], 
   y] - D[ρ[x, y, t]*v[x, y, t]*v[x, y, t] + 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]),

v[0, y, t] == v[12*10^6, y, t],
u[0, y, t] == u[12*10^6, y, t],
T[0, y, t] == T[12*10^6, y, t],
ρ[0, y, t] == ρ[12*10^6, y, t],

(D[u[x, y, t], x] /. x -> 0) == (D[u[x, y, t], x] /. x -> 12000000),
(D[v[x, y, t], x] /. x -> 0) == (D[v[x, y, t], x] /. x -> 12000000),
D[u[0, y, t], y] == D[u[12000000, y, t], y],
D[v[0, y, t], y] == D[v[12000000, y, t], y],
D[u[0, y, t], t] == D[u[12000000, y, t], t],
D[v[0, y, t], t] == D[v[12000000, y, t], t],

(D[T[x, y, t], x] /. x -> 0) == (D[T[x, y, t], x] /. x -> 12000000),
(D[ρ[x, y, t], x] /. x -> 0) == (D[ρ[x, y, t], x] /. 
x -> 12000000),
D[T[0, y, t], y] == D[T[12000000, y, t], y],
D[ρ[0, y, t], y] == D[ρ[12000000, y, t], y],
D[T[0, y, t], t] == D[T[12000000, y, t], t],
D[ρ[0, y, t], t] == D[ρ[12000000, y, t], t],

(D[e[x, y, t], x] /. x -> 0) == (D[e[x, y, t], x] /. x -> 12000000),
D[e[0, y, t], y] == D[e[12000000, y, t], y],
D[e[0, y, t], t] == D[e[12000000, y, t], t],


e[x, 0, t] == 3.83767261162,
v[x, 4000000, t] == 0,
v[x, 0, t] == 0,
(D[u[x, y, t], y] /. y -> 0) == 0,
(D[u[x, y, t], y] /. y -> 4000000) == 0,

v[x, y, 0] == 0,
u[x, y, 0] == 0,
T[x, y, 0] == 5770 + 0.00835414960707927 y,
ρ[x, y, 0] == 
1.42*10^-7*1.408*10^3 + 7.3561137493644*10^-10 y
},
{u, v, T, ρ}, {x, 0, 12000000}, {y, 0, 4000000}, {t, 0, 100},
Method -> "FiniteElement" ]
Hidde Rinsema
  • 43
  • 4
  • @MichaelSeifert Yes, there is supposed to be a break between the second and third line of the code and running it on a clean kernel gives the same result as described in the question itself. And there shouldn't be any need for a differnential equation regarding just T. – Hidde Rinsema May 21 '18 at 13:41
  • The first instance of NDSolve runs fine, and its results are not used later. So, I recommend that you delete it from your question to simplify it a bit. Readers tend to avoid questions in which the length of the code obscures the fundamental issues to be resolved. Your second instance of NDSolve fails because T does not appear in any of its PDEs. Perhaps, you miscopied your code. – bbgodfrey May 21 '18 at 22:10
  • Further to bbgodfrey's point, you have 4 equations, but T only appears in the boundary conditions. – SPPearce May 22 '18 at 10:11
  • Aha, I saw the problem, T wasn't added in the definition for e, so it wasn' t in the differential equations, I've updated it, so now it should work as desctribed – Hidde Rinsema May 22 '18 at 11:22
  • I cannot use second order derivatives, so I have updated the problem the be about the problem I have without second order derivatives – Hidde Rinsema May 22 '18 at 13:46
  • In general, specifying boundary conditions of the same order as the differential equations themselves either leads to inconsistencies in the equations (i.e., no solutions exist) or doesn't actually constrain your solution (i.e., a large class of solutions exists.) Neither option is conducive to numerical solving. This is a general property of PDEs, and not really a Mathematica issue; you might try asking about this question over at Math.SE or Physics.SE instead. – Michael Seifert May 22 '18 at 14:52
  • @MichaelSeifert Yes I realized that, but then I updated the boundary conditions, to include the ones in T and [Rho] but then it does nothing. – Hidde Rinsema May 22 '18 at 16:16
  • Last time I checked, Mathematica couldn't solve non-linear PDEs using the finite element method, which is probably why it's coming back unevaluated. It's possible this has changed in v11, though; I haven't upgraded yet. – Michael Seifert May 22 '18 at 19:40

0 Answers0