9
$\begingroup$

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 flow 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 coefficient $c_D$, lift coefficient $c_L$ and pressure difference $\Delta P$ in the frontal and aft part of the cylinder as functions of time, maximum drag coefficient $c_{Dmax}$, maximum lift coefficient $c_{Lmax}$, Strouhal number $St$ and pressure difference $\Delta P(t)$ at $t = t0 +1/2f$. The frequency $f$ is determined by the period of oscillations of lift coefficient $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[0][x_, y_] := 0;
VY[0][x_, y_] := 0;
P0[0][x_, y_] := 0;
Do[
  {UX[i], VY[i], P0[i]} = 
   NDSolveValue[{{Inactive[
           Div][({{-μ, 0}, {0, -μ}}.Inactive[Grad][
             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[
           Div][({{-μ, 0}, {0, -μ}}.Inactive[Grad][
             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};

Fig1

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};

Fig2

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

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["NDSolve`FEM`"];
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[
       Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
         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[
       Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
         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.

$\endgroup$
  • 5
    $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Aug 31 '18 at 7:10
  • $\begingroup$ The question about the code itself - to test it and, perhaps, improve. $\endgroup$ – Alex Trounev Aug 31 '18 at 7:27
  • $\begingroup$ Oh, this runs very slowly; 20 steps need more than 8 seconds. That's actually a well-known issue and is caused by evaluation of InterpolationFunction (i) is implemented somewhat inefficiently and (ii) would always be slower that using the degrees of freedom directly. Moreover, I would have discretized the domain only once with Needs["NDSolveFEM"]; Ω = ToElementMesh[ RegionDifference[Rectangle[{0, 0}, {L, H}], Disk[{x0, y0}, r0]], "MaxCellMeasure" -> 0.001 ]; but apparently, this doesn't make a difference. (There must be a clever caching under the hood.) $\endgroup$ – Henrik Schumacher Aug 31 '18 at 7:55
  • $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Aug 31 '18 at 8:04
  • $\begingroup$ @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. $\endgroup$ – Alex Trounev Aug 31 '18 at 8:31
7
$\begingroup$

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.

To time integrate the transient Navier-Stokes equations one typically needs a special index-2 DAE solver. Currently we do not have that. However, you can trick the DAE solver into solving the equations never the less. For this to work you need to slowly accelerate the inflow. I'll show that in a second.

Set up the problem:

$HistoryLength = 0;
Needs["NDSolve`FEM`"]
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"]

enter image description here

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]

enter image description here

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]

enter image description here

$\endgroup$
  • $\begingroup$ @Nasser, thanks for the cleanup. $\endgroup$ – user21 May 1 at 5:03
  • $\begingroup$ 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. $\endgroup$ – Alex Trounev May 1 at 11:34
  • $\begingroup$ @AlexTrounev, could you be specific which example you want solved (page number and example name) $\endgroup$ – user21 May 2 at 12:33
  • $\begingroup$ 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$. $\endgroup$ – Alex Trounev May 2 at 13:27

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.