Wolfram Blog
Koji Maruyama

非線形偏微分方程式への有限要素法の適用

2020年4月29日 -Koji Maruyama, Sales Engineer

非線形偏微分方程式への有限要素法の適用

数学12has 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 theWolfram Languagefor 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 “有限元程序”,基于本文的大部分。


1. はじめに

manbet万博app沃尔夫勒姆研究社の旗艦製品であるMathematicaは,5,000 を超える組み込み関数を有するWolfram Languageを駆動する.数理モデリング,解析の基本となる常・偏微分方程式の分野においては,これらをシンボリックに,あるいは数値的に解くための強力なソルバを搭載している.最近は有限要素法(FEM) を利用した数値的求解機能が大幅に強化され,偏微分方程式(PDE)を任意の領域上で解いたり,固有値・固有関数を求めたりすることが可能となった.ここでは,最新のバージョン12における非線形偏微分方程式のFEMによる求解を中心に,現実的な問題に応用する上での流れを例とともに紹介する.なお,有限要素法を用いて非線形PDEを解くワークフローの詳細,コードはすべて公開されている.MathematicaのWolframドキュメント内で,チュートリアル“FiniteElementProgramming“を参照いただきたい。

2. 微分方程式の数値求解ワークフロー

Wolfram Languageにおいて微分方程式を数値的に解く際の関数(コマンド)は,NDSolve,あるいはNDSolveValueのふたつある。これら二つは出力のフォーマットが若干异なるだけで,中の处理はまったく同じものなので,以下,本文中では表记が短い"NDSolve”, コード例中では出力の扱いが簡便な"NDSolveValue"で表記する.Mathematica上でFEMを利用するには,パッケージをロードする :

Needs[

Needs["NDSolve`FEM`"]

あとは,NDSolveにPDE,領域,初期・境界条件を与えるだけである.たとえば,単位円上のポアソン方程式2u= 1を例にとると,境界条件としてx0にある境界でu= 0とすると,

eqn = -Inactive

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

Plot3D[usol[x, y], {x, y} \[Element] Subscript[\[CapitalOmega], D]]

とプロットできる.

2.1 式の入力

NDSolveの中の有限要素法が現在適用可能な偏微分方程式は次の形をもたなければならない :

ここで,解くべき従属変数un上の1次元关数のときは,m,d,a,fはスカラー,α,γ,βn次元ベクトル,cn*n行列であるま。た,c,a,f,α,γ,βt,xn,u,uなどに依存してよいが,m,dxに対してのみの依存性をもつ.複数の従属変数udに対する式を連立させる場合は,式 (1)でγ,fd次元のベクトル,その他の系数はベクトルを成分にもつ行列となる。ただし,微分演算子の作用は,系数行列やベクトルuとの演算の際は,行列・ベクトルの成分内での演算に受け継がれる.たとえば,(u1, …,ud)T= (u1, …,ud)Tや,

といった具合である。

自然科学から工学的応用において现れる多くのPDEは,式(1)の特别な场合である。たとえば,波动方程式は,mc以外の系数をすべてゼロにした场合であるし,非圧缩性の流体の速度,圧力场u= (v,p)T4を記述するNavier–Stokesの式

は,

により,式(1)の形に表現できる.以下,しばらくは時間に依存しない問題を考え,空間次元に関するFEMを扱う.時間に依存する問題については 3節の最後で簡単に説明し,4.3節と4.4節で例をあげる.

重要なのは,FEMが適用可能と判断されるためには,NDSolveに与える式が各係数の依存性を含めて式(1)の形(“Coefficient Form”)として認識されることである.簡単な例として,

を考えてみよう.これは,式(1)でc= –u,f= –4として他の系数をすべてゼロとしたものに相当する。しかし,NDSolveに与えるPDEとして,

Div

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)の形のまま渡すには,待用,あるいは灭活を用いて,

待用

待用[Div][{{Derivative[1][u][x]}}.Inactive[Grad][ u[x], {x}], {x}] + 4 == 0

のように∇の評価を保留しておけばよい.

2.2 領域の指定

任意の次元の任意の領域が指定可能である.上のポアソン方程式の例のように単純な図形であれば,Disk多边形などを組み合わせて作ることもできるし,等式や不等式で表される領域であれば,ParametricRegion,ImplicitRegionなどが利用できる。さらに,写真などから作成した领域指定画像をImageMeshによりNDSolveで利用可能な領域データに変換することも可能である.

