# Solver for unsteady flow with the use of the Navier-Stokes and Mathematica FEM

There are many commercial and open code for solving the problems of unsteady flows. We are interested in the possibility of solving these problems using Mathematica FEM. Previously proposed solvers for stationary incompressible isothermal flows:

Solving 2D Incompressible Flows using Finite Elements http://community.wolfram.com/groups/-/m/t/610335

FEM Solver for Navier-Stokes equations in 2D http://community.wolfram.com/groups/-/m/t/611304

Nonlinear FEM Solver for Navier-Stokes equations in 2D Nonlinear FEM Solver for Navier-Stokes equations in 2D

We give several examples of the successful application of the finite element method for solving unsteady problem including nonisothermal and compressible flows. We will begin with two standard tests that were proposed to solve this class of problems by

M. Schäfer and S. Turek, Benchmark computations of laminar ﬂow around a cylinder (With support by F. Durst, E. Krause and R. Rannacher). In E. Hirschel, editor, Flow Simulation with High-Performance Computers II. DFG priority research program results 1993-1995, number 52 in Notes Numer. Fluid Mech., pp.547–566. Vieweg, Weisbaden, 1996. https://www.uio.no/studier/emner/matnat/math/MEK4300/v14/undervisningsmateriale/schaeferturek1996.pdf

Let us consider the flow in a flat channel around a cylinder at Reynolds number = 100, when self-oscillations occur leading to the detachment of vortices in the aft part of cylinder. In this problem it is necessary to calculate drag coeﬃcient $$c_D$$, lift coeﬃcient $$c_L$$ and pressure diﬀerence $$\Delta P$$ in the frontal and aft part of the cylinder as functions of time, maximum drag coeﬃcient $$c_{Dmax}$$, maximum lift coeﬃcient $$c_{Lmax}$$, Strouhal number $$St$$ and pressure diﬀerence $$\Delta P(t)$$ at $$t = t0 +1/2f$$. The frequency $$f$$ is determined by the period of oscillations of lift coeﬃcient $$f=f(c_L)$$. The data for this test, the code and the results are shown below.

H = .41; L = 2.2; {x0, y0, r0} = {1/5, 1/5, 1/20};
Ω = RegionDifference[Rectangle[{0, 0}, {L, H}], Disk[{x0, y0}, r0]];
RegionPlot[Ω, AspectRatio -> Automatic]

