I started to try to solve PDEs in WM and I go gradually according to tutorial. Unfortunately, I couldn't get any solution other than that copied from the tutorial. So I ask someone more experienced to advise me where I made a mistake. Take, for example, Navier-Sotkes PDEs to solve fluid flow around a cylindrical obstacle (the first part is a stationary solution, the second part is a time-varying case).
My question is general: How to construct initial and boundary conditions for the same time-varying task (second part) so that the initial conditions are not constants but functions and that there is a solution and WM finds it?
My attempt: I tried to set the initial conditions so that at time t = 0, the functions of velocities u and v were equal to the stationary solution. And the boundary condition at one end the value of velocity in the x-direction increased (that is, a variation of the original solution, instead zero initial velocities is a stationary solution) all see the code. However, WM throws an error: "LinearSolve::nosol: Linear equation encountered that has no solution." Do you have any ideas what can be wrong and how to fix it?
CODE:
(*Stationary Navier–Stokes equation*)
ClearAll["Global`*"]
rules = {length -> 22/10, hight -> 41/100};
Ω =
RegionDifference[Rectangle[{0, 0}, {length, hight}],
Disk[{1/5, 1/5}, 1/20]] /. rules;
region = RegionPlot[Ω, AspectRatio -> Automatic]
op = {Inactive[
Div][({{-μ, 0}, {0, -μ}}.Inactive[Grad][
u[x, y], {x, y}]), {x,
y}] + ρ {{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[
Div][({{-μ, 0}, {0, -μ}}.Inactive[Grad][
v[x, y], {x, y}]), {x,
y}] + ρ {{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]} /. {μ -> 10^-3, ρ -> 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[Ω], 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*)
{sxVel, syVel, spressure} =
NDSolveValue[{pde, bcs}, {u, v, p},
Element[{x, y}, Ω],
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[sxVel[x, y]^2 + syVel[x, y]^2], {x,
y} ∈ Ω, ColorFunction -> "TemperatureMap",
Contours -> 4],
ContourPlot[spressure[x, y] == #, {x, y} ∈ Ω,
ContourStyle -> Directive[Black, Thickness[0.002]]] & /@ {0.01,
0.02, 0.05, 0.08},
VectorPlot[{sxVel[x, y],
syVel[x, y]}, {x, y} ∈ Ω, VectorPoints -> 9,
VectorScale -> .05, VectorStyle -> {LightGray}],
AspectRatio -> Automatic, ImageSize -> Large]
(*Transient Navier–Stokes equation*)
ClearAll[ρ, μ]
op = {ρ*D[u[t, x, y], t] +
Inactive[
Div][({{-μ, 0}, {0, -μ}}.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}] +
D[p[t, x, y], x], ρ*D[v[t, x, y], t] +
Inactive[
Div][({{-μ, 0}, {0, -μ}}.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}] + D[p[t, x, y], y],
D[u[t, x, y], x] + D[v[t, x, y], y]} /. {μ -> 10^-3, ρ ->
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;(*original bc*)
bcs = {DirichletCondition[u[t, x, y] == (sf[t] + 1) sxVel[x, y],
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;(*my bc*)
(*initial conditions*)
ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0} /.
rules;(*original ic*)
ic = {u[0, x, y] == sxVel[x, y], v[0, x, y] == syVel[x, y],
p[0, x, y] == spressure[x, y]} /. rules;(*my ic*)
(*solution*)
Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[{xVel, yVel, pressure} =
NDSolveValue[{op == {0, 0, 0}, bcs, ic}, {u, v,
p}, {x, y} ∈ Ω, {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} ∈ 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}];]
bigwig in programming in WMI assumed you meantbeginningthere. Normally when one is starting to learn new language/environment, they pick very small and simple problems first. You seem to jump to the most complex problems right away, and then ask how to solve them using Mathematica when you say yourself you are just a beginner in it. This for me does not sound like a logical way to learn new language. If you want to learn FEM in Mathematica, it is better to start with simpler examples. Just my 2 cents ofcourse as you are free to do what you wish. – Nasser Apr 20 '20 at 04:29