# The FEMDampingElements operator failed

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.

<<NDSolveFEM

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[ρ[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]

• 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 . Sep 29, 2019 at 19:49
• Unfortunately, you code has syntax errors when I copy it. If you fix those I can have another look. Ping me with @user21. Sep 30, 2019 at 12:40
• @CarlVoss There are some errors with signs in op. What does this model describe? Oct 10, 2019 at 17:53
• @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) Oct 11, 2019 at 2:45

FEM code for isentropic viscous flows. We 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 with Mach number of about $$M=3.77$$.

<< NDSolveFEM
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[
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[
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]


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]


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]


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

<< NDSolveFEM

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)

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):

<< NDSolveFEM

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

The $$\chi$$ function calculated according to equation (1.74) from this book as follows $$\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]


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 do 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.

<< NDSolveFEM

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.

• Unfortunately, even after fixing all the errors, this code does not work. Oct 9, 2019 at 19:40
• @AlexTrounev, Works as I have written (will hit the stiff spot) with Version 12.0 - what errors do you see? Oct 10, 2019 at 5:37
• 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. Oct 10, 2019 at 11:55
• See update to my answer. Oct 10, 2019 at 23:17
• @AlexTrounev, great! Oct 11, 2019 at 4:25