Skip to main content

differential equations - Solver for unsteady flow with the use of the Navier-Stokes and Mathematica FEM


There are many commercial and open code for solving the problems of unsteady flows. We are interested in the possibility of solving these problems using Mathematica FEM. Previously proposed solvers for stationary incompressible isothermal flows:


Solving 2D Incompressible Flows using Finite Elements http://community.wolfram.com/groups/-/m/t/610335


FEM Solver for Navier-Stokes equations in 2D http://community.wolfram.com/groups/-/m/t/611304


Nonlinear FEM Solver for Navier-Stokes equations in 2D Nonlinear FEM Solver for Navier-Stokes equations in 2D



We give several examples of the successful application of the finite element method for solving unsteady problem including nonisothermal and compressible flows. We will begin with two standard tests that were proposed to solve this class of problems by


M. Schäfer and S. Turek, Benchmark computations of laminar flow around a cylinder (With support by F. Durst, E. Krause and R. Rannacher). In E. Hirschel, editor, Flow Simulation with High-Performance Computers II. DFG priority research program results 1993-1995, number 52 in Notes Numer. Fluid Mech., pp.547–566. Vieweg, Weisbaden, 1996. https://www.uio.no/studier/emner/matnat/math/MEK4300/v14/undervisningsmateriale/schaeferturek1996.pdf


Let us consider the flow in a flat channel around a cylinder at Reynolds number = 100, when self-oscillations occur leading to the detachment of vortices in the aft part of cylinder. In this problem it is necessary to calculate drag coefficient $c_D$, lift coefficient $c_L$ and pressure difference $\Delta P$ in the frontal and aft part of the cylinder as functions of time, maximum drag coefficient $c_{Dmax}$, maximum lift coefficient $c_{Lmax}$, Strouhal number $St$ and pressure difference $\Delta P(t)$ at $t = t0 +1/2f$. The frequency $f$ is determined by the period of oscillations of lift coefficient $f=f(c_L)$. The data for this test, the code and the results are shown below.


H = .41; L = 2.2; {x0, y0, r0} = {1/5, 1/5, 1/20}; 
Ω = RegionDifference[Rectangle[{0, 0}, {L, H}], Disk[{x0, y0}, r0]];
RegionPlot[Ω, AspectRatio -> Automatic]

K = 2000; Um = 1.5; ν = 10^-3; t0 = .004;
U0[y_, t_] := 4*Um*y/H*(1 - y/H)
UX[0][x_, y_] := 0;