2.3境界条件の指定

境界∂Ω上での关数値を直接与えるディリクレ境界条件は,

NDSolveValue

NDSolveValue[{PDE (s) for f[x, y], DirichletCondition[f[x, y] == bc, predicate]}, f, {x, y} \[Element] \[CapitalOmega]]

のように,PDEとともに指定するだけである.ここでbcは境界での値を与える関数,谓语f(x,y)=bcが満たされるべき境界を指定する。谓语Trueのみにすれば,∂Ω全体が指定される.

一般化されたノイマン境界条件(ロビン(Robin)境界条件)は,NeumannValueで指定する.ロビン境界条件は境界を外向き垂直に貫く流束(flux)の成分を次の形で規定する:

である.Ω上の外向き法線(単位)ベクトル,右辺のgがユーザが与える値である。ただし,NeumannValueDirichletConditionとは指定方法が异なるので注意が必要である。これは,有限要素近似において,PDEにテスト关数ϕをかけて領域Ωで積分して弱形式を得ることに由来している。式(1)第1項にϕをかけて積分すると,·(–cuαu+γ)の項は,

となる.境界∂Ω上の積分の被積分関数が,ちょうどロビン境界条件で指定されるべきものに相当している.このため,この項をgの积分で置き换えることで,NDSolveがこの境界条件を正しく扱うことができる.

たとえば,単位円境界のx0の領域においてu(x,y) = 0,x1/2においては·u=xy2というディリクレ条件とノイマン条件が課せられたラプラス方程式2u= 0を解くには,

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と認識されるので微分演算子を待用にする必要はないが,待用[Div],待用[Grad]としてももちろんかまわない.)

Divの前の负号は式(1)中のcを1とするためのもので,これによりノイマン条件·cu=·u=xy2がそのままNeumannValueに入力できる.解をプロットすれば,次のようになる.

Plot3D

Plot3D[usol[x, y], {x, y} \[Element] \[CapitalOmega]]

留意すべきは,NeumannValueが式(1)のPDEをベースにした式(4)の形で決められるため,ノイマン条件を「手で」調整する必要がある場合があることである.たとえば,ポアソン方程式2u+ 1/5 = 0に対して,ノイマン条件·u=xy2(x1/2)をディリクレ条件u(x,y) = 0(x0)とともに課すにあたり,NDSolveに与えるPDE自体を-52u+ 1 = 0としたとしよう.すると,NDSolve内では式(1)においてc= 5と認識するため,·cuに相当するNeumannValue5xy2としなければならない。言われてみれば当然のことにも见えるが,ディリクレ条件と异なり,ノイマン(ロビン)条件をPDEから独立して指定するわけではないことに注意が必要である。2u+ 1/5 = 0の场合と-52u+ 1 = 0の場合を以下に示しておく.

入力式が2u+ 1/5== 0の場合

Disk[]

\[CapitalOmega] = Disk[]; Subscript[\[CapitalGamma], D] = DirichletCondition[u[x, y] == 0, x <= 0];

USOL = NDSolveValue
USOL = NDSolveValue

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

Plot3D[usol[x, y], {x, y} \[Element] \[CapitalOmega]]

入力式が-52u+ 1 == 0の場合

USOL = NDSolveValue

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

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], x1/2]とする。

3. Wolfram Languageにおける非線形FEM

FEMを適用するには対象領域内にメッシュを生成する必要があるが,これについてはここでは深くは立ち入らない.興味ある方のために簡単にまとめておくと,2 次元領域についてはTriangle,3 次元領域にはTetGenとよばれるツールを利用している.TriangleはDelaunay三角,约束Delaunay三角剖分,贴合Delaunay三角を行い,TetGenは约束Delaunay四面体,边界贴合德洛奈メッシュなどの3次元领域の四面体メッシュを生成することができる.manbet万博appWolfram语言はこれらを必要なときに自动的に利用するが,ユーザが柔软にカスタマイズして使うこともできる。详细は,三角についてはここ,そしてTetGenについてはここの説明を参照されたい.

