# Is WM able to solve PDEs (non-linear with variable coefficients) with variable initial conditions?

I'm learning how to solve PDE from tutorial. My question is about passage: "The transient Navier–Stokes equation" (code below) at the end of tutorial. First, there is a stable flow through the chanelle and around the cylinder, then a time-dependent flow. I am able to reproduce these results. My questions about time-varying solution:

1) If I change the initial pressure value from $$p(0, x, y) = 0$$ to another constant eg $$p(0, x, y) = 1$$, everything works and the solution is similar. How come there is no problem when the boundary condition says $$p (0, x, y) = 0, x == xend$$, wehre xend is end of region of integration in x coordinate?

2) If I change the initial condition to any function $$p (0, x, y) = x$$, $$p (0, x, y) = \sin(x * pi / xend)$$ - it satisfies the boundary condition $$p (0, x, y ) = 0, x == xend$$, the solution fails (PDE solution time increment decreases until it reaches the critical value about 10 ^ -10). Why does it not work?

3) General question: "A stable solution can be found if the velocities are interpolated with higher order than the pressure." is mentioned when dealing with stationary problem. If I solve a whole new problem, how can I figure it out what to do with higher interpolation? There are a lot of these combinations.

Code:

(*Stationary Navier\[Dash]Stokes equation*)
ClearAll["Global*"]
rules = {length -> 22/10, hight -> 41/100};
\[CapitalOmega] =
RegionDifference[Rectangle[{0, 0}, {length, hight}],
Disk[{1/5, 1/5}, 1/20]] /. rules;
region = RegionPlot[\[CapitalOmega], AspectRatio -> Automatic]

op = {Inactive[
u[x, y], {x, y}]), {x,
y}] + \[Rho] {{u[x, y], v[x, y]}}.Inactive[Grad][
u[x, y], {x, y}] +

\!$$\*SuperscriptBox[\(p$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y],
Inactive[
v[x, y], {x, y}]), {x,
y}] + \[Rho] {{u[x, y], v[x, y]}}.Inactive[Grad][
v[x, y], {x, y}] +

\!$$\*SuperscriptBox[\(p$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y],

\!$$\*SuperscriptBox[\(u$$,
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] +

\!$$\*SuperscriptBox[\(v$$,
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y]} /. {\[Mu] -> 10^-3, \[Rho] -> 1}
pde = op == {0, 0, 0};

(*Boudnary conditions*)
bcs = {DirichletCondition[u[x, y] == 4*0.3*y*(hight - y)/hight^2,
x == 0], DirichletCondition[v[x, y] == 0, x == 0],
DirichletCondition[{u[x, y] == 0., v[x, y] == 0.},
x > 0 && x < length],
DirichletCondition[p[x, y] == 0., x == length]} /. rules;
(*Create a refinement function based on the region of refinement, \
smaller mesh grid in this region - more precises*)
refinementRegion =
ImplicitRegion[
a^4 (23 (y - 1/5))^2 +
b^2 (5/2 x - 2.05)^3 (2 a + (5/2 x - 2.05)) <
0, {{x, 0, length}, {y, 0, hight}}] /.
Flatten[{rules, a -> 1, b -> 5/2}];
Show[RegionPlot[\[CapitalOmega]], RegionPlot[refinementRegion],
AspectRatio -> Automatic]
mrf = With[{rmf = RegionMember[refinementRegion]},
Function[{vertices, area}, Block[{x, y}, {x, y} = Mean[vertices];
If[rmf[{x, y}], area > 0.00025, area > 0.0025]]]];
(*Solution, u has to have higher order of of interpolation than p*)
{xVel, yVel, pressure} =
NDSolveValue[{pde, bcs}, {u, v, p},
Element[{x, y}, \[CapitalOmega]],
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},
"MeshOptions" -> {"IncludePoints" -> {{0.15, 0.2}, {0.25, 0.2}},
AccuracyGoal -> 5, PrecisionGoal -> 5,
MeshRefinementFunction -> mrf}}];
Show[ContourPlot[
Sqrt[xVel[x, y]^2 + yVel[x, y]^2], {x,
y} \[Element] \[CapitalOmega], ColorFunction -> "TemperatureMap",
Contours -> 4],
ContourPlot[pressure[x, y] == #, {x, y} \[Element] \[CapitalOmega],
ContourStyle -> Directive[Black, Thickness[0.002]]] & /@ {0.01,
0.02, 0.05, 0.08},
VectorPlot[{xVel[x, y],
yVel[x, y]}, {x, y} \[Element] \[CapitalOmega], VectorPoints -> 9,
VectorScale -> .05, VectorStyle -> {LightGray}],
AspectRatio -> Automatic, ImageSize -> Large]

