9
$\begingroup$

I am using Mathematica's (v12.0) NDSolveValue method to solve a finite element method problem (Navier-Stokes equations for a compressible gas). During the initialization process at t==0, I repeatedly get the following error- "The FEMDampingElements operator failed". What does the error mean and does any one have any suggestions as how to resolve it.

<<NDSolve`FEM`

ifctn = Interpolation[{{2.01,0.6},{1.5,0.6},{0.5,0.35},{-1.0,1.0},{-2.01,1.2}}];
Ω = ImplicitRegion[Abs[y] <= Abs[ ifctn[ Max[ -2.01, Min[ 2.01, x ] ] ] ], {{x,-2,2}, {y,-1.2,1.2}} ];

softrampfunction[min_,max_,c_,r_] := Function[t,Max[0.,(min*Exp[c*r]+max*Exp[r*t])/(Exp[c*r]+Exp[r*t])]]
sf = softrampfunction[0,1,3,3];

af = softrampfunction[0,1,3,4];
aperaturefunction[w_] := Function[y,af[(10*y/w)+5.5]-af[(10*y/w)+0.5]]
aperature = aperaturefunction[1.2];

χ [t_,x_,y_] := D[u[t,x,y],{t,0},{x,1},{y,0}]*D[u[t,x,y],{t,0},{x,1},{y,0}] + D[u[t,x,y],{t,0},{x,0},{y,1}]*D[v[t,x,y],{t,0},{x,1},{y,0}] - 2/3 D[u[t,x,y],{t,0},{x,1},{y,0}]*D[v[t,x,y],{t,0},{x,0},{y,1}] +
                 D[v[t,x,y],{t,0},{x,1},{y,0}]*D[v[t,x,y],{t,0},{x,1},{y,0}] + D[v[t,x,y],{t,0},{x,1},{y,0}]*D[u[t,x,y],{t,0},{x,0},{y,1}] - 2/3 D[v[t,x,y],{t,0},{x,0},{y,1}]*D[u[t,x,y],{t,0},{x,1},{y,0}];

parameters = {μ->Sqrt[Pr/Ra],ϵ->1/Sqrt[Pr*Ra],Μ->M,Subscript[R,m]->Ra/M,Subscript[C,v]->1,λ->1.2,κ->1}/.{M->18.01528,Pr->7.1,Ra->2*10^5}; 
parameters0 = {μ->1,Μ->1,Subscript[R,m]->1,Subscript[C,v]->10,->1.4,κ->1};

op = Flatten[{
       D[{u[t,x,y],v[t,x,y]},{t,1},{x,0},{y,0}]ρ[t,x,y]+Grad[ρ[t,x,y],{x,y}]+μ(Div[Grad[{u[t,x,y],v[t,x,y]},{x,y}],{x,y}]+1/3 Grad[Div[{u[t,x,y],v[t,x,y]},{x,y}],{x,y}]),

       D[ρ[t,x,y],{t,1},{x,0},{y,0}] + ρ[t,x,y]*Div[{u[t,x,y],v[t,x,y]},{x,y}],

       1/(λ-1)*(D[p[t,x,y],{t,1},{x,0},{y,0}] - λ*Subscript[R,m]*Τ[t,x,y]D[ρ[t,x,y],{t,1},{x,0},{y,0}]) - χ[t,x,y] - κ Inactive[Laplacian][Τ[t,x,y],{x,y}],

       ρ[t,x,y]*Subscript[C,v]/M D[Τ[t,x,y],{t,1},{x,0},{y,0}] - Subscript[R,m]*ρ[t,x,y]*Τ[t,x,y]*Div[{u[t,x,y],v[t,x,y]},{x,y}] + χ[t,x,y] + κ*Inactive[Laplacian][Τ[t,x,y],{x,y}
       ]
     }]/.parameters;

ic = { u[0,x,y]==0,v[0,x,y]==0,ρ[0,x,y]==1 ,p[0,x,y]==0 ,Τ[0,x,y]==0 };