K = 2000; Um = 1.5; ν = 10^-3; t0 = .004;
U0[y_, t_] := 4*Um*y/H*(1 - y/H)
UX[x_, y_] := 0;
VY[x_, y_] := 0;
P0[x_, y_] := 0;
Do[
{UX[i], VY[i], P0[i]} =
NDSolveValue[{{Inactive[
u[x, y], {x, y}]), {x, y}] +
\!$$\*SuperscriptBox[\(p$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] + (u[x, y] - UX[i - 1][x, y])/t0 +
UX[i - 1][x, y]*D[u[x, y], x] +
VY[i - 1][x, y]*D[u[x, y], y],
Inactive[
v[x, y], {x, y}]), {x, y}] +
\!$$\*SuperscriptBox[\(p$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] + (v[x, y] - VY[i - 1][x, y])/t0 +
UX[i - 1][x, y]*D[v[x, y], x] +
VY[i - 1][x, y]*D[v[x, y], y],
\!$$\*SuperscriptBox[\(u$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] +
\!$$\*SuperscriptBox[\(v$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y]} == {0, 0, 0} /. μ -> ν, {
DirichletCondition[{u[x, y] == U0[y, i*t0], v[x, y] == 0},
x == 0.],
DirichletCondition[{u[x, y] == 0., v[x, y] == 0.},
0 <= x <= L && y == 0 || y == H],
DirichletCondition[{u[x, y] == 0,
v[x, y] == 0}, (x - x0)^2 + (y - y0)^2 == r0^2],
DirichletCondition[p[x, y] == P0[i - 1][x, y], x == L]}}, {u, v,
p}, {x, y} ∈ Ω,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},
"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}], {i, 1, K}];
{ContourPlot[UX[K/2][x, y], {x, y} ∈ Ω,
AspectRatio -> Automatic, ColorFunction -> "BlueGreenYellow",
FrameLabel -> {x, y}, PlotLegends -> Automatic, Contours -> 20,
PlotPoints -> 25, PlotLabel -> u, MaxRecursion -> 2],
ContourPlot[VY[K/2][x, y], {x, y} ∈ Ω,
AspectRatio -> Automatic, ColorFunction -> "BlueGreenYellow",
FrameLabel -> {x, y}, PlotLegends -> Automatic, Contours -> 20,
PlotPoints -> 25, PlotLabel -> v, MaxRecursion -> 2,
PlotRange -> All]} // Quiet
{DensityPlot[UX[K][x, y], {x, y} ∈ Ω,
AspectRatio -> Automatic, ColorFunction -> "BlueGreenYellow",
FrameLabel -> {x, y}, PlotLegends -> Automatic, PlotPoints -> 25,
PlotLabel -> u, MaxRecursion -> 2],
DensityPlot[VY[K][x, y], {x, y} ∈ Ω,
AspectRatio -> Automatic, ColorFunction -> "BlueGreenYellow",
FrameLabel -> {x, y}, PlotLegends -> Automatic, PlotPoints -> 25,
PlotLabel -> v, MaxRecursion -> 2, PlotRange -> All]} // Quiet
dPl = Interpolation[
Table[{i*t0, (P0[i][.15, .2] - P0[i][.25, .2])}, {i, 0, K, 1}]];

cD = Table[{t0*i, NIntegrate[(-ν*(-Sin[θ] (Sin[θ]
\!$$\*SuperscriptBox[\(UX[i]$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]] + Cos[θ]
\!$$\*SuperscriptBox[\(UX[i]$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]]) + Cos[θ] (Sin[θ]
\!$$\*SuperscriptBox[\(VY[i]$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]] + Cos[θ]
\!$$\*SuperscriptBox[\(VY[i]$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]]))*Sin[θ] -
P0[i][x0 + r Cos[θ], y0 + r Sin[θ]]*
Cos[θ]) /. {r -> r0}, {θ, 0, 2*Pi}]}, {i,
1000, 2000}]; // Quiet

cL = Table[{t0*i, -NIntegrate[(-ν*(-Sin[θ] (Sin[θ]
\!$$\*SuperscriptBox[\(UX[i]$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]] + Cos[θ]
\!$$\*SuperscriptBox[\(UX[i]$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]]) +
Cos[θ] (Sin[θ]
\!$$\*SuperscriptBox[\(VY[i]$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]] + Cos[θ]
\!$$\*SuperscriptBox[\(VY[i]$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]]))*Cos[θ] +
P0[i][x0 + r Cos[θ], y0 + r Sin[θ]]*
Sin[θ]) /. {r -> r0}, {θ, 0, 2*Pi}]}, {i,
1000, 2000}]; // Quiet
{ListLinePlot[cD,
AxesLabel -> {"t", "\!$$\*SubscriptBox[\(c$$, $$D$$]\)"}],
ListLinePlot[cL,
AxesLabel -> {"t", "\!$$\*SubscriptBox[\(c$$, $$L$$]\)"}],
Plot[dPl[x], {x, 0, 8}, AxesLabel -> {"t", "ΔP"}]}

f002 = FindFit[cL, a*.5 + b*.8*Sin[k*16*t + c*1.], {a, b, k, c}, t]

Plot[Evaluate[a*.5 + b*.8*Sin[k*16*t + c*1.] /. f002], {t, 4, 8},
Epilog -> Map[Point, cL]]

k0=k/.f002;
Struhalnumber = .1*16*k0/2/Pi

cLm = MaximalBy[cL, Last]

sol = {Max[cD[[All, 2]]], Max[cL[[All, 2]]], Struhalnumber,
dPl[cLm[[1, 1]] + Pi/(16*k0)]}


In Fig. 1 shows the components of the flow velocity and the required coefficients. Our solution of the problem and what is required in the test

{3.17805, 1.03297, 0.266606, 2.60427}

lowerbound= { 3.2200, 0.9900, 0.2950, 2.4600};
upperbound = {3.2400, 1.0100, 0.3050, 2.5000}; Note that our results differ from allowable by several percent, but if you look at all the results of Table 4 from the cited article, then the agreement is quite acceptable.The worst prediction is for the Strouhal number. We note that we use the explicit Euler method, which gives an underestimate of the Strouhal number, as follows from the data in Table 4.
The next test differs from the previous one in that the input speed varies according to the Sin[Pi*t/8]. It is necessary to determine the time dependence of the drag and lift parameters for a half-period of oscillation, as well as the pressure drop at the last moment of time. In Fig. 2 shows the components of the flow velocity and the required coefficients. Our solution of the problem and what is required in the test

sol = {3.0438934441256595,
0.5073345082785012, -0.11152933279750943};

lowerbound = {2.9300, 0.4700, -0.1150};
upperbound = {2.9700, 0.4900, -0.1050}; For this test, the agreement with the data in Table 5 is good. Consequently, the two tests are almost completely passed.

I wrote and debugged this code using Mathematics 11.01. But when I ran this code using Mathematics 11.3, I got strange pictures, for example, the disk is represented as a hexagon, the size of the area is changed. In addition, the numerical solution of the problem has changed

{3.17805, 1.03297, 0.266606, 2.60427} v11.0.1
{3.15711, 1.11377, 0.266043, 2.54356} v11.3.0


Update 1 After the release of version 12, I began testing linear and non-linear FEM on the first problem. The good news is that linear FEM is working, and the result is slightly different from versions 11.01 and 11.3, but in general the test is passed. I will give for comparison the output data for the test 2D2 in three versions of the linear FEM

{3.17805, 1.03297, 0.266606, 2.60427} v11.0.1
{3.15711, 1.11377, 0.266043, 2.54356} v11.3.0
{3.15711, 1.11377, 0.266043, 2.54356} v12.0.0.0


I also managed to find the parameters for nonlinear FEM to run the 2D2 test. In this code, I used GaussianQuadratureWeights[] to calculate 600 integrals.

Needs["NDSolveFEM"];
Get["NumericalDifferentialEquationAnalysis"];
np = 60; points = weights = Table[Null, {np}]; Do[
points[[i]] = GaussianQuadratureWeights[np, 0, 2*Pi][[i, 1]], {i, 1,
np}]
Do[weights[[i]] = GaussianQuadratureWeights[np, 0, 2*Pi][[i, 2]], {i,
1, np}]
GaussInt[f_, z_] :=
Sum[(f /. z -> points[[i]])*weights[[i]], {i, 1, np}]
{x0, y0, r0, \[Nu]} = {1/5, 1/5, 1/20, 10^-3};
rules = {length -> 22/10, hight -> 41/100};
\[CapitalOmega] =
RegionDifference[Rectangle[{0, 0}, {length, hight}],
Disk[{1/5, 1/5}, 1/20]] /. rules;
region = RegionPlot[\[CapitalOmega], AspectRatio -> Automatic]

op = {
\[Rho]*D[u[t, x, y], t] +
Inactive[
u[t, x, y], {x, y}]), {x,
y}] + \[Rho] *{{u[t, x, y], v[t, x, y]}}.Inactive[Grad][
u[t, x, y], {x, y}] + D[p[t, x, y], x],
\[Rho]*D[v[t, x, y], t] +
Inactive[
v[t, x, y], {x, y}]), {x,
y}] + \[Rho] *{{u[t, x, y], v[t, x, y]}}.Inactive[Grad][
v[t, x, y], {x, y}] + D[p[t, x, y], y],
D[u[t, x, y], x] + D[v[t, x, y], y]} /. {\[Mu] -> 10^-3, \[Rho] ->
1};

