News, Views & Insights
非線形偏微分方程式への有限要素法の適用
Mathematica 12 has powerful functionality for solving partial differential equations (PDEs) both symbolically and numerically. This article focuses on, among other things, the finite element method (FEM)–based solver for nonlinear PDEs that has been newly implemented in Version 12. After briefly reviewing basic syntax of the Wolfram Language for PDEs, including how to designate Dirichlet and Neumann boundary conditions, we will delineate how Mathematica 12 finds the solution of a given nonlinear problem with FEM. We then show some examples in physics and chemistry, such as the Gray–Scott model and the time-dependent Navier–Stokes equation. More information can be found in the Wolfram Language tutorial “Finite Element Programming,” on which most of this article is based.
1. はじめに
Wolfram Research社の旗艦製品であるMathematicaは，5,000 を超える組み込み関数を有するWolfram Languageを駆動する．数理モデリング，解析の基本となる常・偏微分方程式の分野においては，これらをシンボリックに，あるいは数値的に解くための強力なソルバを搭載している．最近は有限要素法(FEM) を利用した数値的求解機能が大幅に強化され，偏微分方程式(PDE)を任意の領域上で解いたり，固有値・固有関数を求めたりすることが可能となった．ここでは，最新のバージョン12における非線形偏微分方程式のFEMによる求解を中心に，現実的な問題に応用する上での流れを例とともに紹介する．なお，有限要素法を用いて非線形PDEを解くワークフローの詳細，コードはすべて公開されている．MathematicaのWolframドキュメント内で，チュートリアル“FiniteElementProgramming”を参照いただきたい．
2. 微分方程式の数値求解ワークフロー
Wolfram Languageにおいて微分方程式を数値的に解く際の関数(コマンド)は，NDSolve，あるいはNDSolveValueのふたつある．これら二つは出力のフォーマットが若干異なるだけで，中の処理はまったく同じものなので，以下，本文中では表記が短い"NDSolve", コード例中では出力の扱いが簡便な"NDSolveValue" で表記する．Mathematica上でFEMを利用するには，パッケージをロードする :
![Needs["NDSolve`FEM`"]](https://content.wolfram.com/sites/39/2020/04/km0416img1.png "Needs[\"NDSolve`FEM`\"]"){kind=small}
✕
Needs["NDSolve`FEM`"] |
あとは，NDSolveにPDE，領域，初期・境界条件を与えるだけである．たとえば，単位円上のポアソン方程式–∇2u = 1を例にとると，境界条件としてx ≥ 0にある境界でu = 0とすると，
✕
eqn = -Inactive[Div][Inactive[Grad][u[x, y], {x, y}], {x, y}] == 1;
Subscript[\[CapitalOmega], D] = Disk[];
Subscript[\[CapitalGamma], D] =
DirichletCondition[u[x, y] == 0, x >= 0];
usol = NDSolveValue[{eqn, Subscript[\[CapitalGamma], D]},
u, {x, y} \[Element] Subscript[\[CapitalOmega], D]]
|
で解が得られ，
✕
Plot3D[usol[x, y], {x, y} \[Element] Subscript[\[CapitalOmega], D]]
|
とプロットできる．
2.1 式の入力
NDSolveの中の有限要素法が現在適用可能な偏微分方程式は次の形をもたなければならない :
ここで，解くべき従属変数uがn上の1次元関数のときは，m, d, a, fはスカラー，α, γ, βはn次元ベクトル，cはn*n行列である．また，c, a, f, α, γ, βはt, x ∈ n, u, ∇uなどに依存してよいが，m, dはxに対してのみの依存性をもつ．複数の従属変数u ∈ dに対する式を連立させる場合は，式 (1)でγ, fはd次元のベクトル，その他の係数はベクトルを成分にもつ行列となる．ただし，微分演算子の作用は，係数行列やベクトルuとの演算の際は，行列・ベクトルの成分内での演算に受け継がれる．たとえば，∇(u1, …, ud)T = (∇u1, …, ∇ud)T や，
といった具合である．
自然科学から工学的応用において現れる多くのPDEは，式(1)の特別な場合である．たとえば，波動方程式は，mとc以外の係数をすべてゼロにした場合であるし，非圧縮性の流体の速度，圧力場u = (v, p)T ∈ 4を記述するNavier–Stokesの式
は，
により，式(1)の形に表現できる．以下，しばらくは時間に依存しない問題を考え，空間次元に関するFEMを扱う．時間に依存する問題については 3節の最後で簡単に説明し，4.3節と4.4節で例をあげる．
重要なのは，FEMが適用可能と判断されるためには，NDSolveに与える式が各係数の依存性を含めて式(1)の形(“Coefficient Form”)として認識されることである．簡単な例として，
を考えてみよう．これは，式(1)でc = –∇u, f = –4として他の係数をすべてゼロとしたものに相当する．しかし，NDSolveに与えるPDEとして，
✕
Div[{{Derivative[1][u][x]}}.Grad[u[x], {x}], {x}] + 4 == 0
|
を入力すると，NDSolveが処理を開始する前にPDEが評価され，結果として式(1)第1項は2u´ (x)u´´ (x)とみなされてしまう．これは式(1)の Coefficient Formの形ではないため，FEMにより解くことができない(u´´ (x)がu´ (x)の係数と認識され，係数が2階導関数に依存することになってしまう)．NDSolveに式(1)の形のまま渡すには，Inactive，あるいはInactivateを用いて，
✕
Inactive[Div][{{Derivative[1][u][x]}}.Inactive[Grad][
u[x], {x}], {x}] + 4 == 0
|
のように∇の評価を保留しておけばよい．
2.2 領域の指定
任意の次元の任意の領域が指定可能である．上のポアソン方程式の例のように単純な図形であれば，DiskやPolygonなどを組み合わせて作ることもできるし，等式や不等式で表される領域であれば，ParametricRegion, ImplicitRegionなどが利用できる．さらに，写真などから作成した領域指定画像をImageMeshによりNDSolveで利用可能な領域データに変換することも可能である．
2.3 境界条件の指定
境界∂Ω上での関数値を直接与えるディリクレ境界条件は，
✕
NDSolveValue[{PDE (s) for f[x, y],
DirichletCondition[f[x, y] == bc, predicate]}, f, {x,
y} \[Element] \[CapitalOmega]]
|
のように，PDEとともに指定するだけである．ここでbcは境界での値を与える関数，predicateはf(x, y)=bcが満たされるべき境界を指定する．predicateをTrueのみにすれば，∂Ω全体が指定される．
一般化されたノイマン境界条件(ロビン(Robin)境界条件)は，NeumannValueで指定する．ロビン境界条件は境界を外向き垂直に貫く流束(flux)の成分を次の形で規定する：
である．
は∂Ω上の外向き法線(単位)ベクトル，右辺のg–quがユーザが与える値である．ただし，NeumannValueはDirichletConditionとは指定方法が異なるので注意が必要である．これは，有限要素近似において，PDEにテスト関数ϕをかけて領域Ωで積分して弱形式を得ることに由来している．式(1)第1項にϕをかけて積分すると，
·(–c∇u – αu + γ)の項は，
となる．境界∂Ω上の積分の被積分関数が，ちょうどロビン境界条件で指定されるべきものに相当している．このため，この項をg–quの積分で置き換えることで，NDSolveがこの境界条件を正しく扱うことができる．
たとえば，単位円境界のx ≤ 0の領域においてu(x,y) = 0, x ≥ 1/2においては ·∇u = xy2というディリクレ条件とノイマン条件が課せられたラプラス方程式–∇2u = 0を解くには，
![Disk[]](https://content.wolfram.com/sites/39/2020/04/km0416img7.png "Disk[]")
✕
\[CapitalOmega] = Disk[];
Subscript[\[CapitalGamma], D] =
DirichletCondition[u[x, y] == 0, x <= 0];
usol = NDSolveValue[{-Div[Grad[u[x, y], {x, y}], {x, y}] ==
NeumannValue[x*y^2, x >= 1/2], Subscript[\[CapitalGamma], D]},
u, {x, y} \[Element] \[CapitalOmega]]
|
とすればよい．(この場合のPDEはあいまいさなくCoefficient Formと認識されるので微分演算子をInactiveにする必要はないが，Inactive[Div], Inactive[Grad]としてももちろんかまわない．)
Divの前の負号は式(1)中のcを1とするためのもので，これによりノイマン条件 ·c∇u = ·∇u = xy2がそのままNeumannValueに入力できる．解をプロットすれば，次のようになる．
✕
Plot3D[usol[x, y], {x, y} \[Element] \[CapitalOmega]]
|
留意すべきは，NeumannValueが式(1)のPDEをベースにした式(4)の形で決められるため，ノイマン条件を「手で」調整する必要がある場合があることである．たとえば，ポアソン方程式–∇2u + 1/5 = 0に対して，ノイマン条件 ·∇u = xy2(x ≥ 1/2)をディリクレ条件u(x,y) = 0(x ≤ 0)とともに課すにあたり，NDSolveに与えるPDE自体を–5 ∇2u + 1 = 0としたとしよう．すると，NDSolve内では式(1)においてc = 5と認識するため，·c∇uに相当するNeumannValueは5xy2としなければならない．言われてみれば当然のことにも見えるが，ディリクレ条件と異なり，ノイマン(ロビン)条件をPDEから独立して指定するわけではないことに注意が必要である．–∇2u + 1/5 = 0の場合と–5 ∇2u + 1 = 0の場合を以下に示しておく．
入力式が–∇2u + 1/5 == 0の場合
![Disk[]](https://content.wolfram.com/sites/39/2020/04/km0416img9.png "Disk[]")
✕
\[CapitalOmega] = Disk[]; Subscript[\[CapitalGamma], D] = DirichletCondition[u[x, y] == 0, x <= 0]; |
 
✕
usol = NDSolveValue[{-Div[Grad[u[x, y], {x, y}], {x, y}] + 1/5 ==
NeumannValue[x y^2, x >= 1/2], Subscript[\[CapitalGamma], D]},
u, {x, y} \[Element] \[CapitalOmega]]
|
✕
Plot3D[usol[x, y], {x, y} \[Element] \[CapitalOmega]]
|
入力式が–5 ∇2u + 1 == 0の場合
✕
usol = NDSolveValue[{-5 Div[Grad[u[x, y], {x, y}], {x, y}] + 1 ==
NeumannValue[5 x y^2, x >= 1/2], Subscript[\[CapitalGamma], D]},
u, {x, y} \[Element] \[CapitalOmega]]
|
✕
Plot3D[usol[x, y], {x, y} \[Element] \[CapitalOmega]]
|
のようにすればよい．ロビン条件3u + ·∇u = xy2のような場合も同様で，NDSolve内のPDE左辺が–∇2u + 1/5であれば，右辺は5*NeumannValue[1/2*x*y2 – 3/2*u[x, y], x≥1/2]とする．
3. Wolfram Languageにおける非線形FEM
FEMを適用するには対象領域内にメッシュを生成する必要があるが，これについてはここでは深くは立ち入らない．興味ある方のために簡単にまとめておくと，2 次元領域についてはTriangle，3 次元領域にはTetGenとよばれるツールを利用している．TriangleはDelaunay triangulations, constrained Delaunay triangulations, conforming Delaunay triangulationsを行い，TetGenはconstrained Delaunay tetrahedralization, boundary conforming Delaunayメッシュなどの3次元領域の四面体メッシュを生成することができる．Wolfram Languageはこれらを必要なときに自動的に利用するが，ユーザが柔軟にカスタマイズして使うこともできる．詳細は，Triangleについてはここ，そしてTetGenについてはここの説明を参照されたい．
線形PDEの場合は，PDEの弱形式から離散化を経て連立1次方程式を解くが，非線形PDEを解く際にもこれを利用する．基本的な流れは
- シード(種)となる解候補の近辺で非線形PDEを線形化
- 線形化された式を離散化して解く
- シードと得られた解の差が許容する誤差以内なら終了
- 得られた解を新たなシードとして1の線形化作業へ戻る
である．つまり，非線形代数方程式のNewton–Raphson法による求解と同様のプロセスをたどる．詳細は上述のWolframドキュメント内のチュートリアルに記載があるが，簡単にまとめると以下のようになる．
まず式(1)の時間微分に関する部分を取り除くと，
であるが，
とすると，
と単純な形となる．この非線形PDEを線形化するわけだが，1 変数非線形方程式の解を数値的に求めるときのように，ある適当な関数u0をシードとして，そこから∇2·Γ (u∗) – F (u∗) = 0となる真の解u∗ヘ漸近的に近づいていく．u∗とu0の差をr = u∗ – u0とすると，
Γ, Fをu0のまわりでテイラー展開して，O(r2)の高次の項を無視する近似をすると，
となる．Γ 'やF 'などの微分は，∂Γ/∂∇u, ∂Γ/∂u, ∂F/∂∇u, ∂F/∂u などを計算して得られる．これらをu0において評価すると，式(9)は離散化した各点(節点)でのuに対する連立1次方程式となる．ここで初期・境界条件も合わせて連立させることで閉じた連立方程式となり，これによりrが得られる．
シードu0としてはデフォルトで u(x) = 0, ∀ x ∈ Ωとされるが，これはNDSolveのオプションのひとつとして，たとえばInitialSeeding→{u[x,y]==x+Exp[-Abs[y]]}などと指定することができる．線形化による漸近的求解で意図せぬ局所解に陥る可能性を考えると，問題の背景からシードをうまく与えることも有用である．また，式(13)から残差rを求める際には，左辺に現れるヤコビアン ∇·Γ '(u0) – F '(u0)の計算量が大きく，これが全体の計算時間に大きく影響する．このため，Wolfram Languageでは非線形FEMを適用する際にはNewton–Raphson法そのものではなく，Affine covariant Newton法を用いた上で，許容できる範囲内で前のステップで用いたヤコビアンを再利用するBroyden法を使い，ヤコビアンの計算回数を大幅に減らす策をとっている．
時間についての積分は，空間次元について離散化して連立方程式(行列)を得たあと，これを時間に関する常微分方程式とみなすことで種々の計算法が適用できる．線の方法(method of lines)や，場合によっては時間方向にもFEMを適用することも可能である．
4. 実行例
4.1 非線形透磁率のもとの磁場分布
電流
のまわりには磁場
が生じる．モーターのような構成で，電流をコイルに流したときどのような磁場分布が生まれるか，特に，固定子と回転子を構成する強磁性体の透磁率が磁場に依存するような非線形材料の場合で計算してみる．基礎となるモデル式は，磁場とベクトルポテンシャルの関係
と，Ampereの法則
である．これらをひとつにまとめて，さらにクーロンゲージ
をとれば，
となる．電流
はz方向のみに成分をもつとし，また問題を簡単にするためベクトルポテンシャルのx成分，y成分は定数であり，透磁率はz方向に対して一定であると仮定する．すると，式(10)で意味があるのはz成分のみとなり，スカラー量u = AzについてのPDEとなる;
透磁率μ(B)は，下の実測データよりフィッティングした式を用いた．
✕
ListPlot[BHData,PlotLabel->"Measured magnetic susceptibility"] |
✕
Clear[a1, a2, b1, b2, c1, c2];
model = a1 Exp[-(x - b1)^2/(c1^2)] + a2/(1 + Exp[(x - b2)^2/(c2^2)])^2;
fitData = FindFit[BHData, model, {a1, b1, c1, a2, b2, c2}, x];
fit = Function[x, Evaluate[model /. fitData]]
|
✕
Show[ListPlot[BHData], Plot[fit[x], {x, 0, 3}],
PlotLabel -> "Fitted curve for magnetic susceptibility"]
|
次の図はモータの断面を模式的に表したもので，黄色とオレンジ色の部分に画面に垂直方向に電流を流すとの仮定のもとで，磁場の強度分布を非線形PDE式(11)により計算してみる．
✕
mesh["Wireframe"[
"MeshElementStyle" -> (Directive[FaceForm[#],
EdgeForm[]] & /@ {Blue, Red, Gray, Orange, LightOrange,
Yellow})]]
|
透磁率の設定，電流を流す要素の指定．NDSolveでは，モーター固定子の外側では磁場がゼロになるというディリクレ境界条件を課している．
✕
B2Norm = Sqrt[
Total[Grad[
u[x, y], {x,
y}]^2] + $MachineEpsilon];(* norm of grad u *)
\[Mu]Air =
4 \[Pi]*10^-7;
\[Nu] = Piecewise[{
{-1/fit[B2Norm], ElementMarker == 2 || ElementMarker == 3}
}, -1/\[Mu]Air]*IdentityMatrix[2];
|
✕
jz = Piecewise[{
{10, ElementMarker == 4},
{-10, ElementMarker == 6}
}, 0];
|
 
✕
usol = NDSolveValue[{Inactive[Div][(\[Nu]\!\(\*
TagBox[".",
"InactiveToken",
BaseStyle->"Inactive",
SyntaxForm->"."]\)Inactive[Grad][u[x, y], {x, y}]), {x, y}] == jz,
DirichletCondition[u[x, y] == 0, x^2 + y^2 >= 0.95]},
u, {x, y} \[Element] mesh]
|
得られた結果をモータ構造のワイヤフレームとともに表示．
✕
{minsol, maxsol} = MinMax[usol["ValuesOnGrid"]];
Show[
ContourPlot[usol[x, y], {x, y} \[Element] mesh, PlotRange -> All,
ColorFunction -> "TemperatureMap",
Contours -> Range[minsol, maxsol, (maxsol - minsol)/15]],
ToBoundaryMesh[mesh]["Wireframe"],
VectorPlot[
Evaluate[{{0, 1}, {-1, 0}}.Grad[usol[x, y], {x, y}]], {x,
y} \[Element] mesh, StreamPoints -> Coarse]]
|
4.2 Driven Cavity
定常状態にある非圧縮性流体を記述するNavier–Stokes方程式
は，第1式2項めの対流項により本質的に非線形である(ここでは密度ρ = 1，外力場はゼロとした)．ここで，uはベクトルであるので，2 次元であれば，第1式は微分演算子∇が作用するのがuxかuyかの2式で構成される(下のコード参照)．2次元キャビティ内の速度場を計算してみよう．キャビティに充満した流体の上辺が右方向に一定の速度uxで駆動されること，残りの辺における流束はゼロであること，図の左下隅での圧力がゼロであることをディリクレ条件として与えた．Wolfram Languageコードは以下のとおりである．
✕
\[Nu] =.; navierstokes = {\[Rho]*{{u[x, y], v[x, y]}}.Inactive[Grad][
u[x, y], {x, y}] +
Inactive[
Div][{{-\[Nu], 0}, {0, -\[Nu]}}.Inactive[Grad][
u[x, y], {x, y}], {x, y}] + Derivative[1, 0][p][x, y],
\[Rho]*{{u[x, y], v[x, y]}}.Inactive[Grad][v[x, y], {x, y}] +
Inactive[
Div][{{-\[Nu], 0}, {0, -\[Nu]}}.Inactive[Grad][
v[x, y], {x, y}], {x, y}] + Derivative[0, 1][p][x, y],
Derivative[0, 1][v][x, y] + Derivative[1, 0][u][x, y]};
|
✕
bcs = {DirichletCondition[u[x, y] == 2, y == 1],
DirichletCondition[u[x, y] == 0, y != 1],
DirichletCondition[v[x, y] == 0, True],
DirichletCondition[p[x, y] == 0, x == 0 && y == 0]};
|
✕
op = navierstokes /. {\[Rho] -> 1, \[Nu] -> 1/1000};
{uVel, vVel, pressure} =
NDSolveValue[{op == {0, 0, 0}, bcs}, {u, v, p}, {x, y} \[Element]
Rectangle[{0, 0}, {1, 1}],
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1},
"MeshOptions" -> {"MaxCellMeasure" -> 0.0001}}];
|
得られる速度場を可視化してみよう．
✕
Show[
StreamPlot[{uVel[x, y], vVel[x, y]}, {x, 0, 1}, {y, 0, 1},
Axes -> None, Frame -> None,
StreamPoints -> {Automatic, Scaled[0.02]}],
ToBoundaryMesh[uVel["ElementMesh"]]["Wireframe"]]
|
圧力分布は以下のようである．
✕
Plot3D[pressure[x, y], {x, 0, 1}, {y, 0, 1}, PlotRange -> {-0.5, 1.5},
PlotPoints -> 80, Boxed -> True]
|
4.3 Gray–Scottモデル
反応拡散系とよばれる，化学反応と物質の拡散による複数の物質の濃度変化をモデル化した非線形連立PDE(Gray–Scottモデル)を計算してみるのが次の例である．外部から原料となる化学物質Uが，別の物質Vで満たされた反応容器内に連続的に導入され，自己触媒反応
を経て，Uが最終生成物Pに変化し，Pは系外へ排出されるとする．UとVの濃度u, vの時間変化は
によって記述されるとするのがこのモデルである．Du, Dvはそれぞれの拡散係数，Fは物質Uの補充率，Kは反応V→Pの速さを規定するパラメータである．以下に，式の入力から可視化(アニメ化)までのコードを示しておく．
✕
eqn = {
D[u[t, x, y], t] +
Inactive[Plus][
Inactive[
Div][{{-c1, 0}, {0, -c1}}.Inactive[Grad][
u[t, x, y], {x, y}], {x, y}], (v[t, x, y]^2 + f)*
u[t, x, y]] == f,
D[v[t, x, y], t] +
Inactive[Plus][
Inactive[
Div][{{-c2, 0}, {0, -c2}}.Inactive[Grad][
v[t, x, y], {x, y}], {x,
y}], (-u[t, x, y]*v[t, x, y] + f + k)*v[t, x, y]] == 0
} //. {c1 -> 2.*10^-5, c2 -> c1/4, f -> 0.04, k -> 0.06};
|
✕
ics = {u[0, x, y] == 1/2,
v[0, x, y] == If[x^2 + y^2 <= 0.025, 1., 0.]};
bcs = {DirichletCondition[u[t, x, y] == 0, True],
DirichletCondition[v[t, x, y] == 0, True]};
|
✕
{ufun, vfun} =
NDSolveValue[{eqn, bcs, ics}, {u, v}, {x, y} \[Element] Disk[], {t,
0, 3000},
Method -> {"TimeIntegration" -> {"IDA",
"MaxDifferenceOrder" -> 2},
"PDEDiscretization" -> {"MethodOfLines",
"DifferentiateBoundaryConditions" -> True,
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.002}}}}];
|
✕
{vmin, vmax} = MinMax[vfun["ValuesOnGrid"]];
frames = Table[
Rasterize[
ContourPlot[vfun[t, x, y], {x, -1, 1}, {y, -1, 1},
Contours -> Range[vmin, vmax, (vmax - vmin)/4],
PlotRange -> All], RasterSize -> Large], {t, 100, 2000, 20}];
|
✕
ListAnimate[frames] |
4.4 円柱のまわりの流れ
時間に依存する非圧縮性流体の流れを記述するのは，式(12)に時間変化の項を付け加えた
である．二枚の無限に広い平行平板の間の空間を流体が流れる状況で，この空間に無限に長い円柱を流れと垂直方向に置いたときの流速の分布を計算してみる．平板と円柱に垂直な面をxy平面として，未知数は速度(u,v)と圧力pとなる．Wolfram Languageコードは以下のとおりである．
Navier–Stokes方程式
✕
transientnavierstokes = {\[Rho]*{{u[t, x, y], v[t, x, y]}}.Inactive[
Grad][u[t, x, y], {x, y}] +
Inactive[
Div][{{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
u[t, x, y], {x, y}], {x, y}] +
Derivative[0, 1, 0][p][t, x, y] + \[Rho]*
Derivative[1, 0, 0][u][t, x,
y], \[Rho]*{{u[t, x, y], v[t, x, y]}}.Inactive[Grad][
v[t, x, y], {x, y}] +
Inactive[
Div][{{-\[Mu], 0}, {0, -\[Mu]}}.Inactive[Grad][
v[t, x, y], {x, y}], {x, y}] +
Derivative[0, 0, 1][p][t, x, y] + \[Rho]*
Derivative[1, 0, 0][v][t, x, y],
Derivative[0, 0, 1][v][t, x, y] + Derivative[0, 1, 0][u][t, x, y]};
|
流域のサイズの設定と，流入口での速度分布の設定．ある時刻で不連続に流速がゼロから非ゼロに変化しないように，滑らかな流速変化を与えるrampFunctionを定義する．流域のサイズと流速等から，ここでの流れのレイノルズ数は約200である．
✕
rules = {length -> 2.2, height -> 0.41};
\[CapitalOmega] =
RegionDifference[Rectangle[{0, 0}, {length, height}],
Disk[{1/5, 1/5}, 1/20]] /. rules;
rmf = RegionMember[\[CapitalOmega]];
rampFunction[min_, max_, c_, r_] :=
Function[t, (min*Exp[c*r] + max*Exp[r*t])/(Exp[c*r] + Exp[r*t])]
ramp = rampFunction[0, 1, 4, 5];
GraphicsRow[{
Show[BoundaryDiscretizeRegion[\[CapitalOmega]],
VectorPlot[{4*1.5*y*(0.41 - y)/0.41^2, 0}, {x, 0, 2.2}, {y, 0,
0.41}, VectorScale -> Small, VectorStyle -> Red,
VectorMarkers -> Placed["Arrow", "Start"],
VectorPoints -> Table[{0, y}, {y, 0.05, 0.35, 0.075}],
ImageSize -> Large]]
, Plot[ramp[t], {t, -1, 10}, PlotRange -> All, ImageSize -> Large,
AspectRatio -> 1/5]
}]
|
境界条件と初期条件
✕
op = transientnavierstokes /. {\[Mu] -> 10^-3, \[Rho] -> 1};
bcs = {
DirichletCondition[
u[t, x, y] == ramp[t]*4*1.5*y*(height - y)/height^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};
|
t = 0 から10までの速度分布の変化を，tをモニターしながらNDSolveにより計算．一般的なPC(3.1GHz Intel Core i5，メモリ16GB)で要6分程度．
  
✕
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.0002}}}},
EvaluationMonitor :> (currentTime = t;)];]
|
各点での速度の絶対値をもとに色付けし，アニメーションを作成．
✕
{minX, maxX} =
MinMax[Sqrt[xVel["ValuesOnGrid"]^2 + yVel["ValuesOnGrid"]^2]];
mesh = xVel["ElementMesh"];
frames = Table[
Rasterize[
ContourPlot[
Norm[{xVel[t, x, y], yVel[t, x, y]}], {x, 0, 2.2}, {y, 0, 0.41},
PlotRange -> All, AspectRatio -> Automatic,
ColorFunction -> "TemperatureMap",
Contours -> Range[minX, maxX, (maxX - minX)/7], Axes -> False,
Frame -> None,
RegionFunction -> Function[{x, y, z}, rmf[{x, y}]]],
RasterSize -> 4*{360, 68}, ImageSize -> 2*{360, 68}], {t, 4, 10,
1/20}];
|
✕
ListAnimate[frames] |
領域に対して生成されたメッシュは以下のように確認できる．
✕
ToElementMesh[\[CapitalOmega], "MaxCellMeasure" -> 0.0002]["Wireframe"] |
5. まとめ
ここまで見てきたように，Mathematica 12(Wolfram Language 12)では有限要素法の適用範囲が大幅に広がり，Navier–Stokes方程式を始めとする多くの非線形偏微分方程式の求解が可能となった．シンボリックな計算に強みを持つWolfram Languageにより，個々のPDEの形によらず，高度な一般性を保ったまま統一的な処理を効率的に実行することが可能となっている．FEM関連の内部処理の詳細が公開されていることは前述したが，同時に弾性解析，音響解析，熱・振動伝導解析など多くの分野における応用例も，Mathematicaチュートリアルで詳しく解説されているので，参考にしていただければ幸いである．
| Get full access to the latest Wolfram Language functionality with a Mathematica 12.1 or Wolfram|One trial. |
Engage with the code in this post by downloading the Wolfram Notebook
Wonderfull ! At last, Mathematica becomes numerical solver of Nonlinear PDE.