(*Transient Navier\[Dash]Stokes equation*)
ClearAll[\[Rho], \[Mu]]
op = {\[Rho]*D[u[t, x, y], t] +
Inactive[
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], \[Rho]*D[v[t, x, y], t] +
Inactive[
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};

(*ramp fucntion*)
rampFunction[min_, max_, c_, r_] :=
Function[t, (min*Exp[c*r] + max*Exp[r*t])/(Exp[c*r] + Exp[r*t])]
sf = rampFunction[0, 1, 4, 5];
Plot[sf[t], {t, -1, 10}, PlotRange -> All]

(*boundary conditions*)
bcs = {DirichletCondition[
u[t, x, y] == sf[t]*4*1.5*y*(hight - y)/hight^2, x == 0],
DirichletCondition[u[t, x, y] == 0., 0 < x < length],
DirichletCondition[v[t, x, y] == 0, 0 <= x < length],
DirichletCondition[p[t, x, y] == 0., x == length]} /. rules;

(*initial conditions*)
ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0};

(*solution*)
Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[{xVel, yVel, pressure} =
NDSolveValue[{op == {0, 0, 0}, bcs, ic}, {u, v,
p}, {x, y} \[Element] \[CapitalOmega], {t, 0, 10},
Method -> {"TimeIntegration" -> {"IDA",
"MaxDifferenceOrder" -> 2},
"PDEDiscretization" -> {"MethodOfLines",
"DifferentiateBoundaryConditions" -> True,
"SpatialDiscretization" -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},
"MeshOptions" -> {"MaxCellMeasure" -> 0.0005}}}},
EvaluationMonitor :> (currentTime = t;)];]

(*visualization*)
{minX, maxX} = MinMax[xVel["ValuesOnGrid"]]
mesh = xVel["ElementMesh"]
AbsoluteTiming[
frames = Table[
Rasterize[
ContourPlot[xVel[t, x, y], {x, y} \[Element] mesh,
PlotRange -> All, AspectRatio -> Automatic,
ColorFunction -> "TemperatureMap",
Contours -> Range[minX, maxX, (maxX - minX)/7], Axes -> False,
Frame -> None], RasterSize -> 2*{360, 68},
ImageSize -> {360, 68}], {t, 0, 10, 1/12}];]
ListAnimate[frames, SaveDefinitions -> True]

• @Nasser Thank you for your comment. I added whole code. I thought my questions are more or less theoretical. I am not sure what you mean by second part. Shall I create "new question" here for any question? Apr 17 '20 at 5:18

Question 1 : pressure is arbitrary for any uniform shift because the gradient of uniform function is zero. There is little relation between the downstream BC and the initial conditions on pressure. They just have to be compatible at time zero. So when you set initial pressure to 1, probably the downstream BC is going to be stronger. Actually, as you also impose the initial speed to zero, then the initial pressure must be uniform (as there is no gravity), hence initial pressure is set everywhere by the boundary condition. Pressure is "special".

Question 2 : you can't choose the initial condition for pressure freely, it should be a solution (in the tutorial : fluid at rest) or maybe not far from a solution. That is why you start from fluid at rest, which has a well known hydrostatic solution.

Question 3 : "A stable solution can be found if the velocities are interpolated with a higher order than the pressure." I can understand this qualitatively. If pressure and velocity could be separated, pressure would obey a Poisson equation and speed would obey a non-linear diffusion equation (Ginzbug-Landau like). Consequently pressure perturbations travel instantly while velocity is on convective or diffusive scales. So pressure can be more dangerous than velocity, hence the greater care for pressure. Also FEM makes local cells and is going to have trouble in capturing the global pressure problem.

