# FEM: how to find 2D flow under external forcing?

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


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


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


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


• 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]}; Apr 25, 2020 at 22:04
• @RodionStepanov Yes, I know, see update to my answer. Apr 25, 2020 at 22:44
• 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`? Apr 27, 2020 at 13:56
• @RodionStepanov There is no periodic boundary condition for pressure in your conditions. So what kind of channel flow you suppose to simulate? Apr 27, 2020 at 16:11
• @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. Apr 28, 2020 at 12:11