6
$\begingroup$

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]

enter image description here

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]

enter image description here

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]

enter image description here

Is it possible to overcome it and obtain something like that? enter image description here

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]

enter image description here

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 enter image description hereenter image description here

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"}]

enter image description here

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.

$\endgroup$
15
  • 1
    $\begingroup$ In this problem there is no symmetry between left and right side. Why you try to get periodic solution? $\endgroup$ – Alex Trounev Apr 29 '20 at 12:44
  • $\begingroup$ This problem is superposition of symmetry and antisymmetry. In fact it doesn't matter. The obstacle can be moved in center of symmetry. $\endgroup$ – Rodion Stepanov Apr 29 '20 at 14:00
  • $\begingroup$ 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 $\endgroup$ – Alex Trounev Apr 29 '20 at 14:11
  • $\begingroup$ I've updated my post in order the geometry will be symmetric. But it does not help. $\endgroup$ – Rodion Stepanov Apr 29 '20 at 15:29
  • 1
    $\begingroup$ 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? $\endgroup$ – Hugh Apr 29 '20 at 15:42
2
$\begingroup$

There is periodic solution with zero pressure drop:

Needs["NDSolve`FEM`"]
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]

Figure 1

$\endgroup$
2
  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$ – Mr.Wizard May 1 '20 at 11:12
  • $\begingroup$ 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. $\endgroup$ – Rodion Stepanov May 1 '20 at 11:58

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.