Skip to main content

finite element method - Solving Navier-Stokes equations for a steady-state compressible viscous flow in a 2D axisymmetric step


Note: you may apply or follow the edits on the code here in this GitHub Gist


I'm trying to follow this post to solve Navier-Stokes equations for a compressible viscous flow in a 2D axisymmetric step. The geometry is :


enter image description here



lc = 0.03;
rc = 0.01;
xp = 0.01;
c = 0.005;
rp = rc - c;
lp = lc - xp;
Subscript[T, 0] = 300;
Subscript[\[Eta], 0] = 1.846*10^-5;
Subscript[P, 1] = 6*10^5 ;
Subscript[P, 0] = 10^5;

Subscript[c, P] = 1004.9;
Subscript[c, \[Nu]] = 717.8;
Subscript[R, 0] = Subscript[c, P] - Subscript[c, \[Nu]];
\[CapitalOmega] = RegionDifference[
   Rectangle[{0, 0}, {lc, rc}], 
   Rectangle[{xp, 0}, {xp + lp, rp}]];

And meshing:


Needs["NDSolve`FEM`"];
mesh = ToElementMesh[\[CapitalOmega], 

   "MaxBoundaryCellMeasure" -> 0.00001, 
   MaxCellMeasure -> {"Length" -> 0.0008}, 
   "MeshElementConstraint" -> 20, MeshQualityGoal -> "Maximal"][
  "Wireframe"]

enter image description here


Where the model is axisymmetric around the x axis, the governing equations including conservation equations of mass, momentum and heat can be written as:


$$ \frac{\partial}{\partial x}\left( \rho \nu_x \right)+\frac{1}{r}\frac{\partial }{\partial r}\left(r \rho \nu_r\right)=0 \tag{1}$$


$$\frac{\partial}{\partial x}\left( \rho \nu_x^2+\mathring{R} \rho T \right)+\frac{1}{r}\frac{\partial}{\partial r}\left( r \left( \rho \nu_r \nu_x + \eta \frac{\partial \nu_x}{\partial r} \right)\right) \tag{2}$$


$$ \frac{\partial}{\partial x}\left( \rho \nu_x \nu_r+\eta \frac{\partial \nu_r}{\partial x} \right)+ \frac{1}{r}\frac{\partial}{\partial r}\left( r \left( \rho \nu_r ^2 +\mathring{R} \rho T \right) \right)=0 \tag{3}$$



$$\rho c_\nu\left( \nu_x \frac{\partial T}{\partial x} + \nu_r \frac{\partial T}{\partial r} \right)+ \mathring{R} \rho T \left( \frac{1}{r}\frac{\partial}{\partial r} \left( r \nu_r \right)+ \frac{\partial \nu_x}{\partial x} \right)+ \eta \left( 2 \left( \frac{\partial \nu_x}{\partial x} \right)^2+ 2 \left( \frac{\partial \nu_r}{\partial r} \right)^2+ \left( \frac{\partial \nu_r}{\partial x}+ \frac{\partial \nu_x}{\partial r} \right)^2 \\ -\frac{2}{3}\left( \frac{1}{r} \frac{\partial}{\partial r}\left( r \nu_r \right) + \frac{\partial \nu_x}{\partial x} \right)^2 \right)=0 \tag{4}$$


eqn1 = D[\[Rho][x, r]*Subscript[\[Nu], x][x, r], x] + 
    D[r*\[Rho][x, r]*Subscript[\[Nu], r][x, r], r]/r == 0 ;
eqn2 = D[\[Rho][x, r]*Subscript[\[Nu], x][x, r]^2 + 
      Subscript[R, 0] \[Rho][x, r]*T[x, r], x] + 
    D[r*(\[Rho][x, r]*Subscript[\[Nu], x][x, r]*
          Subscript[\[Nu], r][x, r] + 
         Subscript[\[Eta], 0]*D[Subscript[\[Nu], x][x, r], r]), r]/
     r == 0 ;
eqn3 = D[\[Rho][x, r]*Subscript[\[Nu], x][x, r]*

       Subscript[\[Nu], r][x, r] + 
      Subscript[\[Eta], 0]*D[Subscript[\[Nu], r][x, r], x], x] + 
    D[r*(\[Rho][x, r]*Subscript[\[Nu], r][x, r]^2 + 
         Subscript[R, 0] \[Rho][x, r]*T[x, r]), r]/r == 0;
