Message Boards Message Boards

GROUPS:

Solver for unsteady flow with the use of Mathematica FEM

Posted 3 months ago
1194 Views
|
7 Replies
|
12 Total Likes
|

fig7

I started the discussion here but I also want to repeat on this forum. 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: https://mathematica.stackexchange.com/questions/94914/nonlinear-fem-solver-for-navier-stokes-equations-in-2d/96579#96579

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

fig8

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, lift coefficient and pressure difference in the frontal and aft part of the cylinder as functions of time, maximum drag coefficient, maximum lift coefficient , Strouhal number 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 U0[y_, t_] := 4*Um*y/H*(1 - y/H)*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, for example, test 2D2

{3.17805, 1.03297, 0.266606, 2.60427} v11.01
{3.15711, 1.11377, 0.266043, 2.54356} v11.03

The attached file contains the working code for test 2D3 describing the flow around the cylinder in a flat channel with a change in the flow velocity.

7 Replies

The Euler method is a poor way to solve a differential equation. It can be numerically unstable.

We almost passed the test, therefore, the solution is stable. Explicit Euler is widely used in solving problems in hydrodynamics. The stability of the Euler method depends on the step of integration. In this case, the difference in solutions in v11.01 and v11.3 is due to integration over the surface of the cylinder when calculating the drag coefficient and lift coefficient. In our example, 2000 integrals over the surface are calculated. From these data, the maximum values of the coefficients and the Strouhal number are calculated. This is probably the main source of errors.

enter image description here - Congratulations! This post is now a Staff Pick as distinguished by a badge on your profile! Thank you, keep it coming!

Would be helpful to see the mathematical equations and conditions in traditional textbook format. Also would be helpful to see the animation of the fluid from t=0…100s

I express my gratitude to the Moderation Team for appreciating my work. I'll add an animation for the 2D2 test here and the Navier-Stokes for an unsteady flow with an explanation of the statement of the problem for tests 2D2 and 2D3 from the cited article of M. Schäfer and S. Turek

fig7

fig8

Thanks for it! It is pretty amazing and fun to see that some, for example textbook-scope, fluid mechanics example simulations can be carried out in Mathematica, i.e. by way of pure coding and then hitting the Shift+Enter button.

Needless to comment that for fluid mechanics simulations in any slightly more complicated 2-D geometries or RL 3-D geometries, one would resort to a GUI-driven software where the user interacts with the mouse and context menus and does not have to code anything. Often enough I miss GUI interactivity (similar to Word, Excel. PowerPoint, Comsol) in our software, but then I also realize that Mathematica grants more flexibility and direct control over other details which can't be manipulated in those specialized user-friendly applications.

The bottom line is that for each task at hand one should choose the proper tool wisely. Me too, i enjoy trying to use the software for all and everything at least once. Good for my coding learning experience.

I tested this method of solution on the problems of natural convection. Here we give three tests for the problem of natural convection of air in a rectangular cavity with a Rayleigh number of $Ra=10^3,10^4, 10^5$. On the side walls of the cavity a constant temperature difference is maintained, the upper and lower walls are thermally insulated. We use the Boussinesq approximation. As a standard, we use the data from the article:

Vahl Davis, G.de (1983) : Natural convection of air in a square cavity : A bench mark numerical solution. Int. J. Numer. Methods Fluids 3, 249-264.

Coincidence in all three tests is generally good. Thus, we showed the application of the method for solving the problems of thermal convection in the case of laminar flows. A system of equations describing free convection

$\nabla .\vec u=0, \frac{d\vec u}{dt}=Pr\nabla ^2\vec u -RaPrT\frac {\vec g}{g}$

$\frac {dT}{dt}=\nabla^2T$

$ d/dt =\partial/\partial t +(\vec u .\nabla )$

$\vec u$ is velocity field vector, $T$ - temperature, $\vec g $ - gravity vector, $Pr$ - the Prandtl number, $Ra$ - the Rayleigh number.

Code for the case $Pr=0.71, Ra=10^3$:

