
[Sponsors] 
Mathematica: steadystate compressible air viscous flow in a 2D axisymmetric step 

LinkBack  Thread Tools  Search this Thread  Display Modes 
February 21, 2018, 12:39 
Mathematica: steadystate compressible air viscous flow in a 2D axisymmetric step

#1  
Member
Foad
Join Date: Aug 2017
Posts: 58
Rep Power: 8 
I have also posted my question here in StackExchange and also my code here in this GitHub Gist.
I'm trying to follow this post using Mathematica to solve NavierStokes equations for a compressible viscous flow in a 2D axisymmetric step. The geometry is in the attachment. Defining variables, geometry and meshing: Code:
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"]; Where the model is axisymmetric around the x axis, the governing equations including conservation equations of mass, momentum and heat can be written as: forth equation is big so I add it as a image in the attachments. Implementing the equation is Mathematica: Code:
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}; And the boundary conditions are:
Implemented in Mathematica as: Code:
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}; Code:
{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}] Quote:
Quote:


Tags 
air, axisymmetric, compressible, mathematica, viscous 
Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Steady state free surface flow in rotational closed system  hilde  CFX  5  January 21, 2015 16:47 
Constant velocity of the material  Sas  CFX  15  July 13, 2010 08:56 
Could anybody help me see this error and give help  liugx212  OpenFOAM Running, Solving & CFD  3  January 4, 2006 18:07 
steady state creeping flow  dominik  Main CFD Forum  4  March 29, 2004 09:29 
About the difference between steady and unsteady problems  Lisa  Main CFD Forum  11  July 5, 2000 14:37 