• Pierre ALBARÈDE, as you suggested in your answer, I tried to set initial condition as results of steady flow. What do you think? :-) Solution (NDSolve part of the code) blows up with error message: "LinearSolve::nosol: Linear equation encountered that has no solution.". What does it say to me? That initial conditions are not fine? Or that it is impossible to solve such equations with initial conditions defined as a functions? Apr 17 '20 at 12:23

Very interesting but not clear. "initial pressure value from 𝑝(0,𝑥,𝑦)=0" so first coordinate must be time.

Then "the boundary condition says 𝑝(0,𝑥,𝑦)=0,𝑥==𝑒𝑛𝑑". If the first coordinate is time, then your boundary condition is not complete. It should be "𝑝(t,𝑥,𝑦)=0,𝑥==𝑒𝑛𝑑".

The boundary condition must be expressed for all time.

I guess x==end is upstream?

As you only have x,y, I guess you solve Navier-Stokes 2D so you can't study the flow around a cylinder. You just get the flow around some shape in 2D. This problem is quite distinct from the 3D flow around a cylinder.

If so you should write "flow around a disk" to avoid confusion.

• Sorry I did not mean answer but comment. Apr 17 '20 at 7:39
• Pierre ALBARÈDE, Thank you for your comment. All code I presented is coming from tutorial (link in my question), I just copied it and learn from it and try to change. I believe all boundary and initial conditions are defined right. My change is only cosmetic in sense of all calculation (only value of initial condition p nothing more). end (I rewrite it now as xend) is end of integration region. Apr 17 '20 at 7:48
• Pierre ALBARÈDE THANK YOU!! That opened my eyes. Apr 17 '20 at 11:20
• "set initial condition as results of steady flow" that is your very good idea. It should work. Either you or Mathematica goofed. I have also encountered this kind of behaviour of LinearSolve when I know there is a solution. You could try to Monitor the Det. Also you should determine the Reynolds number and keep it low to begin with. Anyway, if you keep asking why on such a problem it can take a life time to answer. Apr 17 '20 at 12:31
• You approximate a non-linear problem with potentially many solutions, not all of them stable, by a linear system that has a unique solution or an infinity of solutions depending on the Det. You can imagine it is not going to be easy. Apr 17 '20 at 12:39

[Sorry this is comment again but does not fit in the box so I put answer again.]

I have had a look at the tutorial and the program. There is inconsistency in your question but maybe not in the program. You do not simulate the flow around a cylinder but the flow inside a channel and around a disk. It is not possible to discuss the problem in terms of height and length. The proper directions are streamwise and crossstream (and spanwise in 3D but you are 2D). There is an upstream BC on speed and downstream BC on pressure and crossstream no-slip (on the side walls). Because of the sidewalls, the upstream boundary condition is parabolic. This is essentially correct. There is stationary and transient cases. In the transient case the BC is time-dependent and the flow is turned on by a step function. Everything looks globally mathematically correct even though it has little relation with real fluid mechanics. Now you have a better stand to reformulate your question about the numerical scheme.

• "A stable solution can be found if the velocities are interpolated with a higher order than the pressure." I can understand this qualitatively. If pressure and velocity could be separated, pressure would obey a Poisson equation and speed would obey a non-linear diffusion equation (Ginzbug-Landau like). Consequently pressure perturbations travel instantly while velocity is on convective or diffusive scales. So pressure can be more dangerous than velocity, hence the greater care for pressure. Also FEM makes local cells and is going to have trouble in capturing the global pressure problem. Apr 17 '20 at 10:35
• Pierre ALBARÈDE, thank you again. I think that your last comment is a bit answer to question 3. On the other hand, my other questions seems to be quite clear. Can Mathematica deal with initial condition as (variable) function? I just tried to help and divided it to two questions and showed a proof it works for other constant. Apr 17 '20 at 10:45
• That is clear now. Pressure and velocity are coupled, You can't set pressure when velocity has been set already. This is quite obvious at rest : if velocity is null then pressure is uniform. In your code u[0, x, y] == 0, v[0, x, y] == 0`. Apr 17 '20 at 11:07
• If vorticity is zero the coupling is Bernoulli theorem p+1/2rvˆ2=uniform. You can see interesting things by applying rot and grad to NS. Apr 17 '20 at 12:16