K = 100; Pr0 = .71; Ra = 10^3; R = Ra*Pr0; t0 = 1/100; a = 0;
\[CapitalOmega] = ImplicitRegion[0 <= x <= 1 && 0 <= y <= 1, {x, y}];;
UX[0][x_, y_] := 0;
VY[0][x_, y_] := 0;
\[CapitalRho][0][x_, y_] := 0;
TK[0][x_, y_] := 0; TX[0][x_, y_] := 0;
Do[
  {UX[i], VY[i], \[CapitalRho][i], TK[i], TX[i]} = 
   NDSolveValue[{{Inactive[
           Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
             u[x, y], {x, y}]), {x, y}] + D[p[x, y], x] + 
         UX[i - 1][x, y]*D[u[x, y], x] + 
         VY[i - 1][x, y]*D[u[x, y], y] + (u[x, y] - UX[i - 1][x, y])/
          t0 - R*T[x, y]*Sin[a], 
        Inactive[
           Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
             v[x, y], {x, y}]), {x, y}] + D[p[x, y], y] + 
         UX[i - 1][x, y]*D[v[x, y], x] + 
         VY[i - 1][x, y]*D[v[x, y], y] - 
         R*T[x, y]*Cos[a] + (v[x, y] - VY[i - 1][x, y])/t0, 
        D[u[x, y], x] + D[v[x, y], y]} == {0, 0, 0} /. \[Mu] -> Pr0, {
      DirichletCondition[{p[x, y] == 0}, y == 1 && x == 1], 
      DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, 
       y == 0 || y == 1]}, {UX[i - 1][x, y]*D[T[x, y], x] + 
        VY[i - 1][x, y]*D[T[x, y], y] + (T[x, y] - TK[i - 1][x, y])/
         t0 - (D[T[x, y], x, x] + D[T[x, y], y, y]) == 
       NeumannValue[0, y == 0 || y == 1], 
      DirichletCondition[{u[x, y] == 0, v[x, y] == 0, T[x, y] == 1}, 
       x == 0], 
      DirichletCondition[{u[x, y] == 0, v[x, y] == 0, T[x, y] == 0}, 
       x == 1]}}, {u, v, p, T, 
\!\(\*SuperscriptBox[\(T\), 
TagBox[
RowBox[{"(", 
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)}, {x, y} \[Element] \[CapitalOmega], 
    Method -> {"FiniteElement", 
      "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1, T -> 2}, 
      "MeshOptions" -> {"MaxCellMeasure" -> 0.001}}], {i, 1, K}];
{ContourPlot[TK[K][x, y], {x, y} \[Element] \[CapitalOmega], 
PlotLegends -> Automatic, Contours -> 20, PlotPoints -> 25, 
ColorFunction -> "TemperatureMap", AxesLabel -> {"x", "y"}, 
PlotLabel -> T], 
StreamDensityPlot[{UX[K][x, y], 
VY[K][x, y]}, {x, y} \[Element] \[CapitalOmega], 
StreamPoints -> Fine, StreamStyle -> LightGray, 
PlotLegends -> Automatic, VectorPoints -> Fine, 
ColorFunction -> "TemperatureMap", PlotLabel -> Row[{"Ra=", Ra}]]}

MaximalBy[Table[{x, VY[K][x, .5]}, {x, 0, 1, .0001}], Last]

 {{0.1787, 3.69762}}

(*BenchmarkCase={x=0.178,vmax=3.697}*)

 MaximalBy[Table[{y, UX[K][0.5, y]}, {y, 0, 1, .0001}], Last]

 {{0.8132, 3.64937}}

 (*BenchmarkCase={y=0.813,umax=3.649}*)



 MaximalBy[Table[{y, -TX[K][0, y]}, {y, 0, 1, 0.0001}], Last]

 {{0.0909, 1.50724}}

 (*BenchmarkCase={y=0.092,Numax=1.505}*)

 NuAv = 
 NIntegrate[
  UX[K][x, y]*TK[K][x, y] - TX[K][x, y], {x, 0, 1}, {y, 0, 1}, 
  WorkingPrecision -> 5]

 1.1159

 (*BenchmarkCase=1.118*)


 Nu0 =  NIntegrate[-TX[K][0, y], {y, 0, 1}, WorkingPrecision -> 5]

 1.1182

(*BenchmarkCase=1.117*)

 Nu05 = 
 NIntegrate[UX[K][.5, y]*TK[K][.5, y] - TX[K][0.5, y], {y, 0, 1}, 
  WorkingPrecision -> 5]

1.1179

 (*BenchmarkCase=1.118*)

MinimalBy[Table[{y, -TX[K][0, y]}, {y, 0, 1, 0.0001}], Last]

 {{1., 0.691182}}

 (*BenchmarkCase={y=1,Numin=0.692}*)

fig1

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract