Skip to content

Instantly share code, notes, and snippets.

@Foadsf
Last active February 21, 2018 17:32
Show Gist options
  • Save Foadsf/0743534b0d766cabcd3eab49a79a1b8b to your computer and use it in GitHub Desktop.
Save Foadsf/0743534b0d766cabcd3eab49a79a1b8b to your computer and use it in GitHub Desktop.
Solving Navier-Stokes equations for a steady-state compressible viscous flow in a 2D axisymmetric step
lc = 0.03;
rc = 0.01;
xp = 0.01;
c = 0.005;
rp = rc - c;
lp = lc - xp;
T0 = 300;
eta0 = 1.846*10^-5;
P1 = 6*10^5 ;
P0 = 10^5;
cP = 1004.9;
cnu = 717.8;
R0 = cP - cnu;
Omega = RegionDifference[Rectangle[{0, 0}, {lc, rc}],
Rectangle[{xp, 0}, {xp + lp, rp}]];
Needs["NDSolve`FEM`"];
mesh = ToElementMesh[Omega, "MaxBoundaryCellMeasure" -> 0.00001,
MaxCellMeasure -> {"Length" -> 0.0008},
"MeshElementConstraint" -> 20, MeshQualityGoal -> "Maximal"];
eqn1 = D[rho[x, r]*nux[x, r], x] + D[r*rho[x, r]*nur[x, r], r]/r == 0;
eqn2 = D[rho[x, r]*nux[x, r]^2 + R0*rho[x, r]*T[x, r], x] +
D[r*(rho[x, r]*nux[x, r]*nur[x, r] + eta0*D[nux[x, r], r]), r]/
r == 0;
eqn3 = D[rho[x, r]*nux[x, r]*nur[x, r] + eta0*D[nur[x, r], x], x] +
D[r*(rho[x, r]*nur[x, r]^2 + R0*rho[x, r]*T[x, r]), r]/r == 0;
eqn4 = cnu*
rho[x, r]*(nux[x, r]*D[T[x, r], x] + nur[x, r]*D[T[x, r], r]) +
R0*rho[x, r]*
T[x, r]*(D[nux[x, r], x] +
D[r*nur[x, r], x]/r) + (2*D[nux[x, r], x]^2 +
2*D[nur[x, r], r]^2 + (D[nux[x, r], r] +
D[nur[x, r],
x])^2 - ((D[nux[x, r], x] + D[r*nur[x, r], x]/r)^2)*2/3)*
eta0 == 0;
eqns = {eqn1, eqn2, eqn3, eqn4};
bc1 = R0* rho[0, r]*T0 == P1;
bc2 = R0 *rho[lc, r]*T0 == P0;
bc3 = DirichletCondition[{nur[x, 0] == 0, D[nur[x, r], r] == 0,
D[nux[x, r], r] == 0, D[rho[x, r], r] == 0, D[T[x, r], r] == 0},
r == 0 && (0 <= x <= xp)];
bc4 = DirichletCondition[{nur[x, r] == 0,
nux[x, r] ==
0}, (0 <= r <= rp && x == xp) || (r == rp &&
xp <= x <= xp + lp) || (r == rc && 0 <= x <= lc)];
bcs = {bc1, bc2, bc3, bc4};
{nuxsol, nursol, rhosol, Tsol} =
NDSolveValue[{eqns, bcs}, {nux, nur, rho, T}, {x, r} \[Element] mesh,
Method -> {"FiniteElement",
"InterpolationOrder" -> {nux -> 2, nur -> 2, rho -> 1, T -> 1},
"IntegrationOrder" -> 5}]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment