Steady state fluid flow for downward flow past an obstacle

Question

I was attempting to simulate fluid flow past a circular obstacle. The following is the code which I used

(*define the variables*)
rules3 = {\[Mu] -> 1.63*10^-2, \[Rho] -> 1.1136, g -> -9.81, 
  epr -> 2*10^-3, ExitV -> 1, PeakV -> 1.5, PipePos -> 0}

(*Define the flow region*)
\[CapitalOmega] = 
  RegionUnion[Rectangle[{-epr, 0}, {epr, 0.05}], 
    RegionDifference[Rectangle[{-epr - 10^-3, -0.1}, {0.075, 0}], 
     Disk[{6.5*10^-3, -0.05}, 6*10^-3]]] /. rules3;

(*Define the Navier-Stokes Equation*)
snsef1 = {
   -\[Mu] \!\(
\*SubsuperscriptBox[\(\[Del]\), \({x, y}\), \(2\)]\(u[x, 
       y]\)\) + \[Rho] {u[x, y], v[x, y]} . \!\(
\*SubscriptBox[\(\[Del]\), \({x, y}\)]\(u[x, y]\)\) + \!\(
\*SubscriptBox[\(\[PartialD]\), \(x\)]\(p[x, y]\)\),
   -\[Mu] \!\(
\*SubsuperscriptBox[\(\[Del]\), \({x, y}\), \(2\)]\(v[x, 
       y]\)\) + \[Rho] {u[x, y], v[x, y]} . \!\(
\*SubscriptBox[\(\[Del]\), \({x, y}\)]\(v[x, y]\)\) + \!\(
\*SubscriptBox[\(\[PartialD]\), \(x\)]\(p[x, y]\)\) - \[Rho] g,
   \!\(
\*SubscriptBox[\(\[Del]\), \({x, y}\)] . \({u[x, y], v[x, y]}\)\)
   } /. rules3

(*Modelling the Entry velocity as a Poiseulle Flow*)
InflowBC = 
  DirichletCondition[{u[x, y] == 0, 
     v[x, y] == -PeakV (1 - ((x - PipePos)/epr)^2)}, y == 0.05] /. 
   rules3;

(*Outflow Boundary Condition*)
OutflowBC = DirichletCondition[p[x, y] == 0, y == -0.1];

(*Wall Boundary Condition (No-slip) *)
WallBC = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, 0.05 > y > -0.25];

(*Solving the equation*)
{uVel, vVel, pressure} = NDSolve[{
   snsef1 == {0, 0, 0},
   bcs
   }, {u, v, p}, {x, y} \[Element] \[CapitalOmega], Method -> {
    "PDEDiscretization" -> {"FiniteElement", 
      "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}, 
      "MeshOptions" -> {"MaxCellMeasure" -> 0.0000001},
      "PDESolveOptions" -> {"FindRootOptions" -> {Method -> \
{"AffineCovariantNewton", "MinimalDampingFactor" -> 10^-8}}}}
    }]

This returns the error:

FindRoot::dfmin: The minimal damping factor of 1/100000000 has been reached.

FindRoot::sszero: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the function value is still greater than the tolerance prescribed by the AccuracyGoal option.

NDSolve::fempsf: PDESolve could not find a solution.

Set::shape: Lists {uVel,vVel,pressure} and <<1>> are not the same shape.

What do these errors mean and how can I fix this? Thanks in advance

Clarifications: The Inflow BC represents flow in a pipe, so I used poiseulle flow equation to model this

Edits: The Inflow velocity should be negative (downward), changed in the code

Question Answered:

Thanks to user Alex Trounev for the solution! Appreciate it and thank you!

Answer 1

This problem can be solved with linear FEM as it described here. First, please note, that we can't fixed in one time the inflow, the outflow pressure and force $\rho g$. With a given inflow and pressure we can include potential force in pressure. Second, to resolve flow in a channel and around cylinder, we use scale = 2*10^-3. Third, we use the method of false transient. Finally we have

Needs["NDSolve`FEM`"];