VY[0][x_, y_] := 0;
P0[0][x_, y_] := 0;
Do[
{UX[i], VY[i], P0[i]} =
NDSolveValue[{{Inactive[
Div][({{-μ, 0}, {0, -μ}}.Inactive[Grad][
u[x, y], {x, y}]), {x, y}] +
\!\(\*SuperscriptBox[\(p\),
TagBox[
RowBox[{"(",

RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] + (u[x, y] - UX[i - 1][x, y])/t0 +
UX[i - 1][x, y]*D[u[x, y], x] +
VY[i - 1][x, y]*D[u[x, y], y],
Inactive[
Div][({{-μ, 0}, {0, -μ}}.Inactive[Grad][
v[x, y], {x, y}]), {x, y}] +
\!\(\*SuperscriptBox[\(p\),
TagBox[

RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] + (v[x, y] - VY[i - 1][x, y])/t0 +
UX[i - 1][x, y]*D[v[x, y], x] +
VY[i - 1][x, y]*D[v[x, y], 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]} == {0, 0, 0} /. μ -> ν, {
DirichletCondition[{u[x, y] == U0[y, i*t0], v[x, y] == 0},
x == 0.],

DirichletCondition[{u[x, y] == 0., v[x, y] == 0.},
0 <= x <= L && y == 0 || y == H],
DirichletCondition[{u[x, y] == 0,
v[x, y] == 0}, (x - x0)^2 + (y - y0)^2 == r0^2],
DirichletCondition[p[x, y] == P0[i - 1][x, y], x == L]}}, {u, v,
p}, {x, y} ∈ Ω,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},
"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}], {i, 1, K}];
{ContourPlot[UX[K/2][x, y], {x, y} ∈ Ω,

AspectRatio -> Automatic, ColorFunction -> "BlueGreenYellow",
FrameLabel -> {x, y}, PlotLegends -> Automatic, Contours -> 20,
PlotPoints -> 25, PlotLabel -> u, MaxRecursion -> 2],
ContourPlot[VY[K/2][x, y], {x, y} ∈ Ω,
AspectRatio -> Automatic, ColorFunction -> "BlueGreenYellow",
FrameLabel -> {x, y}, PlotLegends -> Automatic, Contours -> 20,
PlotPoints -> 25, PlotLabel -> v, MaxRecursion -> 2,
PlotRange -> All]} // Quiet
{DensityPlot[UX[K][x, y], {x, y} ∈ Ω,
AspectRatio -> Automatic, ColorFunction -> "BlueGreenYellow",

FrameLabel -> {x, y}, PlotLegends -> Automatic, PlotPoints -> 25,
PlotLabel -> u, MaxRecursion -> 2],
DensityPlot[VY[K][x, y], {x, y} ∈ Ω,
AspectRatio -> Automatic, ColorFunction -> "BlueGreenYellow",
FrameLabel -> {x, y}, PlotLegends -> Automatic, PlotPoints -> 25,
PlotLabel -> v, MaxRecursion -> 2, PlotRange -> All]} // Quiet
dPl = Interpolation[
Table[{i*t0, (P0[i][.15, .2] - P0[i][.25, .2])}, {i, 0, K, 1}]];



cD = Table[{t0*i, NIntegrate[(-ν*(-Sin[θ] (Sin[θ]
\!\(\*SuperscriptBox[\(UX[i]\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]] + Cos[θ]
\!\(\*SuperscriptBox[\(UX[i]\),
TagBox[

RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]]) + Cos[θ] (Sin[θ]
\!\(\*SuperscriptBox[\(VY[i]\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],

MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]] + Cos[θ]
\!\(\*SuperscriptBox[\(VY[i]\),
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]]))*Sin[θ] -
P0[i][x0 + r Cos[θ], y0 + r Sin[θ]]*

Cos[θ]) /. {r -> r0}, {θ, 0, 2*Pi}]}, {i,
1000, 2000}]; // Quiet

cL = Table[{t0*i, -NIntegrate[(-ν*(-Sin[θ] (Sin[θ]
\!\(\*SuperscriptBox[\(UX[i]\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],

y0 + r Sin[θ]] + Cos[θ]
\!\(\*SuperscriptBox[\(UX[i]\),
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]]) +
Cos[θ] (Sin[θ]
\!\(\*SuperscriptBox[\(VY[i]\),

TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]] + Cos[θ]
\!\(\*SuperscriptBox[\(VY[i]\),
TagBox[
RowBox[{"(",
RowBox[{"1", ",", "0"}], ")"}],

Derivative],
MultilineFunction->None]\)[x0 + r Cos[θ],
y0 + r Sin[θ]]))*Cos[θ] +
P0[i][x0 + r Cos[θ], y0 + r Sin[θ]]*
Sin[θ]) /. {r -> r0}, {θ, 0, 2*Pi}]}, {i,
1000, 2000}]; // Quiet
{ListLinePlot[cD,
AxesLabel -> {"t", "\!\(\*SubscriptBox[\(c\), \(D\)]\)"}],
ListLinePlot[cL,
AxesLabel -> {"t", "\!\(\*SubscriptBox[\(c\), \(L\)]\)"}],

Plot[dPl[x], {x, 0, 8}, AxesLabel -> {"t", "ΔP"}]}

f002 = FindFit[cL, a*.5 + b*.8*Sin[k*16*t + c*1.], {a, b, k, c}, t]


Plot[Evaluate[a*.5 + b*.8*Sin[k*16*t + c*1.] /. f002], {t, 4, 8},
Epilog -> Map[Point, cL]]

k0=k/.f002;
Struhalnumber = .1*16*k0/2/Pi



cLm = MaximalBy[cL, Last]

sol = {Max[cD[[All, 2]]], Max[cL[[All, 2]]], Struhalnumber,
dPl[cLm[[1, 1]] + Pi/(16*k0)]}

In Fig. 1 shows the components of the flow velocity and the required coefficients. Our solution of the problem and what is required in the test


{3.17805, 1.03297, 0.266606, 2.60427}


lowerbound= { 3.2200, 0.9900, 0.2950, 2.4600};
upperbound = {3.2400, 1.0100, 0.3050, 2.5000};

Fig1


Note that our results differ from allowable by several percent, but if you look at all the results of Table 4 from the cited article, then the agreement is quite acceptable.The worst prediction is for the Strouhal number. We note that we use the explicit Euler method, which gives an underestimate of the Strouhal number, as follows from the data in Table 4.
The next test differs from the previous one in that the input speed varies according to the Sin[Pi*t/8]. It is necessary to determine the time dependence of the drag and lift parameters for a half-period of oscillation, as well as the pressure drop at the last moment of time. In Fig. 2 shows the components of the flow velocity and the required coefficients. Our solution of the problem and what is required in the test


sol = {3.0438934441256595`, 
0.5073345082785012`, -0.11152933279750943`};

lowerbound = {2.9300, 0.4700, -0.1150};

upperbound = {2.9700, 0.4900, -0.1050};

Fig2


For this test, the agreement with the data in Table 5 is good. Consequently, the two tests are almost completely passed.


I wrote and debugged this code using Mathematics 11.01. But when I ran this code using Mathematics 11.3, I got strange pictures, for example, the disk is represented as a hexagon, the size of the area is changed. Fig3


In addition, the numerical solution of the problem has changed


{3.17805, 1.03297, 0.266606, 2.60427} v11.0.1
{3.15711, 1.11377, 0.266043, 2.54356} v11.3.0

Update 1 After the release of version 12, I began testing linear and non-linear FEM on the first problem. The good news is that linear FEM is working, and the result is slightly different from versions 11.01 and 11.3, but in general the test is passed. I will give for comparison the output data for the test 2D2 in three versions of the linear FEM



{3.17805, 1.03297, 0.266606, 2.60427} v11.0.1
{3.15711, 1.11377, 0.266043, 2.54356} v11.3.0
{3.15711, 1.11377, 0.266043, 2.54356} v12.0.0.0

I also managed to find the parameters for nonlinear FEM to run the 2D2 test. In this code, I used GaussianQuadratureWeights[] to calculate 600 integrals.


Needs["NDSolve`FEM`"];
Get["NumericalDifferentialEquationAnalysis`"];
np = 60; points = weights = Table[Null, {np}]; Do[
points[[i]] = GaussianQuadratureWeights[np, 0, 2*Pi][[i, 1]], {i, 1,
np}]

Do[weights[[i]] = GaussianQuadratureWeights[np, 0, 2*Pi][[i, 2]], {i,
1, np}]
GaussInt[f_, z_] :=
Sum[(f /. z -> points[[i]])*weights[[i]], {i, 1, np}]
{x0, y0, r0, \[Nu]} = {1/5, 1/5, 1/20, 10^-3};
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 = {
\[Rho]*D[u[t, x, y], t] +
Inactive[
Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
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[

Div][({{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
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};

bcs = {DirichletCondition[
u[t, x, y] == If[t < 10^-4, 0, 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;
ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0};
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, 8},
Method -> {
"TimeIntegration" -> {"IDA", "MaxDifferenceOrder" -> 2},
"PDEDiscretization" -> {"MethodOfLines",

"DifferentiateBoundaryConditions" -> True,
"SpatialDiscretization" -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},
"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}}},
StartingStepSize -> 5.*10^-11,
EvaluationMonitor :> (currentTime = t;)];]

Calculating integrals


cD = Table[{t, GaussInt[(-\[Nu]*(-Sin[\[Theta]] (Sin[\[Theta]] 
\!\(\*SuperscriptBox[\(xVel\),

TagBox[
RowBox[{"(",
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x0 + r Cos[\[Theta]],
y0 + r Sin[\[Theta]]] + Cos[\[Theta]]
\!\(\*SuperscriptBox[\(xVel\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1", ",", "0"}], ")"}],

Derivative],
MultilineFunction->None]\)[t, x0 + r Cos[\[Theta]],
y0 + r Sin[\[Theta]]]) + Cos[\[Theta]] (Sin[\[Theta]]
\!\(\*SuperscriptBox[\(yVel\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x0 + r Cos[\[Theta]],
y0 + r Sin[\[Theta]]] + Cos[\[Theta]]

\!\(\*SuperscriptBox[\(yVel\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x0 + r Cos[\[Theta]],
y0 + r Sin[\[Theta]]]))*Sin[\[Theta]] -
pressure[t, x0 + r Cos[\[Theta]], y0 + r Sin[\[Theta]]]*
Cos[\[Theta]]) /. {r -> r0}, \[Theta]]}, {t, 5, 8, .01}];


cL = Table[{t, -GaussInt[(-\[Nu]*(-Sin[\[Theta]] (Sin[\[Theta]]
\!\(\*SuperscriptBox[\(xVel\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x0 + r Cos[\[Theta]],
y0 + r Sin[\[Theta]]] + Cos[\[Theta]]
\!\(\*SuperscriptBox[\(xVel\),
TagBox[

RowBox[{"(",
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x0 + r Cos[\[Theta]],
y0 + r Sin[\[Theta]]]) + Cos[\[Theta]] (Sin[\[Theta]]
\!\(\*SuperscriptBox[\(yVel\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "0", ",", "1"}], ")"}],
Derivative],

MultilineFunction->None]\)[t, x0 + r Cos[\[Theta]],
y0 + r Sin[\[Theta]]] + Cos[\[Theta]]
\!\(\*SuperscriptBox[\(yVel\),
TagBox[
RowBox[{"(",
RowBox[{"0", ",", "1", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[t, x0 + r Cos[\[Theta]],
y0 + r Sin[\[Theta]]]))*Cos[\[Theta]] +
pressure[t, x0 + r Cos[\[Theta]], y0 + r Sin[\[Theta]]]*

Sin[\[Theta]]) /. {r -> r0}, \[Theta]]}, {t, 5, 8, .01}];

dPl = Interpolation[
Table[{t, (pressure[t, .15, .2] - pressure[t, .25, .2])}, {t, 0,
8, .05}]];
{ListLinePlot[cD,
AxesLabel -> {"t", "\!\(\*SubscriptBox[\(c\), \(D\)]\)"}],
ListLinePlot[cL,
AxesLabel -> {"t", "\!\(\*SubscriptBox[\(c\), \(L\)]\)"}],
Plot[dPl[x], {x, 0, 8}, AxesLabel -> {"t", "\[CapitalDelta]P"}]}

f002 = FindFit[cL[[190 ;; 290]],
a*.5 + b*.8*Sin[k*16*t + c*1.], {a, b, k, c}, t]
Plot[Evaluate[a*.5 + b*.8*Sin[k*16*t + c*1.] /. f002], {t, 4, 8},
Epilog -> Map[Point, cL]]
Struhalnumber = .1*16*k/2/Pi /. f002
cLm = MaximalBy[cL[[190 ;; 290]], Last]
sol = {Max[cD[[All, 2]]], Max[cL[[190 ;; 290, 2]]], Struhalnumber,
dPl[cLm[[1, 1]] + Pi/(16*k) /. f002]}

As a result, we find



{3.30573, 1.49414, 0.284643, 2.54809}

This is different from the boundaries defined by other method


lowerbound= { 3.2200, 0.9900, 0.2950, 2.4600};
upperbound = {3.2400, 1.0100, 0.3050, 2.5000};

The greatest error (about 50%) has the lift coefficient.



Answer



Here is a result for your second test (update 2). Let me show how I tackled example c) 2D-3 (unsteady), page 4 from the cited paper. The reason I choose this example is that for example 2D-2 the initial conditions described in the text are not clear to me; they do not seem to be u=v=p=0.


To time integrate the transient Navier-Stokes equations one typically needs a special index-2 DAE solver. Currently we do not have that. However, you can trick the DAE solver into solving the equations never the less. For this to work you need to slowly accelerate the inflow. I'll show that in a second.



Set up the problem:


$HistoryLength = 0;
Needs["NDSolve`FEM`"]
H = 41/100; L = 22/10;
{x0, y0, r0} = {1/5, 1/5, 1/20};
Ω =
RegionDifference[Rectangle[{0, 0}, {L, H}], Disk[{x0, y0}, r0]];

We use this as a mesh where we use include points for the later pressure measurement. We refine the boundary a bit and uses a bit finer base mesh as the default.


mesh = ToElementMesh[Ω, 

"IncludePoints" -> {{0.15, 0.2}, {0.25, 0.2}}, AccuracyGoal -> 5,
PrecisionGoal -> 5, "MaxCellMeasure" -> 0.0005,
"MaxBoundaryCellMeasure" -> 0.01]

mesh["Wireframe"]

enter image description here


We set up the equation:


op = {
ρ*D[u[t, x, y], t] +

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}] + D[p[t, x, y], x],
ρ*D[v[t, x, y], t] +
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}] + D[p[t, x, y], y],
D[u[t, x, y], x] + D[v[t, x, y], y]} /. {μ -> 10^-3, ρ ->
1};


Set up the inflow profile


u0[Um_][x_, y_] := 4*Um*y (H - y)/H^2
uMax = 1.5;
uMean = 2*u0[uMax][0, H/2]/3;

And a function that will ramp up the inflow.


rampFunction[min_, max_, c_, r_] := 
Function[t, (min*Exp[c*r] + max*Exp[r*t])/(Exp[c*r] + Exp[r*t])]

The function gets a minimum value (0) a maximum value (1) and the point in time when the ramp should be at 1/2 (the value -2) and a value for the 'steepness' (5)



sf = rampFunction[0, 1, -2, 5];
Plot[sf[t], {t, -10, 10}, PlotRange -> All]

enter image description here


The boundary conditions:


bcs = {DirichletCondition[{u[t, x, y] == 
sf[t]*4*1.5*Sin[Pi*t/8]*y*(H - y)/H^2, v[t, x, y] == 0}, x == 0],
DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, 0 < x < L],
DirichletCondition[p[t, x, y] == 0., x == L]};


What we are doing is starting the time integration at tInit=-6 and let the solver move to 0 smoothly:


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

Note that at t=0 the inflow velocity will be 0 in the x-direction due to the Sin.


Do the integration:


Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[{xVel, yVel, pressure} =
NDSolveValue[{op == {0, 0, 0}, bcs, ic}, {u, v,
p}, {x, y} ∈ mesh, {t, tInit, 8},

Method -> {
"TimeIntegration" -> {"IDA", "MaxDifferenceOrder" -> 2},
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}}},
EvaluationMonitor :> (currentTime = t;)];]

This takes about 2 minutes and then we compute the requested pressure difference and compare with the expected lower and upper bounds for that value.


dP = pressure[8, 0.15, 0.2] - pressure[8, 0.25, 0.2]
-0.115 <= (dP) <= 0.1050

-0.10925764169799679`
True

I have not looked into drag and lift coefficients but those should be doable as well; I did do them for a stationary problem were the worked out.


Hope this gets you started.


Here are some pictures for the pressure:


frames = Table[
ContourPlot[pressure[t, x, y], {x, y} ∈ mesh,
Contours -> 20, PlotRange -> All, AspectRatio -> All], {t, 0, 8,
0.25}];


ListAnimate[frames]

enter image description here


Comments

Popular posts from this blog

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...