FEM: periodic solution of 2D Navier-Stokes equations

Let’s consider a horizontal channel with a round obstacle in the middle.

rules = {length -> 2, hight -> 1/2};
Ω =
RegionDifference[Rectangle[{0, 0}, {length, hight}],
Disk[{1, 1/4}, 1/15]] /. rules;
region = RegionPlot[Ω, AspectRatio -> Automatic]


The flow occurs under the action of horizontal force from initial state at rest. Side boundaries are open on which the periodical condition is specified.

op = {
Derivative[1, 0, 0][u][t, x, y] +
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}] +
Derivative[0, 1, 0][p][t, x, y] + (1 - Exp[-t]),
Derivative[1, 0, 0][v][t, x, y] +
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}] + Derivative[0, 0, 1][p][t, x, y],
Derivative[0, 1, 0][u][t, x, y] + Derivative[0, 0, 1][v][t, x, y]};
ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0};
bcsp = {
PeriodicBoundaryCondition[u[t, x, y],
x == 0 && 0 < y < hight, TranslationTransform[{length, 0}]],
PeriodicBoundaryCondition[v[t, x, y],
x == 0 && 0 < y < hight, TranslationTransform[{length, 0}]],
DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0},
0 < x < length],
DirichletCondition[p[t, x, y] == 0.,
x == length && y == hight]} /. rules;
Monitor[AbsoluteTiming[{xVel1, yVel1, pressure1} =
NDSolveValue[{op == {0, 0, 0}, bcsp, ic}, {u, v,
p}, {x, y} ∈ Ω, {t, 0, 5},
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.0005},
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}}},
EvaluationMonitor :> (currentTime =
Row[{"t = ", CForm[t]}])];], currentTime]


Flow is not appeared. It stays trivial.

StreamDensityPlot[
Evaluate[{xVel1[5, x, y], yVel1[5, x, y]}], {x,
y} ∈ Ω, ColorFunction -> "Rainbow",
PlotLegends -> Placed[Automatic, Top], AspectRatio -> Automatic,
ImageSize -> 800, PlotRange -> All]


The problem is that the pressure doen't know that side walls are open. Is it possible to overcome it?

ContourPlot[
Evaluate[pressure1[5, x, y]], {x, y} ∈ Ω,
ColorFunction -> "Rainbow", PlotLegends -> Placed[Automatic, Top],
AspectRatio -> Automatic, ImageSize -> 800, PlotRange -> All]


Is it possible to overcome it and obtain something like that?

Let me explain why answer below by Alex Trounev does not completely satisfy. This preriodic solution (velocity and pressure fields) can be appeared in an infinite periodic array of same obstacles. Let's take 5 for example

rules = {length -> 2.5, hight -> 1/2};
Ω =
RegionDifference[Rectangle[{0, 0}, {length, hight}],
RegionUnion[Table[Disk[{n 1/2 - 1/4, 1/4}, 1/15], {n, 5}]]] /.
rules;
region = RegionPlot[Ω, AspectRatio -> Automatic]


Let's solve the problem with periodic-like boundary conditions as suggested by Alex Trounev

op = {Derivative[1, 0, 0][u][t, x, y] +
10^-2 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}] +
Derivative[0, 1, 0][p][t, x, y] - (1 - Exp[-t]),
Derivative[1, 0, 0][v][t, x, y] +
10^-2 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}] + Derivative[0, 0, 1][p][t, x, y],
Derivative[0, 1, 0][u][t, x, y] + Derivative[0, 0, 1][v][t, x, y]};
ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0};
bcsp = {PeriodicBoundaryCondition[u[t, x, y], x == 0 && 0 < y < hight,
TranslationTransform[{length, 0}]],
PeriodicBoundaryCondition[v[t, x, y], x == 0 && 0 < y < hight,
TranslationTransform[{length, 0}]],
DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0},
0 < x < length], DirichletCondition[p[t, x, y] == 0., True]} /.
rules;
Monitor[AbsoluteTiming[{xVel1, yVel1, pressure1} =
NDSolveValue[{op == {0, 0, 0}, bcsp, ic}, {u, v,
p}, {x, y} ∈ Ω, {t, 0, 10},
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.0005},
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}}},
EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])];],
currentTime];


We obtain almost steady velocity and pressure fields

As you can see that periodic pressure appers at the middle (y-profiles are the same) and not at the side walls.