Subscript[Γ,B] = { 0.0,0.0,0.0,0.0,0.0 };  (* NeumannValue[-100sf[t],(x>=2)&&(Abs[y]<=0.2)] *)

Subscript[Γ,D] = { DirichletCondition[{u[t,x,y]==0.,v[t,x,y]==0.},(-2.<x<2.)&&(y^2>=ifctn[x]^2)],
                   DirichletCondition[u[t,x,y]==sf[t]*aperature[y],(x>=2.)],
                   DirichletCondition[v[t,x,y]==0.,(x<=-2.)],
                   DirichletCondition[p[t,x,y]==0.,(x<=-2.)],
                   DirichletCondition[Τ[t,x,y]==1000.,(x>=2.)] };

pde = op == Subscript[Γ,B];


Monitor[AbsoluteTiming[{Subscript[U,vel],Subscript[V,vel],Subscript[D,ensity],Subscript[P,ressure],Subscript[T,emperatue]} = NDSolveValue[{pde,Subscript[Γ,D],ic},{u,v,ρ,p,Τ},{x,y}∈ToElementMesh    [Ω,MaxCellMeasure->0.005],{t,0,10},
Method->{"FiniteElement","InterpolationOrder"->{u->2,v->2,ρ->2,p->2,Τ->1}},
EvaluationMonitor:>(currentTime=Row[{"t = ",CForm[t]}])];],currentTime]
$\endgroup$
  • $\begingroup$ Please provide the Mathematica code itself (not an image of it) that produces this error. By the way, I presume that you mean NDSolveValue , not NSolveValue . $\endgroup$ – bbgodfrey Sep 29 '19 at 19:49
  • $\begingroup$ Unfortunately, you code has syntax errors when I copy it. If you fix those I can have another look. Ping me with @user21. $\endgroup$ – user21 Sep 30 '19 at 12:40
  • $\begingroup$ @CarlVoss There are some errors with signs in op. What does this model describe? $\endgroup$ – Alex Trounev Oct 10 '19 at 17:53
  • $\begingroup$ @Alex Trounev Thank you for the help you have been providing me. Yes, there are sign errors in the 'op' variable which I have corrected (I hope). I am attempting to solve Navier-Stokes equations for a compressible gas through a nozzle. I am new at programming computational fluid dynamics problems. If you are interested in the equations that I am using, please see a course that I have found in html format: Fluid Mechanics, Richard Fitzpatrick (Professor of Physics at The University of Texas at Austin) (farside.ph.utexas.edu/teaching/336L/Fluidhtml/Fluidhtml.html) $\endgroup$ – Carl Voss Oct 11 '19 at 2:45
10
$\begingroup$

Using my 40+ years of experience modeling unsteady supersonic flows I can offer a working FEM code for isentropic viscous flows. Here you can add the equation for temperature with $\chi $ function, but this almost does not affect the solution. A flow with the inlet/outlet pressure ratio $p_{in}/p_{out}=25$ is considered here. Mach number at nozzle exit $M=3.77$.

<< NDSolve`FEM`
ifctn = Interpolation[{{2.01, 0.6}, {1.5, 0.6}, {0.5, 0.35}, {-1.0, 
     1.0}, {-2.01, 1.2}}];
\[CapitalOmega] = 
  ImplicitRegion[
   Abs[y] <= 
    Abs[ifctn[Max[-2.01, Min[2.01, x]]]], {{x, -2, 2}, {y, -1.2, 
     1.2}}];
mesh = ToElementMesh[\[CapitalOmega], "MaxCellMeasure" -> 0.002];
mesh["Wireframe"]

q = 1.4;
k = 190; Re0 = 1; U0 = 1; M0 = 1; Re1 = Re0/M0^2; rhoin = 1; rhout = 
 1/10; t0 = 1/25;