線形PDEの場合は,PDEの弱形式から離散化を経て連立1次方程式を解くが,非線形PDEを解く際にもこれを利用する.基本的な流れは

  1. シード(種)となる解候補の近辺で非線形PDEを線形化
  2. 線形化された式を離散化して解く
  3. シードと得られた解の差が許容する誤差以内なら終了
  4. 得られた解を新たなシードとして1の线形化作业へ戻る

である.つまり,非線形代数方程式のNewton–Raphson法による求解と同様のプロセスをたどる.詳細は上述のWolframドキュメント内のチュートリアルに記載があるが,簡単にまとめると以下のようになる.

まず式(1)の時間微分に関する部分を取り除くと,

であるが,

とすると,

と単純な形となる.この非線形PDEを線形化するわけだが,1 変数非線形方程式の解を数値的に求めるときのように,ある適当な関数u0をシードとして,そこから2·Γ(u*) -F(u*) = 0となる真の解u*ヘ漸近的に近づいていく.u*u0の差をr=u*u0とすると,

Γ,Fu0のまわりでテイラー展開して,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法を使い,ヤコビアンの計算回数を大幅に減らす策をとっている.

時間についての積分は,空間次元について離散化して連立方程式(行列)を得たあと,これを時間に関する常微分方程式とみなすことで種々の計算法が適用できる.線の方法(线的方法)や,场合によっては时间方向にもFEMを适用することも可能である。

4. 実行例

4.1 非線形透磁率のもとの磁場分布

電流のまわりには磁場が生じる。モーターのような构成で,电流をコイルに流したときどのような磁场分布が生まれるか,特に,固定子と回転子を构成する强磁性体の透磁率が磁场に依存するような非线形材料の场合で计算してみる。基础となるモデル式は,磁场とベクトルポテンシャルの关系と,安培の法则である。これらをひとつにまとめて,さらにクーロンゲージをとれば,

となる.電流z方向のみに成分をもつとし,また問題を簡単にするためベクトルポテンシャルのx成分,y成分は定数であり,透磁率はz方向に対して一定であると仮定する.すると,式(10)で意味があるのはz成分のみとなり,スカラー量u=AzについてのPDEとなる;

透磁率μ(B)は,下の実测データよりフィッティングした式を用いた。

ListPlot

ListPlot[BHData,PlotLabel->"Measured magnetic susceptibility"]

明确