Plot[{pressure1[5, 1, y], pressure1[5, 1.5, y], pressure1[5, 0, y],
pressure1[5, 2.5, y]}, {y, 0, 1/2}, PlotRange -> All,
PlotStyle -> {Black, {Red, Dashed, Thick}, {Blue, Dashed, Thick},
Green}, Frame -> True,
PlotLegends -> {"x=1", "x=1.5", "x=0", "x=2.5"}]


Note true periodic pressure depends on y. So DirichletCondition[p[t, x, y] == 0., True] is partial solution because of periodicity is slightly broken near side walls. Any other suggestions are welcome.

• In this problem there is no symmetry between left and right side. Why you try to get periodic solution? Apr 29, 2020 at 12:44
• This problem is superposition of symmetry and antisymmetry. In fact it doesn't matter. The obstacle can be moved in center of symmetry. Apr 29, 2020 at 14:00
• You are wrong, just have a look on your picture "something like that". There is no symmetry at all. See also my solution on community.wolfram.com/groups/-/m/t/1433064 Apr 29, 2020 at 14:11
• I've updated my post in order the geometry will be symmetric. But it does not help. Apr 29, 2020 at 15:29
• When you say the side boundaries are open are you saying that you are trying to model the flow in an infinite periodic array of circular cylinders? Also, are you trying to have a moving cylinder or are you pushing flow past the cylinder?
– Hugh
Apr 29, 2020 at 15:42

1 Answer

There is periodic solution with zero pressure drop:

Needs["NDSolveFEM"]
rules = {length -> 2, hight -> 1/2}; reg1 = Disk[{1, 1/4}, 1/15];
reg = RegionDifference[Rectangle[{0, 0}, {length, hight}], reg1] /.
rules;
region = RegionPlot[reg, AspectRatio -> Automatic]

mesh = ToElementMesh[reg, AccuracyGoal -> 5, PrecisionGoal -> 5,
"MaxCellMeasure" -> 0.0005, "MaxBoundaryCellMeasure" -> 0.01]

mesh["Wireframe"]

op = {\[Rho]*D[u[t, x, y], t] +
Inactive[Div][-\[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] - (1 - Exp[-t]), \[Rho]*D[v[t, x, y], t] +
Inactive[Div][-\[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};

tInit = 0; {L, H} = {2, .5};
ic = {u[tInit, x, y] == 0, v[tInit, x, y] == 0, p[tInit, x, y] == 0};
bcs = {DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0},
0 < x < L], DirichletCondition[p[t, x, y] == 0., True]};
bcsp = {PeriodicBoundaryCondition[u[t, x, y], x == 0 && 0 < y < H,
TranslationTransform[{L, 0}]],
PeriodicBoundaryCondition[v[t, x, y], x == 0 && 0 < y < H,
TranslationTransform[{L, 0}]]};
Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[{xVel, yVel, pressure} =
NDSolveValue[{op == {0, 0, 0}, bcs, bcsp, ic}, {u, v,
p}, {x, y} \[Element] mesh, {t, tInit, 1},
Method -> {"TimeIntegration" -> {"IDA",
"MaxDifferenceOrder" -> 2},
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}}},
EvaluationMonitor :> (currentTime = t;)];]


Visualisation

Show[{ContourPlot[xVel[1, x, y], {x, y} \[Element] mesh,
ColorFunction -> "Rainbow",
PlotLegends -> Placed[Automatic, Bottom], AspectRatio -> Automatic,
Contours -> 20, PlotRange -> All, ImageSize -> 400],
StreamPlot[
Evaluate[{xVel[1, x, y], yVel[1, x, y]}], {x, y} \[Element] mesh,
StreamStyle -> LightGray, AspectRatio -> Automatic]}]
ContourPlot[pressure[1, x, y], {x, y} \[Element] mesh,
ColorFunction -> "Rainbow", AspectRatio -> Automatic, Contours -> 20,
PlotRange -> All, PlotLegends -> Placed[Automatic, Bottom],
PlotPoints -> 100]


• Comments are not for extended discussion; this conversation has been moved to chat. May 1, 2020 at 11:12
• This solution sutisfies u(x=0,y)=u(x=L,y) ,v(x=0,y)=v(x=L,y), p(x=0,y)=p(x=L,y) which is not suffient to be periodic. If function is periodic f(x)=f(x+L) then f'(x)=f'(x+L). So equality of partial derivatives normal to the side boundaries is very much in demand. May 1, 2020 at 11:58