UX[0][x_, y_] := 0;
VY[0][x_, y_] := 0;
\[CapitalRho][0][x_, y_] := rhoin;
Do[
   {UX[i], VY[i], \[CapitalRho][i]} = 
    NDSolveValue[{{Inactive[
            Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
              u[x, y], {x, y}]), {x, y}] + 
          Re1*(Abs[\[CapitalRho][i - 1][x, y]]^(q - 1))*
\!\(\*SuperscriptBox[\(\[Rho]\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y]/\[CapitalRho][i - 1][x, y] + 
          Re0*UX[i - 1][x, y]*D[u[x, y], x] + 
          Re0*VY[i - 1][x, y]*D[u[x, y], y] + 
          Re0*(u[x, y] - UX[i - 1][x, y])/t0, 
         Inactive[
            Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
              v[x, y], {x, y}]), {x, y}] + 
          Re1*(Abs[\[CapitalRho][i - 1][x, y]^q])*
\!\(\*SuperscriptBox[\(\[Rho]\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y]/\[CapitalRho][i - 1][x, y] + 
          Re0*UX[i - 1][x, y]*D[v[x, y], x] + 
          Re0*VY[i - 1][x, y]*D[v[x, y], y] + 
          Re0*(v[x, y] - VY[i - 1][x, y])/t0, 
         D[\[CapitalRho][i - 1][x, y]*u[x, y], x] + 
          D[\[CapitalRho][i - 1][x, y]*v[x, y], 
           y] + (\[Rho][x, y] - \[CapitalRho][i - 1][x, y])/t0} == {0,
          0, 0} /. \[Mu] -> 1/100, {
       DirichletCondition[{v[x, y] == 0, \[Rho][x, y] == rhoin}, 
        x == 2.], 
       DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, -2 < x < 2],
       DirichletCondition[{\[Rho][x, y] == rhout, v[x, y] == 0}, 
        x == -2]}}, {u, v, \[Rho]}, {x, y} \[Element] mesh, 
     Method -> {"FiniteElement", 
       "InterpolationOrder" -> {u -> 2, v -> 2, \[Rho] -> 1}}], {i, 1,
     k}]; // Quiet

The flow velocity normalized to the speed of sound at the input and the density at the last step

Show[DensityPlot[
  Norm[{UX[190][x, y], VY[190][x, y]}], {x, y} \[Element] mesh, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  AspectRatio -> Automatic], 
 StreamPlot[{UX[k][x, y], VY[k][x, y]}, {x, y} \[Element] mesh, 
  StreamStyle -> LightGray, StreamPoints -> Fine, 
  AspectRatio -> Automatic]]

DensityPlot[\[CapitalRho][190][x, y], {x, y} \[Element] mesh, 
 ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
 AspectRatio -> Automatic]

Figure 1 Density animation

lst = Table[
   DensityPlot[\[CapitalRho][i][x, y]
    , {x, y} \[Element] mesh, 
    ColorFunction -> "Rainbow", Frame -> None, 
    AspectRatio -> Automatic, PlotLegends -> Automatic, 
    PlotRange -> All], {i, 5, 190, 5}];
ListAnimate[lst]

Figure 2

Animation of flow velocity

lst1 = Table[
   DensityPlot[
    Norm[{UX[i][x, y], VY[i][x, y]}], {x, y} \[Element] mesh, 
    ColorFunction -> "Rainbow", Frame -> None, 
    AspectRatio -> Automatic, PlotLegends -> Automatic, 
    PlotRange -> All], {i, 5, 190, 5}];
ListAnimate[lst1]

Figure 3

Let us show now that a similar numerical solution can be obtained using the method proposed by user21. For this, we modify the author’s code so that the equations fully correspond to my code

<< NDSolve`FEM`

ifctn = Interpolation[{{2.01, 0.6}, {1.5, 0.6}, {0.5, 0.35}, {-1.0, 
     1.0}, {-2.01, 1.2}}];
\[CapitalOmega] = 
  ImplicitRegion[
   Abs[y] <= 
    Abs[ifctn[Max[-2.01, Min[2.01, x]]]], {{x, -2, 2}, {y, -1.2, 
     1.2}}];

