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$
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

| improve this answer | |
$\endgroup$
  • 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$ – Rodion Stepanov Apr 25 at 22:04
  • 2
    $\begingroup$ @RodionStepanov Yes, I know, see update to my answer. $\endgroup$ – Alex Trounev Apr 25 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$ – Rodion Stepanov Apr 27 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$ – Alex Trounev Apr 27 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$ – Alex Trounev Apr 28 at 12:11

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.