明确[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

Show[ListPlot[BHData], Plot[fit[x], {x, 0, 3}], PlotLabel -> "Fitted curve for magnetic susceptibility"]

次の図はモータの断面を模式的に表したもので,黄色とオレンジ色の部分に画面に垂直方向に電流を流すとの仮定のもとで,磁場の強度分布を非線形PDE式(11)により計算してみる.

啮合

啮合["Wireframe"[ "MeshElementStyle" -> (Directive[FaceForm[#], EdgeForm[]] & /@ {Blue, Red, Gray, Orange, LightOrange, Yellow})]]

透磁率の設定,電流を流す要素の指定.NDSolveでは,モーター固定子の外侧では磁场がゼロになるというディリクレ境界条件を课している。

B2Norm = SQRT

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

JZ =分段[{{10,ElementMarker == 4},{-10,ElementMarker == 6}},0];

USOL = NDSolveValue
USOL = NDSolveValue

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

{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式は微分演算子が作用するのがuxuyかの2式で构成される(下のコード参照)0.2次元キャビティ内の速度场を计算してみよう。キャビティに充満した流体の上辺が右方向に一定の速度uxで駆動されること,残りの辺における流束はゼロであること,図の左下隅での圧力がゼロであることをディリクレ条件として与えた.Wolfram Languageコードは以下のとおりである.

navierstokes

\[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

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

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},轴 - >无,框架 - >无,StreamPoints  - > {自动的,缩放[0.02]}],ToBoundaryMesh [uVel [ “ElementMesh”]] [ “线框”]]

圧力分布は以下のようである.

Plot3D

Plot3D[pressure[x, y], {x, 0, 1}, {y, 0, 1}, PlotRange -> {-0.5, 1.5}, PlotPoints -> 80, Boxed -> True]

4.3灰色斯科特モデル

反応拡散系とよばれる,化学反応と物質の拡散による複数の物質の濃度変化をモデル化した非線形連立PDE(Gray–Scottモデル)を計算してみるのが次の例である.外部から原料となる化学物質Uが,別の物質Vで満たされた反応容器内に連続的に導入され,自己触媒反応

を経て,Uが最終生成物Pに変化し,Pは系外へ排出されるとする.UとVの濃度u,vの時間変化は

によって記述されるとするのがこのモデルである.Du,Dvはそれぞれの拡散係数,Fは物質Uの補充率,Kは反応VPの速さを規定するパラメータである.以下に,式の入力から可視化(アニメ化)までのコードを示しておく.

eqn = {D

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

集成电路=【U

集成电路=【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

{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

{VMIN,VMAX =} MINMAX [vfun [ “ValuesOnGrid”]];帧=表[栅格化[ContourPlot [vfun [T,X,Y],{X,-1,1},{Y,-1,1},轮廓 - >范围[VMIN,VMAX,(VMAX  -  VMIN)/4],PlotRange  - >全部],RasterSize  - >大],{T,100,2000,20}];

ListAnimate

ListAnimate[frames]

4.4 円柱のまわりの流れ

時間に依存する非圧縮性流体の流れを記述するのは,式(12)に時間変化の項を付け加えた

である.二枚の無限に広い平行平板の間の空間を流体が流れる状況で,この空間に無限に長い円柱を流れと垂直方向に置いたときの流速の分布を計算してみる.平板と円柱に垂直な面をxy平面として,未知数は速度(u,v)と圧力pとなる.Wolfram Languageコードは以下のとおりである.

Navier–Stokes方程式

transientnavierstokes

transientnavierstokes = {\ [的Rho] * {{U [T,X,Y],V [T,X,Y]}}。无效[梯度] [U [T,X,Y],{X,Y}]+无效[股利] [{{ -  \μ,0},{0, -  \ [穆]}}。无效[梯度] [U [T,X,Y],{X,Y}],{的x,y}] +导数[0,1,0] [p] [T,X,Y] + \ [的Rho] *衍生物[1,0,0] [U] [T,X,Y],\[的Rho] * {{U [T,X,Y],v [T,X,Y]}}。无效[梯度] [v [T,X,Y],{X,Y}] +无效[股利] [{{ -  \μ,0},{0, -  \ [穆]}}。无效[梯度] [v [T,X,Y],{X,Y}],{X,Y}] +衍生物[0,0,1] [p] [T,X,Y] + \ [的Rho] *衍生物[1,0,0] [v] [T,X,Y],衍生物[0,0,1] [v] [T,X,Y] +导数[0,1,0] [U] [T,X,Y]};

流域のサイズの设定と,流入口での速度分布の设定。ある时刻で不连続に流速がゼロから非ゼロに変化しないように,滑らかな流速変化を与えるrampFunctionを定義する.流域のサイズと流速等から,ここでの流れのレイノルズ数は約200である.

规则

规则= {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

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分程度.

动态
动态
动态

动态["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

{其minX,maxX的} = MINMAX [的Sqrt [xVel [ “ValuesOnGrid”] ^ 2 + yVel [ “ValuesOnGrid”] ^ 2]];目= xVel [ “ElementMesh”];帧=表[栅格化[ContourPlot [范数[{xVel [T,X,Y],yVel [T,X,Y]}],{X,0,2.2},{Y,0,0.41},PlotRange  - >所有,ASPECTRATIO  - >自动,ColorFunction的 - > “TemperatureMap”,轮廓 - >范围[其minX,maxX的,(maxX的 - 其minX)/ 11],轴 - >假,帧 - >无,RegionFunction  - >函数[{X,Y,Z},RMF [{X,Y}]]],RasterSize  - > 4 * {360} 68,IMAGESIZE  - > 2 * {360,68}],{T,4,10,1/20}]。

ListAnimate

ListAnimate[frames]

領域に対して生成されたメッシュは以下のように確認できる.

ToElementMesh

MaxCellM ToElementMesh [\ [CapitalOmega]。easure" -> 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数学12。1要么Wolfram|Onetrial.

Leave a Comment

One Comment


Toshiaki Itoh

Wonderfull ! At last, Mathematica becomes numerical solver of Nonlinear PDE.

Posted by Toshiaki Itoh May 20, 2020 at 10:01 pm


Leave a comment