softrampfunction[min_, max_, c_, r_] := 
 Function[t, 
  Max[0., (min*Exp[c*r] + max*Exp[r*t])/(Exp[c*r] + Exp[r*t])]]
sf = softrampfunction[0, 1, 3, 3];



op = Flatten[{D[{u[t, x, y], v[t, x, y]}, {t, 
        1}] + {{u[t, x, y], v[t, x, y]}.Grad[
         u[t, x, y], {x, y}], {u[t, x, y], v[t, x, y]}.Grad[
         v[t, x, y], {x, y}]} + 
      Grad[p[t, x, y], {x, y}]/\[Rho][t, x, 
        y] - \[Mu] (Div[
          Grad[{u[t, x, y], v[t, x, y]}, {x, y}], {x, y}] + 
         0/3 Grad[Div[{u[t, x, y], v[t, x, y]}, {x, y}], {x, y}]), 
     D[\[Rho][t, x, y], t] + 
      Div[\[Rho][t, x, y]*{u[t, x, y], v[t, x, y]}, {x, y}], 
     D[p[t, x, y], {t, 1}] + {u[t, x, y], v[t, x, y]}.Grad[
        p[t, x, y], {x, y}] + 
      q  (p[t, x, 
          y] Div[{u[t, x, y], v[t, x, y]}, {x, y}])}] /. {\[Mu] -> 
     1/100, q -> 1.4};
ic = {u[0, x, y] == 0, v[0, x, y] == 0, \[Rho][0, x, y] == 1, 
   p[0, x, y] == 1};

Subscript[\[CapitalGamma], B] = {0.0, 0.0, 0.0, 0.}; 
Subscript[\[CapitalGamma], 
  D] = {DirichletCondition[{u[t, x, y] == 0., v[t, x, y] == 0.}, -2. <
     x < 2.], 
  DirichletCondition[{p[t, x, y] == 1 - .96 sf[t], v[t, x, y] == 0.}, 
   x == -2.], 
  DirichletCondition[{\[Rho][t, x, y] == 1, p[t, x, y] == 1}, 
   x == 2.]};

pde = op == Subscript[\[CapitalGamma], B]; mesh = 
 ToElementMesh[\[CapitalOmega], "MaxCellMeasure" -> 0.004];
mesh["Wireframe"]
Monitor[AbsoluteTiming[{U, V, rho, Psol} = 
    NDSolveValue[{pde, Subscript[\[CapitalGamma], D], ic}, {u, 
      v, \[Rho], p}, {x, y} \[Element] mesh, {t, 0, 7.83}, 
     EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}]), 
     Method -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement", 
         "InterpolationOrder" -> {u -> 2, v -> 2, \[Rho] -> 2, 
           p -> 1}, 
         "InitializePDECoefficientsOptions" -> {"VerificationData" -> \
{"DependentVariables" -> {1, 1, 1, 1}}}}}];], currentTime]

With[{t = 7.83}, {Show[
   DensityPlot[Norm[{U[t, x, y], V[t, x, y]}], {x, y} \[Element] mesh,
     ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
    AspectRatio -> Automatic], 
   StreamPlot[{U[t, x, y], V[t, x, y]}, {x, y} \[Element] mesh, 
    StreamStyle -> LightGray, StreamPoints -> Fine, 
    AspectRatio -> Automatic]],
  DensityPlot[Psol[t, x, y], {x, y} \[Element] mesh, 
   ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
   AspectRatio -> Automatic],
  DensityPlot[rho[t, x, y], {x, y} \[Element] mesh, 
   ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
   AspectRatio -> Automatic]}]

Fig. 4 shows the flow velocity and pressure, which is in agreement with Fig. 1 (there should not be complete coincidence, since different boundary conditions)

Figure 4