bcs = {DirichletCondition[
u[t, x, y] == If[t < 10^-4, 0, 4*1.5*y*(hight - y)/hight^2],
x == 0], DirichletCondition[u[t, x, y] == 0., 0 < x < length],
DirichletCondition[v[t, x, y] == 0, 0 <= x < length],
DirichletCondition[p[t, x, y] == 0., x == length]} /. rules;
ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0};
Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[{xVel, yVel, pressure} =
NDSolveValue[{op == {0, 0, 0}, bcs, ic}, {u, v,
p}, {x, y} \[Element] \[CapitalOmega], {t, 0, 8},
Method -> {
"TimeIntegration" -> {"IDA", "MaxDifferenceOrder" -> 2},
"PDEDiscretization" -> {"MethodOfLines",
"DifferentiateBoundaryConditions" -> True,
"SpatialDiscretization" -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},
"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}}},
StartingStepSize -> 5.*10^-11,
EvaluationMonitor :> (currentTime = t;)];]


Calculating integrals

cD = Table[{t, GaussInt[(-\[Nu]*(-Sin[\[Theta]] (Sin[\[Theta]]
\!$$\*SuperscriptBox[\(xVel$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x0 + r Cos[\[Theta]],
y0 + r Sin[\[Theta]]] + Cos[\[Theta]]
\!$$\*SuperscriptBox[\(xVel$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x0 + r Cos[\[Theta]],
y0 + r Sin[\[Theta]]]) + Cos[\[Theta]] (Sin[\[Theta]]
\!$$\*SuperscriptBox[\(yVel$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x0 + r Cos[\[Theta]],
y0 + r Sin[\[Theta]]] + Cos[\[Theta]]
\!$$\*SuperscriptBox[\(yVel$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x0 + r Cos[\[Theta]],
y0 + r Sin[\[Theta]]]))*Sin[\[Theta]] -
pressure[t, x0 + r Cos[\[Theta]], y0 + r Sin[\[Theta]]]*
Cos[\[Theta]]) /. {r -> r0}, \[Theta]]}, {t, 5, 8, .01}];

cL = Table[{t, -GaussInt[(-\[Nu]*(-Sin[\[Theta]] (Sin[\[Theta]]
\!$$\*SuperscriptBox[\(xVel$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x0 + r Cos[\[Theta]],
y0 + r Sin[\[Theta]]] + Cos[\[Theta]]
\!$$\*SuperscriptBox[\(xVel$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x0 + r Cos[\[Theta]],
y0 + r Sin[\[Theta]]]) + Cos[\[Theta]] (Sin[\[Theta]]
\!$$\*SuperscriptBox[\(yVel$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x0 + r Cos[\[Theta]],
y0 + r Sin[\[Theta]]] + Cos[\[Theta]]
\!$$\*SuperscriptBox[\(yVel$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x0 + r Cos[\[Theta]],
y0 + r Sin[\[Theta]]]))*Cos[\[Theta]] +
pressure[t, x0 + r Cos[\[Theta]], y0 + r Sin[\[Theta]]]*
Sin[\[Theta]]) /. {r -> r0}, \[Theta]]}, {t, 5, 8, .01}];

dPl = Interpolation[
Table[{t, (pressure[t, .15, .2] - pressure[t, .25, .2])}, {t, 0,
8, .05}]];
{ListLinePlot[cD,
AxesLabel -> {"t", "\!$$\*SubscriptBox[\(c$$, $$D$$]\)"}],
ListLinePlot[cL,
AxesLabel -> {"t", "\!$$\*SubscriptBox[\(c$$, $$L$$]\)"}],
Plot[dPl[x], {x, 0, 8}, AxesLabel -> {"t", "\[CapitalDelta]P"}]}
f002 = FindFit[cL[[190 ;; 290]],
a*.5 + b*.8*Sin[k*16*t + c*1.], {a, b, k, c}, t]
Plot[Evaluate[a*.5 + b*.8*Sin[k*16*t + c*1.] /. f002], {t, 4, 8},
Epilog -> Map[Point, cL]]
Struhalnumber = .1*16*k/2/Pi /. f002
cLm = MaximalBy[cL[[190 ;; 290]], Last]
sol = {Max[cD[[All, 2]]], Max[cL[[190 ;; 290, 2]]], Struhalnumber,
dPl[cLm[[1, 1]] + Pi/(16*k) /. f002]}


As a result, we find

{3.30573, 1.49414, 0.284643, 2.54809}


This is different from the boundaries defined by other method

lowerbound= { 3.2200, 0.9900, 0.2950, 2.4600};
upperbound = {3.2400, 1.0100, 0.3050, 2.5000};


The greatest error (about 50%) has the lift coefficient.

• Very interesting. But people will wonder where the question is. Remember that this is a Q&A site and thus it works a bit different from Wolfram Community. It is perfectly legitimate, though, to pose a question and self-answer it afterwards. Aug 31, 2018 at 7:10
• The question about the code itself - to test it and, perhaps, improve. Aug 31, 2018 at 7:27
• Moreover, the Do loop enforces several other computations to be over and over again (e.g. sparsity pattern, LU-factorization of the linear systems etc.). Actually, one does not want to use LU-factorization for solving the linear systems; being a transient pde, iterative solver such as GMRES should perform much better. Well, I guess we have to wait for upcoming versions of Mathematica to handle nonlinear PDE like this with FEM correctly and efficiently. Aug 31, 2018 at 8:04
• @Henrik Schumacher Thank you. I understand the limitations of this approach, but there are always questions about the solver of problems for unsteady flows using Mathematica. Aug 31, 2018 at 8:31
• @user21 Yes, I'm trying, but so far without success. Apr 22, 2019 at 11:13

Here is a result for your second test (update 2). Let me show how I tackled example c) 2D-3 (unsteady), page 4 from the cited paper. The reason I choose this example is that for example 2D-2 the initial conditions described in the text are not clear to me; they do not seem to be u=v=p=0.
$HistoryLength = 0; Needs["NDSolveFEM"] H = 41/100; L = 22/10; {x0, y0, r0} = {1/5, 1/5, 1/20}; Ω = RegionDifference[Rectangle[{0, 0}, {L, H}], Disk[{x0, y0}, r0]];  We use this as a mesh where we use include points for the later pressure measurement. We refine the boundary a bit and uses a bit finer base mesh as the default. mesh = ToElementMesh[Ω, "IncludePoints" -> {{0.15, 0.2}, {0.25, 0.2}}, AccuracyGoal -> 5, PrecisionGoal -> 5, "MaxCellMeasure" -> 0.0005, "MaxBoundaryCellMeasure" -> 0.01] mesh["Wireframe"] We set up the equation: op = { ρ*D[u[t, x, y], t] + Inactive[Div][-μ Inactive[Grad][u[t, x, y], {x, y}], {x, y}] + ρ *{{u[t, x, y], v[t, x, y]}}.Inactive[Grad][ u[t, x, y], {x, y}] + D[p[t, x, y], x], ρ*D[v[t, x, y], t] + Inactive[Div][-μ Inactive[Grad][v[t, x, y], {x, y}], {x, y}] + ρ *{{u[t, x, y], v[t, x, y]}}.Inactive[Grad][ v[t, x, y], {x, y}] + D[p[t, x, y], y], D[u[t, x, y], x] + D[v[t, x, y], y]} /. {μ -> 10^-3, ρ -> 1};  Set up the inflow profile u0[Um_][x_, y_] := 4*Um*y (H - y)/H^2 uMax = 1.5; uMean = 2*u0[uMax][0, H/2]/3;  And a function that will ramp up the inflow. rampFunction[min_, max_, c_, r_] := Function[t, (min*Exp[c*r] + max*Exp[r*t])/(Exp[c*r] + Exp[r*t])]  The function gets a minimum value (0) a maximum value (1) and the point in time when the ramp should be at 1/2 (the value -2) and a value for the 'steepness' (5) sf = rampFunction[0, 1, -2, 5]; Plot[sf[t], {t, -10, 10}, PlotRange -> All] The boundary conditions: bcs = {DirichletCondition[{u[t, x, y] == sf[t]*4*1.5*Sin[Pi*t/8]*y*(H - y)/H^2, v[t, x, y] == 0}, x == 0], DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, 0 < x < L], DirichletCondition[p[t, x, y] == 0., x == L]};  What we are doing is starting the time integration at tInit=-6 and let the solver move to 0 smoothly: tInit = -6; ic = {u[tInit, x, y] == 0, v[tInit, x, y] == 0, p[tInit, x, y] == 0};  Note that at t=0 the inflow velocity will be 0 in the x-direction due to the Sin. Do the integration: Dynamic["time: " <> ToString[CForm[currentTime]]] AbsoluteTiming[{xVel, yVel, pressure} = NDSolveValue[{op == {0, 0, 0}, bcs, ic}, {u, v, p}, {x, y} ∈ mesh, {t, tInit, 8}, Method -> { "TimeIntegration" -> {"IDA", "MaxDifferenceOrder" -> 2}, "PDEDiscretization" -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}}}, EvaluationMonitor :> (currentTime = t;)];]  This takes about 2 minutes and then we compute the requested pressure difference and compare with the expected lower and upper bounds for that value. dP = pressure[8, 0.15, 0.2] - pressure[8, 0.25, 0.2] -0.115 <= (dP) <= 0.1050 -0.10925764169799679 True  I have not looked into drag and lift coefficients but those should be doable as well; I did do them for a stationary problem were the worked out. Hope this gets you started. Here are some pictures for the pressure: frames = Table[ ContourPlot[pressure[t, x, y], {x, y} ∈ mesh, Contours -> 20, PlotRange -> All, AspectRatio -> All], {t, 0, 8, 0.25}]; ListAnimate[frames] • @Nasser, thanks for the cleanup. May 1, 2019 at 5:03 • Thanks for this example. But this is not a test 2D3, it is a completely different problem. However, with this example we see the limitations of the nonlinear FEM. It looks like when solving problems in which there is a clear dependence of boundary conditions on time, we must think about how to change the boundary conditions so that the solution converges. May 1, 2019 at 11:34 • @AlexTrounev, could you be specific which example you want solved (page number and example name) May 2, 2019 at 12:33 • You used the function rampFunction`, so the test conditions are not met. However, your code allows us to correctly determine the parameters$\Delta P, c_D$, but gives an error of 50% when determining$c_L\$. May 2, 2019 at 13:27