(*define the variables*)scale = 2*10^-3;
rules3 = {\[Mu] -> 1.63*10^-2, \[Rho] -> 1.1136, g -> -9.81, 
     epr -> 1, ExitV -> 1, PeakV -> - 1.5, PipePos -> 0};

(*Define the flow region*)
\[CapitalOmega] = 
    RegionUnion[Rectangle[{-1, 0}, {1, 0.05/scale}], 
        RegionDifference[
     Rectangle[{-epr - 10^-3/scale, -0.1/scale}, {0.075/scale, 0}], 
          Disk[{6.5*10^-3/scale, -0.05/scale}, 6*10^-3/scale]]] /. 
   rules3;


mesh = ToElementMesh[\[CapitalOmega], AccuracyGoal -> 5, 
  PrecisionGoal -> 5, "MaxCellMeasure" -> 0.25, 
  "MaxBoundaryCellMeasure" -> 0.2]

{RegionPlot[\[CapitalOmega], AspectRatio -> Automatic], 
 mesh["Wireframe"]}

UX[0][x_, y_] := 0;
VY[0][x_, y_] := 0;
\[CapitalRho][0][x_, y_] := 0;
(*Modelling the Entry velocity as a Poiseulle Flow*)
InflowBC = 
    DirichletCondition[{u[x, y] == 0, 
          v[x, y] == PeakV (1 - ((x - PipePos))^2)}, 
    y == 0.05/scale] /. 
      rules3;

(*Outflow Boundary Condition*)
OutflowBC = DirichletCondition[p[x, y] == 0, y == -0.1/scale];

(*Wall Boundary Condition (No-slip) *)
WallBC = 
  DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, 
   0.05/scale > y > -0.1/scale];


(*Define the Navier-Stokes Equation*)(*Solving the equation*)
dt = 1/100; Do[{UX[i], VY[i], \[CapitalRho][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])/
         dt, 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] + (v[x, y] - VY[i - 1][x, y])/
         dt, D[u[x, y], x] + D[v[x, y], y]} == {0, 0, 0} /. rules3, 
    InflowBC, OutflowBC, WallBC}, {u, v, p}, Element[{x, y}, mesh], 
   Method -> {"FiniteElement", 
     "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}], {i, 1, 30}];

Note, that 30 iterations we used is not sufficient for the flow stabilization, but velocity variations are very small from iteration to iteration, Flow visualization on 3 last iterations

Table[Show[
  DensityPlot[
   Norm[{UX[i][x, y], VY[i][x, y]}], {x, y} \[Element] mesh, 
   PlotLegends -> Automatic, ColorFunction -> "TemperatureMap", 
   PlotLabel -> Row[{"t =", 1. dt i }], PlotRange -> {0, .2}, 
   AspectRatio -> Automatic], 
  StreamPlot[{UX[i][x, y], VY[i][x, y]}, {x, y} \[Element] mesh, 
    StreamPoints -> Fine, StreamStyle -> LightGray] // Quiet], {i, 28,
   30}]

Detailed visualization flow around cylinder and in the channel entrance

{Show[StreamDensityPlot[{UX[30][x, y], VY[30][x, y]}, {x, -1.5, 
    5}, {y, -5, 1}, PlotLegends -> Automatic, 
   ColorFunction -> "TemperatureMap", PlotRange -> All, 
   AspectRatio -> Automatic], 
  RegionPlot[\[CapitalOmega], AspectRatio -> Automatic, 
   PlotStyle -> Directive[Opacity[0.15], Gray], 
   BoundaryStyle -> Directive[Red, Thick]]], 
 Show[StreamDensityPlot[{UX[30][x, y], VY[30][x, y]}, {x, -1.5, 
    10}, {y, -35, -20}, PlotLegends -> Automatic, 
   ColorFunction -> "TemperatureMap", PlotRange -> All, 
   AspectRatio -> Automatic], 
  RegionPlot[\[CapitalOmega], AspectRatio -> Automatic, 
   PlotStyle -> Directive[Opacity[0.15], Gray], 
   BoundaryStyle -> Directive[Red, Thick]]]}