Finally, we use the author code. The following corrections should be made there: 1) the first equation should be of the form (the author made three typos - lost the second term, used $\rho$ instead of $p$ and changed the sign in front of the Laplacian) $$\frac {\partial \vec {v}}{\partial t}+(\vec {v}.\nabla) \vec{v}+\frac {\nabla p}{\rho}-\nu \nabla^2\vec {v}-\frac {\nu}{3}\nabla (\nabla .\vec {v})=0$$

here $\nu=\mu/\rho$; 2) the second equation should be of the form $$\frac {\partial \rho}{\partial t}+\nabla .(\rho \vec {v})=0$$ the author has lost $\vec {v}.\nabla \rho$; 3) the third equation should have the form $$\frac {Dp}{dt}-\gamma \frac {p}{\rho}\frac {D\rho}{Dt}-(\gamma-1)\mu \chi -(\gamma-1) \kappa \nabla ^2(\frac {P}{\rho})=0, \frac {D}{Dt}=\frac {\partial }{\partial t}+\vec {v}.\nabla$$ the author has lost $\vec {v}.\nabla$ and $\mu $; 4) the last equation should be of the form $p=\rho RT$, therefore, we can not use the equation for Τ[t,x,y]. Now the code looks like this (I converted the third equation using the second):

<< NDSolve`FEM`

ifctn = Interpolation[{{2.01, 0.6}, {1.5, 0.6}, {0.5, 0.35}, {-1.0, 
     1.0}, {-2.01, 1.2}}];
\[CapitalOmega] = 
  ImplicitRegion[
   Abs[y] <= 
    Abs[ifctn[Max[-2.01, Min[2.01, x]]]], {{x, -2, 2}, {y, -1.2, 
     1.2}}];

softrampfunction[min_, max_, c_, r_] := 
 Function[t, 
  Max[0., (min*Exp[c*r] + max*Exp[r*t])/(Exp[c*r] + Exp[r*t])]]
sf = softrampfunction[0, 1, 3, 3];

af = softrampfunction[0, 1, 3, 4];

\[Chi][t_, x_, y_] := 4/3 (D[v[t,x,y],y]^2-D[v[t,x,y],y] D[u[t,x,y],x]+D[u[t,x,y],x]^2)+(D[u[t,x,y],y]+D[v[t,x,y],x])^2;
\[Lambda] = 1.4; Pr = .71; Ra = 
 6*10^3; k = (\[Lambda] - 1); k1 = \[Lambda];
parameters = {\[Mu] -> Sqrt[Pr/Ra], \[Epsilon] -> 1/Sqrt[Pr*Ra], 
   Subscript[R, m] -> 1 - 1/\[Lambda], 
   Subscript[C, v] -> 1/\[Lambda], \[Kappa] -> Sqrt[Pr/Ra]/Pr};

op = Flatten[{D[{u[t, x, y], v[t, x, y]}, {t, 
        1}] + {{u[t, x, y], v[t, x, y]}.Grad[
         u[t, x, y], {x, y}], {u[t, x, y], v[t, x, y]}.Grad[
         v[t, x, y], {x, y}]} + 
      Grad[p[t, x, y], {x, y}]/\[Rho][t, x, 
        y] - \[Mu] (Div[
          Grad[{u[t, x, y], v[t, x, y]}, {x, y}], {x, y}] + 
         1/3 Grad[Div[{u[t, x, y], v[t, x, y]}, {x, y}], {x, y}]), 
     D[\[Rho][t, x, y], t] + 
      Div[\[Rho][t, x, y]*{u[t, x, y], v[t, x, y]}, {x, y}], 
     D[p[t, x, y], {t, 1}] + {u[t, x, y], v[t, x, y]}.Grad[
        p[t, x, y], {x, y}] + 
      k1 p[t, x, y] Div[{u[t, x, y], v[t, x, y]}, {x, y}] - 
      k (\[Mu] \[Chi][t, x, y] + \[Kappa]*
          Laplacian[p[t, x, y]/\[Rho][t, x, y], {x, y}])}] /. 
   parameters;
ic = {u[0, x, y] == 0, v[0, x, y] == 0, \[Rho][0, x, y] == 1, 
   p[0, x, y] == 1};

