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.

Improve this question
$\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
    Commented Apr 29, 2020 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$
    – SantaP
    Commented Apr 29, 2020 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
    Commented Apr 29, 2020 at 14:11
  • $\begingroup$ I've updated my post in order the geometry will be symmetric. But it does not help. $\endgroup$
    – SantaP
    Commented Apr 29, 2020 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
    Commented Apr 29, 2020 at 15:42

1 Answer 1

Reset to default
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

Improve this answer
$\endgroup$
2
  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$
    – Mr.Wizard
    Commented May 1, 2020 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$
    – SantaP
    Commented May 1, 2020 at 11:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.