7
$\begingroup$

I'm trying to sovle a 2D Navier-Stokes equation in a rectangular region with no-flow boundary conditions and an external rotating force $f(x,y)=(-y,x)$. I've got the error message:

NDSolveValue::index: The DAE solver failed at t = 0.. The solver is intended for index 1 DAE systems and structural analysis indicates that the DAE index is 2. The option Method->{"IndexReduction"->Automatic} may be used to reduce the index of the system.

I don't see why it does not work. It very similar to the tutorial examples. Could you, please, help me to solve this problem?

This is my code:

Ω = Rectangle[{-1, -1}, {1, 1}];

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] - y, 
   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] + x, 
   Derivative[0, 1, 0][u][t, x, y] + Derivative[0, 0, 1][v][t, x, y]};

bcs = 
  {DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, True], 
   DirichletCondition[p[t, x, y] == 0, x == -1 && y == -1]};

ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0};

Monitor[
  AbsoluteTiming[
    {xVel, yVel, pressure} = 
      NDSolveValue[
        {op == {0, 0, 0}, bcs, ic}, {u, v, p}, {x, y} ∈ Ω, {t, 0, 1},
        Method -> 
          {"PDEDiscretization" -> 
            {"MethodOfLines",
             "SpatialDiscretization" -> 
               {"FiniteElement", 
                "MeshOptions" -> {"MaxCellMeasure" -> 0.01}, 
                "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}}}, 
        EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])];], 
  currentTime]
```
$\endgroup$

1 Answer 1

9
$\begingroup$

This solver works well at low acceleration, so we need to multiply the force by a damping factor $f=(-y,x)(1-e^{-t})$, then we have a solution

Ω = Rectangle[{-1, -1}, {1, 1}];

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] - 
    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] + 
    x (1 - Exp[-t]), 
   Derivative[0, 1, 0][u][t, x, y] + Derivative[0, 0, 1][v][t, x, y]};

bcs = {DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, True], 
   DirichletCondition[p[t, x, y] == 0, x == -1 && y == -1]};

ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0};

Monitor[AbsoluteTiming[{xVel, yVel, pressure} = 
    NDSolveValue[{op == {0, 0, 0}, bcs, ic}, {u, v, 
      p}, {x, y} ∈ Ω, {t, 0, 10}, 
     Method -> {"PDEDiscretization" -> {"MethodOfLines", 
         "SpatialDiscretization" -> {"FiniteElement", 
           "MeshOptions" -> {"MaxCellMeasure" -> 0.01}, 
           "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}}}, 
     EvaluationMonitor :> (currentTime = 
        Row[{"t = ", CForm[t]}])];], currentTime]

Table[StreamDensityPlot[
  Evaluate[{xVel[t, x, y], yVel[t, x, y]}], {x, 
    y} ∈ Ω, 
  PlotLabel -> Row[{"t = ", t}]], {t, 1, 10, 3}]

Figure 1 Update 1. Numerical solution looks quite different with periodic boundary condition on x:

bcsp = {PeriodicBoundaryCondition[u[

t, x, y], x == -1 && -1 < y < 1, 
    TranslationTransform[{2, 0}]], 
   PeriodicBoundaryCondition[v[t, x, y], x == -1 && -1 < y < 1, 
    TranslationTransform[{2, 0}]], 
   DirichletCondition[{u[t, x, y] == 0, 
     v[t, x, y] == 0}, (y == -1 || y == 1) && -1 < x < 1], 
   DirichletCondition[p[t, x, y] == 0, x == -1 && y == -1]};
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.01}, 
           "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}}}, 
     EvaluationMonitor :> (currentTime = 
        Row[{"t = ", CForm[t]}])];], currentTime]
Table[StreamDensityPlot[
  Evaluate[{xVel1[t, x, y], yVel1[t, x, y]}], {x, 
    y} ∈ Ω, PlotLabel -> Row[{"t = ", t}], 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic], {t, 1, 4, 1}]

Figure 2

To simulate Poiseuille flow in a flat channel driven by force $F_x=t$ we put $p=const$ and remove continuity equation from a system, finally we have

op1 = {Derivative[1, 0, 0][u][t, x, y] + 
    Inactive[Div][-I

nactive[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}] - 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}]};

ic1 = {u[0, x, y] == 0, v[0, x, y] == 0};
bcsp1 = {PeriodicBoundaryCondition[u[t, x, y], x == -1 && -1 < y < 1, 
    TranslationTransform[{2, 0}]], 
   PeriodicBoundaryCondition[v[t, x, y], x == -1 && -1 < y < 1, 
    TranslationTransform[{2, 0}]], 
   DirichletCondition[{u[t, x, y] == 0, 
     v[t, x, y] == 0}, (y == -1 || y == 1) && -1 < x < 1]};
Monitor[AbsoluteTiming[{xVel1, yVel1} = 
    NDSolveValue[{op1 == {0, 0}, bcsp1, ic1}, {u, 
      v}, {x, y} \[Element] \[CapitalOmega], {t, 0, 2}, 
     Method -> {"PDEDiscretization" -> {"MethodOfLines", 
         "SpatialDiscretization" -> {"FiniteElement", 
           "MeshOptions" -> {"MaxCellMeasure" -> 0.01}, 
           "InterpolationOrder" -> {u -> 2, v -> 2}}}}, 
     EvaluationMonitor :> (currentTime = 
        Row[{"t = ", CForm[t]}])];], currentTime]

Visualisation

{Table[DensityPlot[
   Evaluate[xVel1[t, x, y]], {x, y} \[Element] \[CapitalOmega], 
   PlotLabel -> Row[{"t = ", t}], ColorFunction -> "Rainbow", 
   PlotLegends -> Automatic], {t, .2, 2, .3}], 
 Plot[Table[xVel1[t, 0, y], {t, .2, 2, .3}], {y, -1, 1}]}

Figure 3

$\endgroup$
12
  • 1
    $\begingroup$ may be you also know how to adjust PeriodicBoundaryCondition? Solver stucks after a few steps. bcs = { PeriodicBoundaryCondition[u[t, x, y], x == -1 && -1 <= y <= 1, TranslationTransform[{2, 0}]], PeriodicBoundaryCondition[v[t, x, y], x == -1 && -1 <= y <= 1, TranslationTransform[{2, 0}]], DirichletCondition[ u[t, x, y] == 0, (y == -1 || y == 1) && -1 < x <= 1], DirichletCondition[ v[t, x, y] == 0, (y == -1 || y == 1) && -1 < x <= 1], DirichletCondition[p[t, x, y] == 0, x == -1 && y == -1]}; $\endgroup$ Apr 25, 2020 at 22:04
  • 2
    $\begingroup$ @RodionStepanov Yes, I know, see update to my answer. $\endgroup$ Apr 25, 2020 at 22:44
  • $\begingroup$ Unfortumatly, channel flow does not appear because of the pressure is still not periodic. Can you get correct channel flow solving equastions with external force Fx=t and Fy=0? $\endgroup$ Apr 27, 2020 at 13:56
  • $\begingroup$ @RodionStepanov There is no periodic boundary condition for pressure in your conditions. So what kind of channel flow you suppose to simulate? $\endgroup$ Apr 27, 2020 at 16:11
  • 1
    $\begingroup$ @RodionStepanov There are no unique solution with this force, since $p=tx, u=0,v=0$ is a solution too. So FEM gets this solution first. $\endgroup$ Apr 28, 2020 at 12:11

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.