Subscript[\[CapitalGamma], B] = {0.0, 0.0, 0.0, 
  0.0};(*NeumannValue[-100sf[t],(x\[GreaterEqual]2)&&(Abs[y]\
\[LessEqual]0.2)]*)
Subscript[\[CapitalGamma], 
  D] = {DirichletCondition[{u[t, x, y] == 0., v[t, x, y] == 0.}, -2. <
     x < 2.], 
  DirichletCondition[{v[t, x, y] == 0., \[Rho][t, x, y] == 
     1 - sf[t] .9, p[t, x, y] == (1 - sf[t] .9)^k1}, x == -2.], 
  DirichletCondition[{\[Rho][t, x, y] == 1, v[t, x, y] == 0., 
    p[t, x, y] == 1}, x == 2.]};

pde = op == Subscript[\[CapitalGamma], B]; mesh = 
 ToElementMesh[\[CapitalOmega], "MaxCellMeasure" -> 0.004];
mesh["Wireframe"]

Monitor[AbsoluteTiming[{U, V, rho, Psol} = 
    NDSolveValue[{pde, Subscript[\[CapitalGamma], D], ic}, {u, 
      v, \[Rho], p}, {x, y} \[Element] mesh, {t, 0, 8}, 
     EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}]), 
     Method -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement", 
         "InterpolationOrder" -> {u -> 2, v -> 2, \[Rho] -> 1, 
           p -> 2}, 
         "InitializePDECoefficientsOptions" -> {"VerificationData" -> \
{"DependentVariables" -> {1, 1, 1, 1}}}}}];], currentTime]
With[{t = 8.}, {Show[
   DensityPlot[Norm[{U[t, x, y], V[t, x, y]}], {x, y} \[Element] mesh,
     ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
    AspectRatio -> Automatic], 
   StreamPlot[{U[t, x, y], V[t, x, y]}, {x, y} \[Element] mesh, 
    StreamStyle -> LightGray, StreamPoints -> Fine, 
    AspectRatio -> Automatic]],
  DensityPlot[Psol[t, x, y], {x, y} \[Element] mesh, PlotRange -> All,
    ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
   AspectRatio -> Automatic],
  DensityPlot[rho[t, x, y], {x, y} \[Element] mesh, PlotRange -> All, 
   ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
   AspectRatio -> Automatic],
  Plot[{rho[t, x, 0], Psol[t, x, 0]}, {x, -2, 2}]}]

Surprisingly, this code works very well. The result is consistent with the first code and the second Figure 5

Unfortunately, the $\chi$ function is not calculated correctly by the author. According to equation (1.74) by http://farside.ph.utexas.edu/teaching/336L/Fluidhtml/node14.html, $$\chi = \frac{\partial v_i}{\partial x_j}\frac{\partial v_i}{\partial x_j}+ \frac{\partial v_i}{\partial x_j}\frac{\partial v_j}{\partial x_i}-\frac {2}{3} \frac{\partial v_i}{\partial x_i}\frac{\partial v_j}{\partial x_j}$$ The corresponding code for this function is

\[Chi][t_, x_, y_] := 
 With[{V = {u[t, x, y], v[t, x, y]}, X = {x, y}}, 
  Sum[D[V[[i]], X[[j]]] D[V[[i]], X[[j]]] + 
     D[V[[i]], X[[j]]] D[V[[j]], X[[i]]] - 
     2/3 D[V[[i]], X[[i]]] D[V[[j]], X[[j]]], {i, 1, 2}, {j, 1, 2}] //
    FullSimplify]

I wrote the correct expression into the code.

$\endgroup$
3
$\begingroup$

The issue here is that currently NDSolve with the finite element method can not handle non-constant coefficients for time derivatives. In order to work around that I reformulated your equations by dividing the equations by that non-constant coefficient. You need to check that I did not make a mistake here.

I started by simplifying your D calls in the equations.