eqn4 = Subscript[
     c, \[Nu]]*\[Rho][x, 
      r]*(Subscript[\[Nu], x][x, r]*D[T[x, r], x] + 
       Subscript[\[Nu], r][x, r]*D[T[x, r], r]) + 
    Subscript[R, 0]*\[Rho][x, r]*
     T[x, r]*(D[Subscript[\[Nu], x][x, r], x] + 

       D[r*Subscript[\[Nu], r][x, r], x]/r) + (2*
        D[Subscript[\[Nu], x][x, r], x]^2 + 
       2*D[Subscript[\[Nu], r][x, r], 
          r]^2 + (D[Subscript[\[Nu], x][x, r], r] + 

          D[Subscript[\[Nu], r][x, r], 
           x])^2 - ((D[Subscript[\[Nu], x][x, r], x] + 
            D[r*Subscript[\[Nu], r][x, r], x]/r)^2)*2/3)*
     Subscript[\[Eta], 0] == 0;
eqns = {eqn1, eqn2, eqn3, eqn4};


And the boundary conditions are:



  1. constant pressure at inlet

  2. constant pressure at outlet

  3. axis of symmetry

  4. no slip


Implemented as


bc1 = Subscript[R, 0] \[Rho][0, r]*Subscript[T, 0] == Subscript[P, 1] 

bc2 = Subscript[R, 0] \[Rho][lc, r]*Subscript[T, 0] == Subscript[P, 0]
bc3 = DirichletCondition[{Subscript[\[Nu], r][x, 0] == 0, 
   D[Subscript[\[Nu], r][x, r], r] == 0, 
   D[Subscript[\[Nu], x][x, r], r] == 0, D[\[Rho][x, r], r] == 0, 
   D[T[x, r], r] == 0}, r == 0 && (0 <= x <= xp  )] 
bc4 = DirichletCondition[{Subscript[\[Nu], r][x, r] == 0, 
    Subscript[\[Nu], x][x, r] == 
     0}, (0 <= r <= rp && x == xp ) || (r == rp && 
      xp <= x <= xp + lp)  || (r == rc && 0 <= x <= lc) ] == 0
bcs = {bc1, bc2, bc3, bc4};


When I try to solve the equations:


{\[Nu]xsol, \[Nu]rsol, \[Rho]sol, Tsol} = 
  NDSolveValue[{eqns, , bcs}, {Subscript[\[Nu], x], Subscript[\[Nu], 
    r], \[Rho], T}, {x, r} \[Element] mesh, 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {Subscript[\[Nu], x] -> 2, 
       Subscript[\[Nu], r] -> 2, \[Rho] -> 1, T -> 1}, 
     "IntegrationOrder" -> 5}];


I get the errors:



NDSolveValue::femnr: {x,r}[Element] is not a valid region specification.



and



Set::shape: Lists {[Nu]xsol,[Nu]rsol,Tsol,[Rho]sol} and NDSolveValue[<<1>>] are not the same shape.



Googling the errors does not offer that much of help (e.g. here). I would appreciate if you could help me know What is the issue and how I can solve it.


P.S.1. The NDSolveValue femnr error was caused by [ "Wireframe"] term at the end of meshing command changing it to



mesh = ToElementMesh[\[CapitalOmega], 
   "MaxBoundaryCellMeasure" -> 0.00001, 
   MaxCellMeasure -> {"Length" -> 0.0008}, 
   "MeshElementConstraint" -> 20, MeshQualityGoal -> "Maximal"];
mesh["Wireframe"]

resolves the issue.


P.S.2. There is an extra ==0 at the end of boundary condition 4 it was edited to:


bc4 = DirichletCondition[{Subscript[\[Nu], r][x, r] == 0, 
    Subscript[\[Nu], x][x, r] == 

     0}, (0 <= r <= rp && x == xp) || (r == rp && 
      xp <= x <= xp + lp) || (r == rc && 0 <= x <= lc)];

at this moment the second error still persists and a new error was added:



NDSolveValue::deqn Equation or list of equations expected instead of Null in the first argument ...



P.S.3 There were multiple issues. So I decided to use this Github Gist to further edit the code.




Comments

Post a Comment

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
Post a Comment
Read more

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more