Other than that I fixed the options such that they will to a time integration - specifying Method->"FiniteElement" will force NDSolve to solve this as a time independent PDE. Also, I changed the InterpolationOrder request. I think you want the pressure on a lower order not the density.

When we do this, we run into another issue: NDSolve needs to verify that the coefficients given to it are of the correct dimensions, evaluate to numerical values, etc. This is not entirely trivial for the nonlinaer case and it makes a bad choice in this specific case. To work around that I have changed the VerificationData for the dependent variables, because the choice resulted in a division by zero.

<< NDSolve`FEM`

ifctn = Interpolation[{{2.01, 0.6}, {1.5, 0.6}, {0.5, 0.35}, {-1.0, 
     1.0}, {-2.01, 1.2}}];
\[CapitalOmega] = 
  ImplicitRegion[
   Abs[y] <= 
    Abs[ifctn[Max[-2.01, Min[2.01, x]]]], {{x, -2, 2}, {y, -1.2, 
     1.2}}];

softrampfunction[min_, max_, c_, r_] := 
 Function[t, 
  Max[0., (min*Exp[c*r] + max*Exp[r*t])/(Exp[c*r] + Exp[r*t])]]
sf = softrampfunction[0, 1, 3, 3];

af = softrampfunction[0, 1, 3, 4];
aperaturefunction[w_] := 
 Function[y, af[(10*y/w) + 5.5] - af[(10*y/w) + 0.5]]
aperature = aperaturefunction[1.2];

\[Chi][t_, x_, y_] := 
  D[u[t, x, y], {t, 0}, {x, 1}, {y, 0}]*
    D[u[t, x, y], {t, 0}, {x, 1}, {y, 0}] + 
   D[u[t, x, y], {t, 0}, {x, 0}, {y, 1}]*
    D[v[t, x, y], {t, 0}, {x, 1}, {y, 0}] - 
   2/3 D[u[t, x, y], {t, 0}, {x, 1}, {y, 0}]*
    D[v[t, x, y], {t, 0}, {x, 0}, {y, 1}] + 
   D[v[t, x, y], {t, 0}, {x, 1}, {y, 0}]*
    D[v[t, x, y], {t, 0}, {x, 1}, {y, 0}] + 
   D[v[t, x, y], {t, 0}, {x, 1}, {y, 0}]*
    D[u[t, x, y], {t, 0}, {x, 0}, {y, 1}] - 
   2/3 D[v[t, x, y], {t, 0}, {x, 0}, {y, 1}]*
    D[u[t, x, y], {t, 0}, {x, 1}, {y, 0}];

parameters = {\[Mu] -> Sqrt[Pr/Ra], \[Epsilon] -> 
     1/Sqrt[Pr*Ra], \[CapitalMu] -> M, Subscript[R, m] -> Ra/M, 
    Subscript[C, v] -> 1, \[Lambda] -> 1.2, \[Kappa] -> 1} /. {M -> 
     18.01528, Pr -> 7.1, Ra -> 2*10^5};
parameters0 = {\[Mu] -> 1, \[CapitalMu] -> 1, Subscript[R, m] -> 1, 
   Subscript[C, v] -> 10, \[DoubledGamma] -> 1.4, \[Kappa] -> 1};
op = Flatten[
    {D[{u[t, x, y], v[t, x, y]}, {t, 1}] + 
      Grad[\[Rho][t, x, y], {x, y}]/\[Rho][t, x, 
        y] + \[Mu]/\[Rho][t, x, 
         y] (Div[Grad[{u[t, x, y], v[t, x, y]}, {x, y}], {x, y}] + 
         1/3 Grad[Div[{u[t, x, y], v[t, x, y]}, {x, y}], {x, y}]),

     D[\[Rho][t, x, y], {t, 1}] + \[Rho][t, x, y]*
       Div[{u[t, x, y], v[t, x, y]}, {x, y}], 
     1/(\[Lambda] - 1)*(D[p[t, x, y], {t, 1}] - \[Lambda]*
          Subscript[R, m]*\[CapitalTau][t, x, 
           y] D[\[Rho][t, x, y], {t, 1}]) - \[Chi][t, x, 
       y] - \[Kappa] Inactive[Laplacian][\[CapitalTau][t, x, y], {x, 
         y}],

     D[\[CapitalTau][t, x, y], {t, 1}] - 
      1/\[Rho][t, x, y]*
       Subscript[C, v]/
        M (Subscript[R, m]*\[Rho][t, x, y]*\[CapitalTau][t, x, y]*
          Div[{u[t, x, y], v[t, x, y]}, {x, y}] + \[Chi][t, x, 
          y] + \[Kappa]*
          Inactive[Laplacian][\[CapitalTau][t, x, y], {x, y}])}
    ] /. parameters;
ic = {u[0, x, y] == 0, v[0, x, y] == 0, \[Rho][0, x, y] == 1, 
   p[0, x, y] == 0, \[CapitalTau][0, x, y] == 0};

Subscript[\[CapitalGamma], B] = {0.0, 0.0, 0.0, 0.0, 
  0.0};(*NeumannValue[-100sf[t],(x\[GreaterEqual]2)&&(Abs[y]\
\[LessEqual]0.2)]*)
Subscript[\[CapitalGamma], 
  D] = {DirichletCondition[{u[t, x, y] == 0., 
    v[t, x, y] == 0.}, (-2. < x < 2.) && (y^2 >= ifctn[x]^2)], 
  DirichletCondition[u[t, x, y] == sf[t]*aperature[y], (x >= 2.)], 
  DirichletCondition[v[t, x, y] == 0., (x <= -2.)], 
  DirichletCondition[p[t, x, y] == 0., (x <= -2.)], 
  DirichletCondition[\[CapitalTau][t, x, y] == 1000., (x >= 2.)]};

pde = op == Subscript[\[CapitalGamma], B];
Monitor[AbsoluteTiming[{Subscript[U, vel], Subscript[V, vel], 
     Subscript[\[Rho], density], Subscript[P, ressure], 
     Subscript[T, emperatue]} = 
    NDSolveValue[{pde, Subscript[\[CapitalGamma], D], ic}, {u, 
      v, \[Rho], p, \[CapitalTau]}, {x, y} \[Element] 
      ToElementMesh[\[CapitalOmega], MaxCellMeasure -> 0.005], {t, 0, 
      10}, EvaluationMonitor :> (currentTime = 
        Row[{"t = ", CForm[t]}]), 
     Method -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement", 
         "InterpolationOrder" -> {u -> 2, v -> 2, \[Rho] -> 2, 
           p -> 1, \[CapitalTau] -> 2}, 
         "InitializePDECoefficientsOptions" -> {"VerificationData" -> \
{"DependentVariables" -> {1, 1, 1, 1, 1}}}}}];], currentTime]

This will now start to time integrate but hit a stiff spot at t=0.125. Most likely this is an issue in the equations. You will have to look at that.

One last thing. If you make updates to your post, you'd need to announce that in some way: Some people use a bold Updates or you can ping people in a comment by using @username. Otherwise it will be hard to know when you made improvements to your code.

$\endgroup$
  • $\begingroup$ Unfortunately, even after fixing all the errors, this code does not work. $\endgroup$ – Alex Trounev Oct 9 '19 at 19:40
  • $\begingroup$ @AlexTrounev, Works as I have written (will hit the stiff spot) with Version 12.0 - what errors do you see? $\endgroup$ – user21 Oct 10 '19 at 5:37
  • $\begingroup$ OK! After fixing all the errors in op, this code works even up to t=10(it takes more time than my code). But I have not debug the test case yet. After all, I will add debugged code to my post. $\endgroup$ – Alex Trounev Oct 10 '19 at 11:55
  • $\begingroup$ See update to my answer. $\endgroup$ – Alex Trounev Oct 10 '19 at 23:17
  • $\begingroup$ @AlexTrounev, great! $\endgroup$ – user21 Oct 11 '19 at 4:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.