有限要素法で偏微分方程式を解く

はじめに	古典的な偏微分方程式
領域	偏微分方程式系

はじめに

このチュートリアルは，有限要素法(FEM)の概要を，NDSolveに実装されているように，紹介するものである．ノートブックは，偏微分方程式を解くための有限要素法の概念を紹介する．まず，典型的なワークフローについて説明する．領域，境界条件，方程式の設定の後には，NDSolveを使った偏微分方程式の解法について述べる．その後，解の可視化とアニメーション化を紹介し，有限要素法の理論的な部分について提示する．ポアソン(Poisson)方程式や熱方程式といった標準的な偏微分方程式についても説明する．また，結合偏微分方程式についても構造力学および流体動力学の例とともに紹介する．

なぜ有限要素なのか

偏微分方程式で明示的な閉形式の解が得られることはほとんどない．有限要素法は，偏微分方程式を数値的に解く方法である．これは，最低でも2つの理由から重要である．まず，有限要素法は，偏微分方程式をほとんどどのような形の領域においてでも解くことができる．さらに，この方法は，いろいろな種類の偏微分方程式に使うのに適している．特定の領域上の特定の偏微分方程式を解くのに，よりよい方法を見付けることは，ほとんどの場合可能であるが，有限要素法は，広範囲の偏微分方程式に対して使用することができる．つまり，有限要素法は，以下を取り扱えるという点で重要である：

任意形状の領域

広範囲の偏微分方程式

有限要素の分析に必要なもの

偏微分方程式を有限要素法を使って解くには，次の3つが必要である：

領域の離散的な表現，つまりメッシュ

偏微分方程式

方程式を領域に関連付ける境界条件

このセクションでは，偏微分方程式とその境界条件について見ていく．有限要素メッシュは，ToElementMeshで生成することができる．有限要素メッシュの生成についての詳しい説明は，「要素メッシュの生成」をご覧いただきたい．

NDSolveに実装される有限要素法のスコープ

有限要素法の実装の現行バージョンでは，以下の機能をサポートする：

一次，二次，三次の空間次元における定常（時間独立）の偏微分方程式

一次，二次，三次の空間次元における非定常（時間従属）の偏微分方程式

二次までの空間，あるいは任意次数の時間での導関数を持つ単一および結合の偏微分方程式

変数係数を持つ線形および非線形の偏微分方程式

線形変数係数のディリクレ(Dirichlet)条件，一般化されたノイマン(Neumann (Robin))境界値（非線形も可能），これらの境界条件を組み合わせたものが可能である

曲線の境界を持つ自動の静的メッシュ生成

手作業での精緻化オプションを含む，一次および二次のメッシュ

解の過程を阻止し，いつでも低レベルデータにアクセスできる

以下のセクションでは，この機能の概要を提供する．特定の物理や工学の分野において，どのような偏微分方程式が解けるかについての詳しいリストは，「PDE Modelsの概要」ページおよびガイドページ「偏微分方程式項」も参照されたい．

有限要素の機能を使いやすくするためには，有限要素法パッケージをロードする．

有限要素法のパッケージをロードする．

Wolfram Language code: Needs["NDSolve`FEM`"]

領域

有限要素の分析を行うためには，偏微分方程式を解く領域を定義しなければならない．空間独立変数のそれぞれについて，{v,v_min,v_max}を使って矩形領域を指定することができる．厳密な引数の指定については，NDSolveValueに詳細がある．

ポアソン方程式を解き，矩形領域での解を可視化する．

Wolfram Language code:

ufun = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == 10, 
	u[x, 0] == 0, u[x, 1] == 1}, 
	u, {x, 0, 1}, {y, 0, 1}];
Plot3D[ufun[x, y], {x, 0, 1}, {y, 0, 1}, ColorFunction -> "TemperatureMap"]

表記 vars∈Ωを使って任意の形の領域を指定することができる．ここでΩは，RegionQ[Ω]がTrueを返す領域である．

矩形領域ではない領域を指定する1つの方法に，ブール述部を使う方法がある．述部は，TrueあるいはFalseを返す関数である．「述部」という用語には数学的にいろいろな意味があるので，「ブール述部」という用語を使い，この述部がブールを返すことを強調する．領域を表すために，ブール述部を例えばRegionPlotおよびRegionPlot3Dのような関数に使うことができる．

領域の設定と可視化．

Wolfram Language code:

Ω = ImplicitRegion[!((x - 5) ^ 2 + (y - 5) ^ 2 ≤ 3 ^ 2), {{x, 0, 5}, {y, 0, 10}}];
RegionPlot[Ω, AspectRatio -> Automatic]

このチュートリアルでは，領域やメッシュの作成については強調せず，有限要素法を使って偏微分方程式を解くことに焦点を当てる．このことについての詳細は，「要素メッシュの生成」をご覧いただきたい．

古典的な偏微分方程式

偏微分方程式の係数形式

NDSolveに実装されるように，有限要素法を使ってどのような方程式を解くことが可能であろうか．での単一の偏微分方程式を考える．

偏微分方程式はで定義される．ここでは，解が求められる従属変数である．係数，，，はスカラー，，，はベクトル，は × 行列である．

次に挙げるのは，よく知られている偏微分方程式のいくつかと，その対応する係数である．(1)の一般性を示すために，特定の方程式に関係のある部分は赤で，関係のない部分は灰色で示す．

ラプラス(Laplace)方程式には，拡散項が含まれるだけである．

ポアソン(Poisson)の方程式をモデル化するには，少し変更し，負荷項を加えるだけでよい．

ヘルムホルツ(Helmholtz)の方程式には，反応項を加える．

対流・拡散・反応のタイプの方程式もよく使われる偏微分方程式である．上の例と比べると，これらには，追加の対流項が含まれる．

ここまで見てきた偏微分方程式は，定常である．つまり，時間に依存しない．熱方程式は，ポアソン方程式に時間依存を加えたものであり，以下の形式を持つ．

同様に，波動方程式は以下で与えられる．

方程式(1)は，次数2までの空間導関数を提供するので，広範囲に渡るさまざまな現象をモデル化するためのコンポーネントを提供する．

ディリクレ(Dirichlet)条件のポアソン(Poisson)方程式

ポアソン方程式は，ディリクレ境界条件を使って，領域上で解かれる．まず，方程式を解く領域を定義しなければならない．それから，方程式と境界条件が定義される．最後に方程式がその領域上で解かれる．

領域を指定する方法に，ブール述部を使う方法がある．

領域の設定と可視化．

Wolfram Language code:

Ω = ImplicitRegion[!((x - 5) ^ 2 + (y - 5) ^ 2 ≤ 3 ^ 2), {{x, 0, 5}, {y, 0, 10}}];
RegionPlot[Ω, AspectRatio -> Automatic]

次に，ポアソン演算子を定義する．

偏微分方程式の演算子の設定．

Wolfram Language code: op = -Laplacian[u[x, y], {x, y}] - 20;

第3のステップとして，境界条件を設定する．最初のディリクレ境界条件は，x0∧8≤y≤10の評価が領域の境界上でTrueを返す場合には必ず，従属変数を0の値に強制する．この場合，これは境界の左上部分である．2つ目のディリクレ境界条件は，が領域の境界上でTrueを返す場合には必ず，従属変数の値を指定する．境界条件については，後のセクションで詳しく説明する．

ディリクレ境界条件の設定．

Wolfram Language code:

Subscript[Γ, D] = {DirichletCondition[u[x, y] == 0, x == 0 && 8 ≤ y ≤ 10], 
	DirichletCondition[u[x, y] == 100, (x - 5) ^ 2 + (y - 5) ^ 2 == 3 ^ 2]};

解のステップでは，NDSolveは偏微分方程式，境界条件，領域を必要とする．

偏微分方程式を解く．

Wolfram Language code: ufun = NDSolveValue[{op == 0, Subscript[Γ, D]}, u, {x, y}∈Ω]

内部では，領域は有限要素メッシュに変換され，有限要素法が自動的に方程式を解くのに使われる．

解を可視化するには，等高線プロットが使える．

解を可視化する．

Wolfram Language code: ContourPlot[ufun[x, y], {x, y}∈Ω, ColorFunction -> "Temperature", AspectRatio -> Automatic]

偏微分方程式と境界条件

NDSolveおよびその関連の関数では，ディリクレ条件，ノイマン値，周期境界条件という3種類の空間境界条件を指定することができる．

ディリクレ境界条件は，境界の一部における値の従属変数上で条件を指示する．

一般化されたノイマン境界値（ロビン(Robin)値としても知られる）は，値を指定する．値は，境界の部分上の外向きの法線での流束を指示する．

有限要素法は，微分方程式の弱形式に基づく．この形式は，方程式(1)を取り，これにいわゆるテスト関数を掛け，領域上で積分することによって得られる．

部分ごとの積分で以下が得られる．

この過程は，内部的に行われる．(11)内の境界積分の被積分関数は， (9) で置換され，以下を返す．

は，弱形式の境界積分の被積分関数(11)の値を指定するので，NeumannValueという名前を持つ．

部分ごとに積分を行って方程式(11)を得る代りに，発散定理とグリーンの恒等式を使うこともできる．

周期境界条件は，従属変数を，境界の2つの独立した部分間に与えらる関係に従って動作させる．ほかの境界条件もあり得るが，現在のところは実装されていない．

ほとんどの場合，ディリクレ境界条件は，特定の方程式に関連付けられている必要はなく，方程式とは独立して指定することができる．ディリクレ境界条件は方程式として指定される．一般化されたノイマン値は，それに対して，値を与えることで指定される．満たされる方程式が値では陰的であるからである．ノイマン値は，数学的に偏微分方程式に結び付けられている．

実用的な理由から，NDSolveおよび関連関数NeumannValueは，方程式の一部として与えられなければならない．結合偏微分方程式系を考えてみよう．その偏微分方程式系のどの指定の単一偏微分方程式に対してでも任意の与えられたノイマン値を明瞭的に指定できると望ましい．偏微分方程式の値が関連付けられているノイマン値から，これを（明瞭的に）導き出すことは不可能である．偏微分方程式のNeumannValue部分を作ることで，明瞭にこの問題を解決できる．偏微分方程式内でDirichletConditionを与えることもできる．

{op[u[x₁,…]] Γ_N,Γ_D,Γ_P}	微分演算子 op，ノイマン境界値Γ_N，ディリクレ境界条件方程式 Γ_Dおよび周期境界条件Γ_Pを持つ偏微分方程式
{op₁[u[x₁,…],v[x₁,…],…] Γ_N₁,op₂[u[x₁,…],v[x₁,…],…] Γ_N₂,…,Γ_D,Γ_P}	結合偏微分方程式の演算子 op_j，ノイマン境界値Γ_{N_j}，ディリクレと周期の境界条件の方程式Γ_DとΓ_Pの系

境界条件で偏微分方程式を指定する．

DirichletCondition，NeumannValue，PeriodicBoundaryConditionはどれも，第2引数として，条件または値が適用される境界上の位置を示す述部が必要である．さらに，PeriodicBoundaryConditionは，境界の2つの部分間の関係を指定する第3引数を持つ．

偏微分方程式は通常，領域のある部分について，一意的に解ける少なくとも1つのディリクレ境界条件が必要である．このため，ディリクレ境界条件は，必須の境界条件とも呼ばれる．

一般化されたノイマン値は，弱形式の積分形式から使用されることになるので，通常自然条件とも呼ばれる．

純粋なノイマン境界方程式は，である一般化されたノイマン境界方程式(9)である．

も指定されない場合には，ノイマンのゼロ境界値が自動的に示唆される()．

そうすると，ノイマン値はゼロであり，項が(12)では削除される．これは，境界に流束や応力がない場合に通常起る．領域境界の一部に対して境界条件が明示的に指定されていない場合には，ノイマンのゼロ値が想定される．

ノイマン(Neumann)値のポアソン(Poisson)方程式

一般化されたノイマン境界値で，ポアソン方程式を解くことを考えてみよう．

同じ領域と偏微分方程式の演算子を使う．

Wolfram Language code:

Ω = ImplicitRegion[!((x - 5) ^ 2 + (y - 5) ^ 2 ≤ 3 ^ 2), {{x, 0, 5}, {y, 0, 10}}];
op = -Laplacian[u[x, y], {x, y}] - 20;

ポアソン方程式は以下で与えられる．

対応するノイマン境界方程式は以下の通りである．

ノイマンタイプの方程式の正規導関数がどのように偏微分方程式に関連付けられているかに注意されたい．ノイマン境界方程式では，右辺のとの値を指定する．左辺は，偏微分方程式を通して暗黙的に指定されている．偏微分方程式の係数が恒等行列として与えられている場合には，ノイマン条件の左辺の係数は，その同じ係数である．

ノイマン境界値は，値とを通して指定される．通常，である場合にはノイマン境界値を使い，である場合には，一般化されたノイマン境界値，つまりロビン境界値を使う．一般化されたノイマンとノイマン境界値の両方をノイマン境界条件と呼ぶことも多い．

一般化されたノイマン値とディリクレ境界条件を指定する．

Wolfram Language code: Subscript[Γ, N] = NeumannValue[-2u[x, y], ((x - 5) ^ 2 + (y - 5) ^ 2 == 3 ^ 2)];

Wolfram Language code: Subscript[Γ, D] = DirichletCondition[u[x, y] == 0, x == 0 && 8 ≤ y ≤ 10];

ノイマン境界方程式の項内の定数が，偏微分方程式(11)の弱形式からの内の定数と厳密にマッチすることを確かめることが重要である．これらの項はマッチしなければならないので，NeumannValueでは，境界方程式の右辺だけを入力すればよい．例えば，境界方程式を入力する場合には，NeumannValue[1,…]と入力するだけで十分である．境界方程式の左辺の正しい式（つまり）は，偏微分方程式(12)の弱形式における項によって示唆される．ノイマン境界値と偏微分方程式の間のこの密接な関係は，方程式のNeumannValue部分を作成することによって，さらに対応することができる．

偏微分方程式を解く．

Wolfram Language code: ufun = NDSolveValue[{op == Subscript[Γ, N], Subscript[Γ, D]}, u, {x, y}∈Ω]

解をプロットする．

Wolfram Language code: Plot3D[ufun[x, y], {x, y}∈Ω]

前の例では，すべての境界に境界条件が明示的に指定されてはいなかった．それらの部分の境界については，ノイマンのゼロ条件が使われる．

における解の値を調べる．

Wolfram Language code: Plot[ufun[0, y], {y, 0, 10}]

ディリクレ条件で必要とされるように，8を超える値はゼロであることに注意する．

周期境界条件を持つポアソン方程式

PeriodicBoundaryConditionは，境界の2つの独立した部分における偏微分方程式の解の値を関連付けるためのものである．PeriodicBoundaryCondition[a+b u,pred,f]は，解が独立したソースからの値に関連付けられるべきである点におけるターゲットを述部 pred が指し，関数 f がターゲット位置（青）から，従属変数関係が評価されるソース位置（緑）までの座標をマップする周期境界条件を指定する．結果の値は，ターゲットにおける条件として設定される．

PeriodicBoundaryConditionがどのように動作するのかを理解するのには，例を見るのが最もよい方法である．

青のターゲットと緑のソースの境界を示した領域を指定し，可視化する．

Wolfram Language code:

Ω = ImplicitRegion[1 ≤ x ^ 2 + y ^ 2 ≤ 2 ^ 2 && x - y ≥ 0 && x ≥ 0 && y ≥ 0, {x, y}];
Show[
	RegionPlot[Ω], 
	ContourPlot[y - x == 0, {x, 1 / Sqrt[2], 2}, {y, -0.1, Sqrt[2]}, ContourStyle -> Blue], 
	ContourPlot[y == 0, {x, 1, 2}, {y, -0.1, Sqrt[2]}, ContourStyle -> Green], AspectRatio -> Automatic]

この例では，偏微分方程式は，で指定される長方形内で1単位の負荷項がアクティブであるラプラシアンである．それ以外の場所では，負荷は0である．

偏微分方程式を指定する．

Wolfram Language code: pde = -Laplacian[u[x, y], {x, y}] == If[1.4 ≤ x ≤ 1.6 && 0 ≤ y ≤ 0.1, 1, 0];

PeriodicBoundaryConditionを使うと，周期境界条件，反周期境界条件，関係境界条件のいずれかをモデル化することができる．この例では，周期境界条件が示されている．周期条件は，第1引数が従属変数 u であり，この従属変数 u が述部（青，ターゲット）が位置（緑，ソース）においてと同じ値を持つようにしなければならない場合に指定される．解を求める過程で，この従属変数 u は，ソースにおいて評価される．このため，述部（青，ターゲット）から（緑，ソース）へのマッピングが必要となる．関数FindGeometricTransformは，座標のリストからそのようなマッピングを構築することができる．

FindGeometricTransformを使って，PeriodicBoundaryConditionのマッピングを計算する．

Wolfram Language code: mapping = Last[FindGeometricTransform[{{1, 0}, {2, 0}}, {{1 / Sqrt[2], 1 / Sqrt[2]}, {Sqrt[2], Sqrt[2]}}]];

マッピングが緑色のソースから青色のターゲットまでの点をマップするかどうかを調べる．

Wolfram Language code: mapping[{1., 0}]

その後，従属変数，ターゲット述部，マッピングの関係を与えることによって，PeriodicBoundaryConditionを設定する．周期条件の場合は，関係はである．

PeriodicBoundaryConditionを指定する．

Wolfram Language code: Subscript[Γ, P] = PeriodicBoundaryCondition[u[x, y], y - x == 0, mapping];

さらに，両方の曲線の辺に，DirichletConditionが置かれる．DirichletConditionは，ターゲットの境界上にDirichletConditionがないように，では除外されることに注意する．そこにディリクレ境界条件があると，それらの場所にマップされた解と整合性がなくなってしまう．

DirichletConditionを

において指定する．

Wolfram Language code: Subscript[Γ, D] = DirichletCondition[u[x, y] == 0, (x^2 + y^2 == 1 || x^2 + y^2 == 4) && !(y - x) > -10 ^ -7];

偏微分方程式を解く．

Wolfram Language code: ufun = NDSolveValue[{pde, Subscript[Γ, P], Subscript[Γ, D]}, u, {x, y}∈Ω];

解を可視化する．

Wolfram Language code: ContourPlot[ufun[x, y], {x, y}∈Ω, ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]

における解が

におけるものと同じであることを調べる．

Wolfram Language code:

imapping = InverseFunction[mapping];
Plot[{ufun[x, 0] - ufun@@imapping[{x, 0}]}, {x, 1, 2}, PlotRange -> All]

反周期境界条件は，PeriodicBoundaryConditionにおける従属変数の負を使うことにより実現できる．次の例では，反周期境界条件以外はすべて上と同じである．

反周期のPeriodicBoundaryConditionを

で指定する．

Wolfram Language code: Subscript[Γ, P] = PeriodicBoundaryCondition[-u[x, y], y - x == 0, mapping];

偏微分方程式を解く．

Wolfram Language code: ufun = NDSolveValue[{pde, Subscript[Γ, P], Subscript[Γ, D]}, u, {x, y}∈Ω];

解を可視化する．

Wolfram Language code: ContourPlot[ufun[x, y], {x, y}∈Ω, ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]

境界の2つの部分における解の総和を調べる．

Wolfram Language code:

imapping = InverseFunction[mapping];
Plot[{ufun[x, 0] + ufun@@imapping[{x, 0}]}, {x, 1, 2}, PlotRange -> All]

変数係数を持つ偏微分方程式

次の例では，偏微分方程式の係数は空間と時間に従属する．ここでは，は時間変数を，は空間変数を示す．

遮断されたマイクロストリップの例は，変数係数の効果を示す．これは，ポアソン方程式でモデルできる静電学の応用である．この方程式では，係数が相対誘電率であり，空間に依存する．

ここで，境界が2つの物質を分離する形状を構築する．それぞれの物質において，は異なる値を持つ．

まず，境界メッシュが作成される．要素メッシュの作成についての詳細は，ToBoundaryMeshを参照いただきたい．

境界メッシュを作成する．

Wolfram Language code:

sh = 2 / 10;sh2 = 2 / 100;sw = 3 / 10;
bmesh = ToBoundaryMesh[
	"Coordinates" -> {{0, 0}, {1, 0}, {1, sh}, {1, 1}, {0, 1}, {0, sh + sh2}, {sw, sh + sh2}, {sw, sh}, {0, sh}}, 
	"BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 5}, {5, 6}, {6, 7}, {7, 8}, {8, 9}, {9, 1}, {3, 8}}]}
	]

境界メッシュを可視化する．

Wolfram Language code: bmesh["Wireframe"]

境界メッシュから要素メッシュへの変換は，ToElementMeshを使って行う．

境界メッシュを完全なメッシュに変換し，結果のメッシュを可視化する．

Wolfram Language code:

mesh = ToElementMesh[bmesh];
mesh["Wireframe"]

次に偏微分方程式を設定する．偏微分方程式の負荷項は，電荷密度を真空誘電率で割ったものである．

形状の下の部分はシリコンであり，上の部分は真空である．シリコンの相対誘電率は約11.7であり，真空の値は1である．

Wolfram Language code: ϵr = If[y ≤ sh, {{11.7, 0.}, {0., 11.7}}, {{1., 0}, {0., 1.}}]

効率的な計算が行える偏微分方程式係数を設定する方法については，このチュートリアルも参照されたい．

電荷密度と真空誘電率は以下で与えられる．

Wolfram Language code:

ρc  = 10 ^ -8.;
ϵ0 = 8.86*^-12;

偏微分方程式の演算子を設定する．

Wolfram Language code: op = Inactive[Div][-ϵr.Inactive[Grad][u[x, y], {x, y}], {x, y}] - ρc / ϵ0

偏微分方程式でInactiveを使うことについては，形式的な偏微分方程式のセクションを参照のこと．

ポテンシャルは，である場合には必ず0に設定される．2つ目のポテンシャルは，内側のストリップで1000に設定される．

境界条件を設定する．

Wolfram Language code:

Subscript[Γ, D] = {DirichletCondition[u[x, y] == 0, x == 1 || y == 1 || y == 0], 
	DirichletCondition[u[x, y] == 10 ^ 3, 0 ≤ x ≤ sw && sh ≤ y ≤ sh + sh2]};

偏微分方程式を解く．

Wolfram Language code: ufun = NDSolveValue[{op == 0, Subscript[Γ, D]}, u, {x, y}∈mesh]

結果を可視化する．

Wolfram Language code: Show[ContourPlot[ufun[x, y], {x, y}∈mesh, ColorFunction -> "Temperature"], bmesh["Wireframe"]]

偏微分方程式係数に明示的な不連続性がある場合には，NDSolveが生成するメッシュにこれらの不連続性を自動的に含むことによって，問題を解決することが多い．における不連続性は，この例では自動的に検知される．

領域の記号的な指定を作成する．

Wolfram Language code:

sh = 2 / 10;sh2 = 2 / 100;sw = 3 / 10;
Ω = RegionDifference[Rectangle[{0, 0}, {1, 1}], Rectangle[{0, sh}, {sw, sh + sh2}]];
RegionPlot[Ω]

偏微分方程式を解く．

Wolfram Language code: ufun = NDSolveValue[{op == 0, Subscript[Γ, D]}, u, {x, y}∈Ω]

自動作成された境界メッシュを調べる．

Wolfram Language code: ToBoundaryMesh[ufun["ElementMesh"]]["Wireframe"]

見て明らかなように，もとの入力領域は不連続性に対して何も示さないが，メッシュ生成過程中に，不連続性に沿って内部に線が加えられている．

複数のものをモデリングする代りに，特定の領域部分にマーカーを使うこともできる．この方法については，「要素メッシュの生成」の「マーカー」セクションに詳しく説明されている．

その他にも，複数マテリアルの応用例が数多く「PDE Modelsの概要」ページに掲載されている．

また，効率的な偏微分方程式係数の設定についてのこの説明も参照されたい．

非線形係数を持つ偏微分方程式

偏微分方程式の係数の中には，空間および時間に加えて，従属変数および一次導関数に依存するものもある．ここでは，は空間変数 , , を指す．係数が従属変数に依存する場合には，その方程式は非線形である．以下の非線形方程式を見てみよう．

例として，極小曲面を計算するために，以下の方程式を使う．

従属変数の導関数が拡散係数内に現れるので，この方程式は非線形である．

円板を領域として設定する．

Wolfram Language code: Ω = Disk[];

偏微分演算子を設定する．

Wolfram Language code:

c = -1 / Sqrt[(1 + Grad[u[x, y], {x, y}].Grad[u[x, y], {x, y}])];
op = Inactive[Div][{{c, 0}, {0, c}} . Inactive[Grad][u[x, y], {x, y}], {x, y}]

非線形の定常方程式を解くために，解の初期シード値が指定する．これは，InitialSeedingsを使って行う．InitialSeedingsは，興味領域上で評価する式を指定する．この式は空間変数に依存することができる．初期シード値が指定されていない場合には，初期シード値は0であると想定される．

非線形方程式を解く．

Wolfram Language code: usol = NDSolveValue[{op == 0, DirichletCondition[u[x, y] == Sin[2π * (x + y)], True]}, u, {x, y}∈Ω]

結果を可視化する．

Wolfram Language code: Plot3D[usol[x, y], {x, y}∈Disk[]]

InitialSeedingの使用例は関数ページにも掲載されている．

その他にも，非線形の応用例が数多く「PDE Modelsの概要」ページに掲載されている．

非線形の変数係数を持つ偏微分方程式

以下のセクションでは，静磁場のアプリケーション例を示す．静磁場のアプリケーションでは通常シミュレーション範囲のさまざまな部分において，いくつかのマテリアル（非線形も可能）を利用する．この例では，磁位の分布と磁束のベクトル密度場が計算することで，モーターの簡単なモデルを見てみる．

上に表されるモーター形状では，固定子セクションが灰色，回転子が赤で表され，青い部分は空気を表している．またいくつかのコイルがグループになっていて，それぞれオレンジ色，薄いオレンジ色，黄色で表されている．

偏微分方程式モデルは，マクスウェル方程式のアンペールの法則から導き出されたものである．

ここでは電流密度，は透磁性である．磁束密度 iは，磁位に以下を通して接続されている．

1を2に代入すると，以下が得られる．

電流が平面と平面に直角にだけ流れると仮定して，モデルを2次元に簡約化し，を得る．さらにが方向に一定で，とが一定であると仮定すると，以下のように書くことができる．

モデルは方向に均一であると仮定されるので，部分を無視して，以下を使うことができる．

透磁性は空間座標と未知の導関数の関数であるので，方程式は非線形である．固定子と回転子では強磁性物質の透磁性が使われるのに対し，残りの領域では，空気の透磁性が使われる．方程式の右辺は，コイルを流れる変数電流密度の成分である．

まず最初にさまざまな範囲に対してマーカーがあるメッシュを作成する．これらのマーカーはその後，例えばシミュレーション範囲のさまざまな領域でをアクティベートするのに使うことができる．マーカーを使うこの方法は，偏微分方程式のさまざまな部分がアクティブである場合に領域について公式をしていするよりもずっと簡単である．マーカーとメッシュでのマーカーの生成についての詳細は，「要素メッシュの生成」のチュートリアルで説明されている．

手作業で生成されたメッシュ領域を使う．

Wolfram Language code: mr = [image];

MeshRegionから境界要素メッシュを生成する．

Wolfram Language code: bmesh = ToBoundaryMesh[mr]

マーカーがさまざまな部分領域に起因するものとするには，これらの部分領域にある座標を指定しなければならない．

マーカーの座標を設定する．

Wolfram Language code:

airCoordinate1 = {0., 0.};
rotorCoordinate = {0., 0.3125};
airCoordinate2 = {0, 0.8};
statorCoordinate1 = {0, 0.95};
fingerCoordinates = {{0.25, 0.6}, {0.6, 0.25}, {0.6, -0.25}, {0.25, -0.6}, {-0.25, -0.6}, {-0.6, -0.25}, {-0.6, 0.25}, {-0.25, 0.6}};
tipCoordinates = {{0.18, 0.44}, {0.44, 0.18}, {0.44, -0.18}, {0.18, -0.44}, {-0.18, -0.44}, {-0.44, -0.18}, {-0.44, 0.18}, {-0.18, 0.44}};
coilCoordinates = {{0.15, 0.64}, {0.35, 0.56}, {0.56, 0.35}, {0.64, 0.15}, {0.64, -0.15}, {0.56, -0.35}, {0.35, -0.56}, {0.15, -0.64}, {-0.15, -0.64}, {-0.35, -0.56}, {-0.56, -0.35}, {-0.64, -0.15}, {-0.64, 0.15}, {-0.56, 0.35}, {-0.35, 0.56}, {-0.15, 0.64}};

境界メッシュとマーカー座標を可視化すると，マーカー座標を指定する過程が理解しやすくなる．可視化のためには，マーカー座標をグループ化して，マーカーに使う色を設定しなければならない．

マーカー座標をグループ化し，マーカーの色を設定する．

Wolfram Language code:

markerColors = {Blue, Red, Gray, Orange};
markerCoordinates = {{airCoordinate1, airCoordinate2}, {rotorCoordinate}, {statorCoordinate1, fingerCoordinates, tipCoordinates}, {coilCoordinates}};

境界メッシュとマーカー座標をさまざまな色で可視化する．

Wolfram Language code:

Show[
	bmesh["Wireframe"], 
	Graphics[MapThread[{PointSize[0.02], #1, Point /@ #2}&, {markerColors, markerCoordinates}]]]

オレンジの点で可視化されているコイルは，さまざまな電流を伝導する．この過程を簡約化するために，同じ電流が流れるコイルを一緒のグループにするように，それぞれのコイルにマーカーが割り当てられる．

コイルのマーカーを指定する．

Wolfram Language code: coilMarkers = {4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 5, 5, 5, 5};

コイル座標と同じ数だけコイルのマーカーがあることを確認する．

Wolfram Language code: Length[coilMarkers] == Length[coilCoordinates]

これで部分領域にある座標のリストが得られ，これらの座標にマーカーが割り当てられる．

マーカーを指定する．

Wolfram Language code:

Short[markerSpecification = Join[{{airCoordinate1, 1}, {airCoordinate2, 1}, {rotorCoordinate, 2}, {statorCoordinate1, 3}}, {#, 3}& /@ fingerCoordinates, {#, 3}& /@ tipCoordinates, MapThread[{#1, #2}&, {coilCoordinates, coilMarkers}]]]

次に完全なメッシュを作成する．

マーカー付きのメッシュを生成する．

Wolfram Language code: mesh = ToElementMesh[bmesh, "RegionMarker" -> markerSpecification]

固定子が灰色，回転子が赤，空気が青，コイルがオレンジ色，薄いオレンジ色，黄色であるメッシュを可視化する．

Wolfram Language code: mesh["Wireframe"["MeshElementStyle" -> FaceForm /@ Join[markerColors, {LightOrange, Yellow}]]]

今度は物質データを設定する．

空気の透磁性を設定する．

Wolfram Language code: μAir = 4π * 10 ^ -7;

強磁性物質の透磁性は通常磁化（B-H）曲線で与えられる．これらの曲線には物質データが含まれる．このデータについて，偏微分方程式の係数として使える補間関数が作成される．より効率的な方法として，偏微分方程式の係数として使えるフィットが生成できる．両方の方法を以下に示す．

磁化曲線データシートから強磁性物質の値を指定する．

Wolfram Language code:

HField = {10, 25, 50, 75, 108, 137, 160, 180, 200, 235, 285, 320, 365, 480, 560, 660, 800, 1000, 1150, 1550, 1840, 2250, 2800, 3450, 4350, 5500, 7707, 11863, 15839, 21076, 27281, 35289, 45642, 59046, 76413, 98907, 128000, 134764, 141880, 149364, 157235};
BField = {0.0025, 0.0066, 0.0151, 0.0264, 0.05, 0.1, 0.15, 0.2, 0.25, 0.35, 0.5, 0.6, 0.7, 0.9, 1., 1.1, 1.2, 1.3, 1.354, 1.45, 1.5, 1.55, 1.6, 1.65, 1.7, 1.75, 1.81, 1.89, 1.945, 2., 2.05, 2.1, 2.15, 2.2, 2.25, 2.3, 2.35, 2.36, 2.37, 2.38, 2.39};

磁化曲線データを設定する．

Wolfram Language code: BHData = Transpose[{BField, BField / HField}];

補間関数が領域外でクエリされた場合には，補間を生成している間に，磁化曲線の最小値と最大値が継続するようにしなければならない．また，領域外でクエリを行う場合に補間関数が警告メッセージを発しないように指示しなければならない．

磁化曲線データから補間関数を生成する．

Wolfram Language code:

μB = With[{extrapolationTest = If[# < BField[[1]], Evaluate[BHData[[1, 2]]], Evaluate[BHData[[-1, 2]]]]}, 
	Interpolation[BHData, InterpolationOrder -> 2
	, ExtrapolationHandler -> {extrapolationTest&, "WarningMessage" -> False}
	]];

評価が速くできる偏微分方程式係数をどのようにして書けばよいかについての説明は，偏微分方程式係数の効率的な評価のセクションを参照されたい．

磁束密度のノルムを計算し，それがゼロにならないようにする．

Wolfram Language code: B2Norm = Sqrt[Total[Grad[u[x, y], {x, y}] ^ 2] + $MachineEpsilon];

透磁性係数が固定子または回転子で評価される場合には，空気の透磁性が使われるすべてのその他の領域で強磁性物質が使われる．

補間データに基づいて，逆透磁性

を指定する．

Wolfram Language code:

ν = Piecewise[{
	{-1 / μB[B2Norm], ElementMarker == 2 || ElementMarker == 3}
	}, -1 / μAir] * IdentityMatrix[2];

補間関数の方法は非常に柔軟なので，係数がそのまま使える．しかし，より効率的に評価できる方法は，磁化曲線データのフィットを求める方法である．

ガウスモデルで磁化曲線データのフィットを求め，関数を生成する．

Wolfram Language code:

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]]

補間された磁化曲線とフィットされた磁化曲線の差を可視化する．

Wolfram Language code: Plot[(μB[b] - fit[b]) / Max[BHData[[All, 2]]], {b, 0., 2.1}, PlotRange -> All]

効率的に評価可能な逆透磁性

を指定する．

Wolfram Language code:

ν = Piecewise[{
	{-1 / fit[B2Norm], ElementMarker == 2 || ElementMarker == 3}
	}, -1 / μAir] * IdentityMatrix[2];

効率的な計算が行える偏微分方程式係数を設定する方法については，このチュートリアルも参照されたい．

コイルの電流密度は，ElementMarker 4の要素内には正電流，ElementMarker 6の要素内には負電流，それ以外の場所には零電流が流れるように指定する．

電流密度を指定する．

Wolfram Language code:

jzCoil = 10;
jz = Piecewise[{{jzCoil, ElementMarker == 4}, {-jzCoil, ElementMarker == 6}}, 0];

偏微分方程式を指定する．

Wolfram Language code: pde = Inactive[Div][Inactive[Dot][ν, Inactive[Grad][u[x, y], {x, y}]], {x, y}] == jz;

固定子の外側の端において磁位がゼロであるように，境界条件を設定する．

境界条件を設定する．

Wolfram Language code: bcs = DirichletCondition[u[x, y] == 0, x ^ 2 + y ^ 2 ≥ 0.95];

非線形の偏微分方程式を解く．

Wolfram Language code: AbsoluteTiming[usol = NDSolveValue[{pde, bcs}, u, {x, y}∈mesh]]

磁位の最小値と最大値を求める．

Wolfram Language code: {minsol, maxsol} = MinMax[usol["ValuesOnGrid"]]

磁位と磁束密度Bを可視化する．

Wolfram Language code:

Show[
	ContourPlot[usol[x, y], {x, y}∈mesh, PlotRange -> All, Contours -> Range[minsol, maxsol, (maxsol - minsol) / 15], ColorFunction -> "TemperatureMap"], 
	bmesh["Wireframe"], 
	VectorPlot[Evaluate[{{0, 1}, {-1, 0}}.Grad[usol[x, y], {x, y}]], {x, y}∈mesh, StreamPoints -> Coarse, PlotTheme -> "Scientific"]
	]

InitialSeedingの使用例は関数ページにも掲載されている．

その他にも，非線形の応用例が数多く「PDE Modelsの概要」ページに掲載されている．

偏微分方程式と非線形境界条件

NDSolveおよび関連する関数では，非線形の一般化されたノイマン値を指定することができる．ディリクレ境界条件と周期境界条件は非線形にはならない．非線形の一般化されたノイマン値は以下で与えられる．

これは，線形の一般化されたノイマン値に非常によく似ているが，従属変数がの関数である点だけ異なる．

非線形の一般化されたノイマン値がよく使われるのは，放射境界をモデル化する際である．一次元の例として，非線形の放射熱移動が使われる．

右辺の境界条件は，非線形の放射条件である．

左辺の境界では，ディリクレ条件を利用する．ここで，は熱移動率係数，は放射にさらされた横断面積，はシュテファン・ボルツマン（Stefan–Boltzmann）定数，は放射率，は周囲温度，は熱源である．

非線形方程式を設定する．

Wolfram Language code: op = -κ * A * D[u[x], {x, 2}] + P * σ * ϵ(u[x]^4 - Subscript[T, ∞]^4) - q;

放射境界条件を設定する．

Wolfram Language code: Subscript[Γ, radiation] = NeumannValue[-ϵ * σ * A(u[x]^4 - Subscript[T, ∞]^4), x == 2 / 5];

ディリクレ条件を設定する．

Wolfram Language code: Subscript[Γ, D] = DirichletCondition[u[x] == Subscript[T, 0], x == 0];

物質パラメータを設定する．

Wolfram Language code:

parameters = {κ -> 40, σ -> 5675 / 1000 * 10 ^ -8, ϵ -> 35 / 100, q -> 0, Subscript[T, 0] -> 500, Subscript[T, ∞] -> 0, A -> π r^2, P -> 0 2π r, r -> 5 / 1000};

物質パラメータが代入された偏微分方程式を設定する．

Wolfram Language code: pde = {op == Subscript[Γ, radiation], Subscript[Γ, D]} //. parameters

方程式を解く．

Wolfram Language code: usol = NDSolveValue[pde, u, {x, 0, 2 / 5}, InitialSeeding -> {u[x] == Subscript[T, 0]} //. parameters]

結果を可視化する．

Wolfram Language code: Plot[usol[x], {x, 0, 2 / 5}]

類似の例を確かめるためには，「非線形有限法の検証テスト」を参照されたい．

熱方程式

熱方程式は，時間依存のポアソン方程式である．

従属変数は，空間座標と時間に依存する．

ディリクレ境界条件もノイマン境界値も時間に依存することがある．

偏微分方程式を解く全体的な手順は同じである．領域を指定し，境界条件を持つ偏微分方程式を設定しなければならない．

空間領域を設定する．

Wolfram Language code: Ω = ImplicitRegion[!((x - 5) ^ 2 + (y - 5) ^ 2 ≤ 3 ^ 2), {{x, 0, 5}, {y, 0, 10}}];

独立変数が時間と空間に依存するようになった．

時間

を

の引数として含むように，偏微分方程式の演算子を設定する．

Wolfram Language code: op = Subscript[∂, t]u[t, x, y] - Subsuperscript[∇, {x, y}, 2]u[t, x, y] - 20;

また，境界条件も空間と時間に依存することがある．

境界条件を設定する．

Wolfram Language code: Subscript[Γ, D] = DirichletCondition[u[t, x, y] == 0, x == 0 && 8 ≤ y ≤ 10];

Wolfram Language code: Subscript[Γ, N] = NeumannValue[-2 * u[t, x, y], (-5 + x) ^ 2 + (-5 + y) ^ 2 == 9];

一階の時間依存問題を解くには，NDSolveは初期条件と時間領域が必要である．

偏微分方程式を解く．

Wolfram Language code: ufunHeat = NDSolveValue[{op == Subscript[Γ, N], Subscript[Γ, D], u[0, x, y] == 0}, u, {t, 0, 100}, {x, y}∈Ω];

解のアニメーションを作成する．

Wolfram Language code:

numberOfFrames = 3;
frames = Table[Plot3D[ufunHeat[t, x, y], {x, y}∈ufunHeat["ElementMesh"], PlotRange -> {0, 160}], {t, 1, 20, (20 - 1) / numberOfFrames}];
ListAnimate[frames, SaveDefinitions -> True]

メッシュを一度変換し，離散化されたメッシュをプロット中に使用する方が，各アニメーションフレームについて領域を再度離散化するよりも，効率的であることに注意されたい．

熱移動方程式と適用可能な境界条件についての情報は，「熱移動」のチュートリアルと「熱移動の偏微分方程式と境界条件」のガイドページに記載されている．

波動方程式

波動方程式は，二階の時間依存方程式である．

円の領域

を定義する．

Wolfram Language code: Ω = RegionDifference[Disk[], Disk[{-1 / 4, 1 / 4}, 1 / 5]];

波動方程式の偏微分方程式演算子を定義する．

Wolfram Language code: op = D[u[t, x, y], {t, 2}] - D[u[t, x, y], {x, 2}] - D[u[t, x, y], {y, 2}];

ディリクレ境界条件は，すべての境界上で，固定値

に設定される．

Wolfram Language code: Subscript[Γ, D] = DirichletCondition[u[t, x, y] == 0, True];

波動方程式は，二階の時間依存偏微分方程式であり，初期条件と初期条件の導関数が指定されなければならない．

波動方程式の初期条件．

Wolfram Language code: ic = u[0, x, y] == Exp[-5((x - 0.2) ^ 2 + y ^ 2)];

Wolfram Language code: dic = Derivative[1, 0, 0][u][0, x, y] == 0;

偏微分方程式を解き，NDSolveに初期値を求めさせる．

Wolfram Language code: ufunWave = NDSolveValue[{op == 0, ic, dic, Subscript[Γ, D]}, u, {t, 0, 2π}, {x, y}∈Ω]

結果のInterpolatingFunctionのアニメーションの作成には，Plot3DとAnimateを使う．

解のアニメーションを作成する．

Wolfram Language code:

frames = Table[Plot3D[ufunWave[t, x, y], {x, y}∈ufunWave["ElementMesh"], PlotRange -> {-1, 1}], {t, 0, 2π, 2π / 20}];
ListAnimate[frames, SaveDefinitions -> True]

境界条件が指定されていない場合（陰的にこれはノイマン零境界条件を意味する）には，波動は境界において反射されている．これは1Dの例でよく見ることができる．

直線の領域

を定義する．

Wolfram Language code: Ω = Line[{{0}, {1}}];

波動演算子を定義する．

Wolfram Language code: op = D[u[t, x], {t, 2}] - D[u[t, x], {x, 2}];

波動方程式の初期条件．

Wolfram Language code:

u0[x_] = D[0.125 Erf[(x - 0.5) / 0.125], x];
ics = {u[0, x] == u0[x], Derivative[1, 0][u][0, x] == 0};

反射境界を持つ波動方程式を解く．

Wolfram Language code: ufun = NDSolveValue[{op == 0, ics}, u, {t, 0, 1}, {x}∈Ω]

InterpolatingFunctionの結果のアニメーションを作成するためには，Plot3DとAnimateを使う．

Wolfram Language code:

frames = Table[Plot[ufun[t, x], {x}∈Ω, PlotRange -> {-0.1, 1.3}], {t, 0, 1, 0.1}];
ListAnimate[frames, SaveDefinitions -> True]

時間依存の系には，境界条件は，方程式の従属変数の時間導関数の次数よりも厳密に小さい次数の従属変数の任意時間導関数に対して指定することができる．このため，時間において一次である方程式については，境界条件は従属変数についてだけ与えることができる．時間において二次である方程式については，境界条件は従属変数と，時間についてのその一次導関数に対して与えることができる．

吸収境界条件を持つ波動方程式をモデル化するためには，ノイマン境界条件の時間導関数を使う．これには，上で使われた波動方程式と全く同じものが使われるが，方程式にノイマン値が加えられる．

吸収境界を持つ波動方程式を解く．

Wolfram Language code: ufun = NDSolveValue[{op == -NeumannValue[Derivative[1, 0][u][t, x], x == 0 || x == 1], ics}, u, {t, 0, 1}, {x}∈Ω]

結果のInterpolatingFunctionのアニメーションを作成するためには，Plot3DとAnimateを使う．

Wolfram Language code:

frames = Table[Plot[ufun[t, x], {x}∈Ω, PlotRange -> {-0.1, 1.3}], {t, 0, 1, 0.1}];
ListAnimate[frames, SaveDefinitions -> True]

波動の周期境界条件をモデル化することは，PeriodicBoundaryConditionを使って実装することができる．この例では，初期条件の導関数はゼロではないが，初期の波動が右に移動するように選ばれる．

波動方程式の初期条件．

Wolfram Language code: ics = {u[0, x] == u0[x], Derivative[1, 0][u][0, x] == -D[u0[x], x]};

右の波動を吸収し，左に再挿入する周期境界条件を指定する．

Wolfram Language code: pbc = PeriodicBoundaryCondition[u[t, x], x == 0, TranslationTransform[{1}]];

吸収境界を持つ波動方程式を解く．

Wolfram Language code:

ufun = NDSolveValue[{op == -NeumannValue[Derivative[1, 0][u][t, x], x == 0 || x == 1], ics, pbc}, u, {t, 0, 2}, {x}∈Ω]

結果のInterpolatingFunctionのアニメーションを作成するためには，Plot3DとAnimateを使う．

Wolfram Language code:

frames = Table[Plot[ufun[t, x], {x}∈Ω, PlotRange -> {-0.1, 1.3}], {t, 0, 2, 0.1}];
ListAnimate[frames, SaveDefinitions -> True]

形式的な偏微分方程式

上の例では，偏微分方程式はLaplacianを利用して設定されている．そのまったく同じ偏微分方程式は，DivとGradを使って形式化することができる．それを行うには，で係数を使う．ここでは，× 行列（は内の次元の数）である．

Wolfram Language code: Div[-{{1, 0}, {0, 1}}.Grad[u[x, y], {x, y}], {x, y}] - 20

Inactiveを使うと，評価をしないようにすることができる．

偏微分方程式の演算子を設定する．

Wolfram Language code:

c = {{1, 0}, {0, 1}};
op = Inactive[Div][-c.Inactive[Grad][u[x, y], {x, y}], {x, y}] - 20

偏微分方程式の形式的なバージョンでは，係数行列を変更することができる．

Inactiveの使用が重要となってくる状況の一つとして，偏微分方程式とノイマン値の関係をモデル化する場合が挙げられる．以下の偏微分方程式を考える．

対応するノイマン境界方程式は以下の通りである．

が恒等行列であり，がで与えられる偏微分方程式を考える．Inactiveを使わずに偏微分方程式を指定すると，式が評価される．

Wolfram Language code:

c = IdentityMatrix[2];
α[x_, y_] = -E ^ (-x - y);
Div[-c.Grad[u[x, y], {x, y}] + Grad[-α[x, y], {x, y}] u[x, y], {x, y}]

ただし，評価によって偏微分方程式が以下のように変化することに注意されたい．

ここでは{-e^(-x-y),-e^(-x-y)}，は2e^(-x-y)である．

Wolfram Language code:

β[x, y] = D[-α[x, y], {{x, y}}];
a[x_, y_] = 2E ^ (-x - y);

これらの偏微分方程式は等しい．

Wolfram Language code:

Div[-c.Grad[u[x, y], {x, y}] + Grad[-α[x, y], {x, y}] u[x, y], {x, y}] == Div[-c.Grad[u[x, y], {x, y}], {x, y}] + β[x, y].Grad[u[x, y], {x, y}] + a[x, y] u[x, y]

偏微分方程式は等しいが，対応するノイマン値は異なる．

2つ目の偏微分方程式のノイマン値は以下のようになる．

評価を妨げ，以下の形式のノイマン値をモデル化するためには，Inactiveを使用することが必要である．

Wolfram Language code: Inactivate[Div[-c.Grad[u[x, y], {x, y}] + Grad[-α[x, y], {x, y}] u[x, y], {x, y}], Div | Grad];

形式的な偏微分方程式についての詳細は，「有限要素法の使用上のヒント」にある．

また，係数を関数として指定することは可能であるが，これを行うと，内部コードが係数積分を最適化することができなくなることに注意する．ここでは，係数は定数であり，最適化で係数は一度だけ積分される．ブラックボックス関数は，領域上で継続的に評価することが必要である．例えば，有限要素法のコードが最適化された拡散演算子を定数の係数に使うことができないので，内部を隠す関数は理想的ではない．

偏微分方程式の係数として避けるべき内部を隠しているブラックボックス関数．

Wolfram Language code: fun[x_ ? Real, y_ ? Real] := IdentityMatrix[2];

偏微分方程式系

偏微分方程式系の係数形式

従属変数とを持つ2つの方程式の系を考える．

内の領域については，各は，× 行列であり，，，のそれぞれはベクトルである．

この場合は，一般化されたノイマン方程式は以下の通りである．

一次元の結合偏微分方程式

次の例は，パラメータが簡単に変更できるように設定されている．まず，領域を指定する．

領域を指定する．

Wolfram Language code:

pos1 = 1;pos2 = 2;
Ω = ImplicitRegion[True, {{x, pos1, pos2}}]

以下は，解くべき偏微分方程式である．

偏微分方程式には，2つの従属変数とがある．結合は，2つの拡散演算子間である．これは一次元の系であるので，それぞれの係数は1×1の行列であり，これは数に等しい．

偏微分方程式の演算子を設定する．

Wolfram Language code:

Subscript[c, 11] = -2;Subscript[c, 12] = 10;Subscript[c, 21] = 3;Subscript[c, 22] = -1;
op1 = D[Subscript[c, 11] u[x], {x, 2}] + D[Subscript[c, 12]v[x], {x, 2}];
op2 = D[Subscript[c, 21] u[x], {x, 2}] + D[Subscript[c, 22] v[x], {x, 2}];

それぞれの方程式は，独自の境界条件を持つ．最も一般的な形式では，ディリクレ境界条件は条件の系である．

現行の実装は，ディリクレ条件の交差結合はサポートされていない．との項は与えることができない．

ノイマン値は，以下の形式の系を指定する．

次に，境界条件が設定される．

1つ目の方程式について境界条件を指定する．

Wolfram Language code:

Subscript[r, 1] = 1;
Subscript[Γ, D1] = DirichletCondition[u[x] == Subscript[r, 1], x == pos1];

Wolfram Language code:

Subscript[g, 1] = -2;Subscript[q, 11] = -1;
Subscript[Γ, N1] = NeumannValue[Subscript[g, 1] - Subscript[q, 11] u[x], x == pos2];

2つ目の方程式について境界条件を指定する．

Wolfram Language code:

Subscript[r, 2] = 1;
Subscript[Γ, D2] = DirichletCondition[v[x] == Subscript[r, 2], x == pos1];

Wolfram Language code:

Subscript[g, 2] = 3;Subscript[q, 22] = -4;
Subscript[Γ, N2] = NeumannValue[Subscript[g, 2] - Subscript[q, 22]v[x], x == pos2];

方程式系を解く．

Wolfram Language code:

{nufun, nvfun} = NDSolveValue[{op1 == Subscript[Γ, N1], op2 == Subscript[Γ, N2], Subscript[Γ, D1], Subscript[Γ, D2]}, {u, v}, {x}∈Ω];

以下は，DSolveを使って計算された，偏微分方程式の系の解析的な解である．

解析的な解．

Wolfram Language code:

{fufun, fvfun} = DSolveValue[{op1 == 0, op2 == 0, u[pos1] == Subscript[r, 1], v[pos1] == Subscript[r, 2], -Subscript[c, 11]u'[pos2] - Subscript[c, 12] v'[pos2] == Subscript[g, 1]  - Subscript[q, 11] u[pos2], -Subscript[c, 21] u'[pos2] - Subscript[c, 22] v'[pos2] == Subscript[g, 2] - Subscript[q, 22] v[pos2]}, {u, v}, x]

方程式系のDSolve形成において，一般化されたノイマン境界方程式がどのように明示的に与えられなければならないか，およびそれらの係数がどのように偏微分方程式の係数にマッチしているかについて注意する．

解析的な結果を数値結果と比べるために，2つの差分のプロットを作成する．

解析的な解と数値解の差分をプロットする．

Wolfram Language code: Plot[{fufun[x] - nufun[x], fvfun[x] - nvfun[x]}, {x, pos1, pos2}]

この例によって，なぜNeumannValueが方程式の一部であるのかも明らかになる．ノイマン値の方程式は以下で与えられる．

境界値を考慮する．

Wolfram Language code: Subscript[Γ, N1]

ノイマン値からは，どの方程式に関連付けられるべきであるかは分からない．ノイマン値を方程式の一部とすることによって，どれにノイマン値が関連付けられているかが明らかになる．

構造力学

構造力学では，負荷をかけられたときのオブジェクトの変形が計算される．構造力学は，よい結合偏微分方程式の例を提供する．領域の例として，5単位の長さで1単位の高さを持つ梁を考える．

物理をモデル化する偏微分方程式は，平面応力関係であり，薄いオブジェクトに有効である．ここで，「薄い」とは，オブジェクトの他の寸法に比べて薄いことを意味する．平面応力形成は，2つの従属変数と1つの独立変数に依存する．2つの従属変数は，変形したオブジェクトの成分規模での寄与を与える．物質データとして，ヤング(Young)の係数とポアソン率が指定されなければならない．

平面応力の演算子を定義する．

Wolfram Language code:

ClearAll[ν, Y]
op = {Inactive[Div][{{0, -((Y*ν)/(1 - ν^2))}, {-((Y*(1 - ν))/(2*(1 - ν^2))), 0}} . 
  Inactive[Grad][v[x, y], {x, y}], {x, y}] + Inactive[Div][{{-(Y/(1 - ν^2)), 0}, {0, -((Y*(1 - ν))/(2*(1 - ν^2)))}} . 
  Inactive[Grad][u[x, y], {x, y}], {x, y}], Inactive[Div][{{0, -((Y*(1 - ν))/(2*(1 - ν^2)))}, {-((Y*ν)/(1 - ν^2)), 0}} . 
  Inactive[Grad][u[x, y], {x, y}], {x, y}] + Inactive[Div][{{-((Y*(1 - ν))/(2*(1 - ν^2))), 0}, {0, -(Y/(1 - ν^2))}} . 
  Inactive[Grad][v[x, y], {x, y}], {x, y}]} /. {Y -> 10 ^ 6, ν -> 33 / 100};

演算子は，Inactive形式の形式的な偏微分方程式として与えられる．

あるいは，関数SolidMechanicsPDEComponentを使って偏微分方程式成分を生成することもできる．SolidMechanicsPDEComponentの関数ページにもテキスト形式の方程式が含まれている．

SolidMechanicsPDEComponentを使って，平面応力演算子を作成する．

Wolfram Language code:

SolidMechanicsPDEComponent[{{u[x, y], v[x, y]}, {x, y}}, <|"ModelForm" -> "PlaneStress", "YoungModulus" -> 10 ^ 6, "PoissonRatio" -> 33 / 100, "Thickness" -> 1|>]

さらに，境界負荷も指定されることがある．この場合は，境界負荷は梁の右端に加えられ，単位の圧力が負の方向に加えられる．

領域の左側では，どちらの従属変数も0の変位で固定されている．

境界条件を設定し，梁が左で固定されるようにする．

Wolfram Language code: Subscript[Γ, D] = DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, x == 0];

方程式を解く．

Wolfram Language code: {ufun, vfun} = NDSolveValue[{op == {0, NeumannValue[-1, x == 5]}, Subscript[Γ, D]}, {u, v}, {x, 0, 5}, {y, 0, 1}];

変形した梁を示すためには，静止している梁について要素メッシュを作成しなければならない．この梁は，与えられた補間関数によって変形される．

静止している梁の境界と，梁の変形したメッシュを可視化する．

Wolfram Language code:

mesh = ufun["ElementMesh"];
Show[{
	mesh["Wireframe"[ "MeshElement" -> "BoundaryElements"]], 
	ElementMeshDeformation[mesh, {ufun, vfun}]["Wireframe"["ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]

ElementMeshDeformationの使用についての詳細は，「要素メッシュの可視化」を参照されたい．

以下の等高線プロットは，梁の変形を示す．従属変数は，方向の変形を表す．梁の右上部分において，正の方向にシフトしたことが分かる．右下の部分は負の方向に動く．

Wolfram Language code:

ContourPlot[ufun[x, y], {x, 0, 5}, {y, 0, 1}, ColorFunction -> "Temperature", AspectRatio -> Automatic, PlotLegends -> Placed[Automatic, Above]]

従属変数の等高線プロットは，方向の変形を表す．梁の右側は，負の方向に変形する．

Wolfram Language code:

ContourPlot[vfun[x, y], {x, 0, 5}, {y, 0, 1}, ColorFunction -> "Temperature", AspectRatio -> Automatic, PlotLegends -> Automatic]

ここで，2つの等高線プロットの色スケールが同じではないことに注意する．

分析を完全なものにするため，物理の平面ひずみをモデル化する偏微分方程式を提示する．平面ひずみの形成は，平面応力の形成とは対照的に，太いオブジェクトに有効である．ここで「太い」とは，オブジェクトの他の寸法に比べて太いことを意味する．平面ひずみの形成は，2つの従属変数と独立変数に依存する．2つの従属変数は，変形したオブジェクトの成分規模での寄与を与える．物質データとして，ヤング(Young)の係数とポアソン率が指定されなければならない．

平面ひずみの演算子を定義する．

Wolfram Language code:

op = {Inactive[Div][{{0, -((Y*ν)/((1 - 2*ν)*(1 + ν)))}, {-(Y/(2*(1 + ν))), 0}} . 
  Inactive[Grad][v[x, y], {x, y}], {x, y}] + Inactive[Div][{{-((Y*(1 - ν))/((1 - 2*ν)*(1 + ν))), 0}, {0, -(Y/(2*(1 + ν)))}} . 
  Inactive[Grad][u[x, y], {x, y}], {x, y}], Inactive[Div][{{0, -(Y/(2*(1 + ν)))}, {-((Y*ν)/((1 - 2*ν)*(1 + ν))), 0}} . 
  Inactive[Grad][u[x, y], {x, y}], {x, y}] + Inactive[Div][{{-(Y/(2*(1 + ν))), 0}, {0, -((Y*(1 - ν))/((1 - 2*ν)*(1 + ν)))}} . 
  Inactive[Grad][v[x, y], {x, y}], {x, y}]} /. {Y -> 10 ^ 6, ν -> 33 / 100};

ここでも，成分は関数SolidMechanicsPDEComponentを使って生成されている．

この場合の境界条件は，両側への指定された変位である．領域の左側では，どちらの従属変数も0の変位で固定されている．右側では，負の軸の方向に，1単位の変位が与えられ，方向の変位は固定されている．

境界条件を設定し，梁が左右で固定されるようにする．

Wolfram Language code:

Subscript[Γ, D] = {DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, x == 0], DirichletCondition[{u[x, y] == 0., v[x, y] == -1.}, x == 5]};

方程式を解く．

Wolfram Language code: {ufun, vfun} = NDSolveValue[{op == {0, 0}, Subscript[Γ, D]}, {u, v}, {x, 0, 5}, {y, 0, 1}];

静止している梁の境界と，梁の変形したメッシュを可視化する．

Wolfram Language code:

mesh = ufun["ElementMesh"];
Show[{
	mesh["Wireframe"[ "MeshElement" -> "BoundaryElements"]], 
	ElementMeshDeformation[mesh, {ufun, vfun}]["Wireframe"["ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]

以下の等高線プロットは，梁の追加的な変形を示す．従属変数は方向の変形を，従属変数は方向の変形を表す．

Wolfram Language code: ContourPlot[ufun[x, y] + vfun[x, y], {x, 0, 5}, {y, 0, 1}, ColorFunction -> "Temperature", AspectRatio -> Automatic]

固体力学と適用可能な境界条件についての詳細は，「固体力学」のモノグラフおよび「固体力学の偏微分方程式と境界条件」のガイドページに記載されている．

液体の流れ

液体の流れの線形モデルは，ストークス方程式である．三次元のストークス方程式を以下に挙げる．

ここで，，はそれぞれ，，，方向の液体の速度である．は圧力，は液体の動的粘度である．

この例では，狭まっている二次元の管を通る流れについて分析する．

狭まっている領域を指定し可視化する．

Wolfram Language code:

Ω = RegionUnion[Rectangle[{0, 0}, {1, 1 / 2}], Rectangle[{1, 1 / 10}, {2, 2 / 5}]];
RegionPlot[Ω, AspectRatio -> Automatic]

ストークスの演算子を粘性

について設定する．

Wolfram Language code:

ClearAll[μ]
op = {Inactive[Div][{{-μ, 0}, {0, -μ}} . Inactive[Grad][u[x, y], {x, y}], {x, y}] + p^(1, 0)[x, y], Inactive[Div][{{-μ, 0}, {0, -μ}} . Inactive[Grad][v[x, y], {x, y}], {x, y}] + p^(0, 1)[x, y], u^(1, 0)[x, y] + v^(0, 1)[x, y]} /. μ -> 1;

演算子は，Inactive形式の形式的な偏微分方程式として与えられる．

ここでは，は方向の液体の粘性を，は方向の粘性を表す．　は液体内の圧力を示す．

偏微分方程式を指定する．

Wolfram Language code: pde = op == {0, 0, 0};

方向の流入の粘性を境界条件を使って規定する．これは，領域左側のの従属変数についての流入の輪郭を指定することによって可能である．領域の底と上の部分では，流れの条件は指定されない．これは，これらの部分のとの成分を0に設定することによって可能である．流出条件は，領域の右側で設定する．ここでは，圧力は，任意に0と設定する．一般に，それぞれの従属変数はディリクレ条件またはロビン型のノイマンが指定されているべきである．ノイマン0のような純粋なノイマン値では不十分である．その場合，方程式の解は定数までのみ正しい．

境界条件の設定．

Wolfram Language code:

bcs = {DirichletCondition[
	{u[x, y] == 4 * 0.3 * y * (0.5 - y) / (0.41) ^ 2, v[x, y] == 0.}, x == 0.], DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, 0 < x < 2], 
	DirichletCondition[p[x, y] == 0., x == 2]};

粘性が圧力よりも高次で補間される場合には，安定した解が求まる．NDSolveを使うと，それぞれの従属変数の補間次数を指定することができる．ほとんどの場合，"InterpolationOrder"を指定する必要はないが，これは液体の流れに有限要素法を使うときには標準的な過程であり，Wolfram言語ではこのように行われる．この設定は，P2P1型あるいはQ2Q1型の要素と呼ばれるものに対応しており，これらの要素はTaylor–Hood（テイラー・フード）要素としても知られる．"InterpolationOrder"を使わないと，すべての変数に対して二次補間を使うことになるが，それでもよい．すべての従属変数に対して一次補間を使うと，誤った結果になることがあり，数値的に安定しないかもしれない．速度に二次，圧力に一次を指定すると，解の確度がよくなり，離散化された偏微分方程式モデルのサイズが少し小さくなる．

デフォルトよりも細かいメッシュを指定して方程式を解く．

Wolfram Language code:

{xVel, yVel, pressure} = NDSolveValue[{pde, bcs}, {u, v, p}, {x, y}∈Ω, Method -> {"FiniteElement", "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}, "MeshOptions" -> {"MaxCellMeasure" -> 0.0005}}];

有限要素法オプションの補間次数の使用については，「有限要素法の使用上のヒント」で説明されている．

方向の流れを可視化する．

Wolfram Language code: ContourPlot[xVel[x, y], {x, y}∈Ω, AspectRatio -> Automatic, ColorFunction -> "TemperatureMap"]

方向の流れを可視化する．

Wolfram Language code: ContourPlot[yVel[x, y], {x, y}∈Ω, PlotRange -> All, AspectRatio -> Automatic, ColorFunction -> "TemperatureMap"]

圧力を可視化する．

Wolfram Language code: ContourPlot[pressure[x, y], {x, y}∈Ω, AspectRatio -> Automatic, ColorFunction -> "TemperatureMap"]

流れ場を可視化する．

Wolfram Language code:

VectorPlot[{xVel[x, y], yVel[x, y]}, {x, y}∈Ω, AspectRatio -> Automatic, StreamPoints -> 6, StreamColorFunction -> "TemperatureMap", StreamColorFunctionScaling -> False]

流量の非線形モデルは，ナビエ・ストークス(Navier–Stokes)方程式である．ナビエ・ストークス方程式を以下に挙げる．

は密度である．

この例では，2Dの横断面をチェックすることによって，円柱の周りの流入を分析する．

円筒の横断面の周りの領域を指定して可視化する．

Wolfram Language code:

rules = {length -> 22 / 10, height -> 41 / 100};
Ω = RegionDifference[Rectangle[{0, 0}, {length, height}], Disk[{1 / 5, 1 / 5}, 1 / 20]] /. rules;
region = RegionPlot[Ω, AspectRatio -> Automatic]

ナビエ・ストークス演算子を粘性

と密度

に対して設定する．

Wolfram Language code:

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

演算子は，Inactive形式の形式的な偏微分方程式として与えられる．

もう一つの方法として，関数FluidFlowPDEComponentを使ってより一般的に使えるナビエ・ストークス偏微分方程式成分のバージョンを生成することも可能である．FluidFlowPDEComponentの関数ページにはテキスト形式での方程式も含まれている．

FluidFlowPDEComponentを使ってナビエ・ストークス演算子を作成する．

Wolfram Language code: FluidFlowPDEComponent[{{u[x, y], v[x, y], p[x, y]}, {x, y}}, <|"DynamicViscosity" -> 10 ^ -3, "MassDensity" -> 1|>]

ここでは液体の方向の速度，は方向の速度を表す．は液体中の圧力を表す．

偏微分方程式を指定する．

Wolfram Language code: pde = op == {0, 0, 0};

方向の流入速度は，境界条件で規定される．これは，境界の左側の従属変数についての流入統計データを指定して行う．領域の下と上の部分，および円筒の周りでは，流出入なしの条件が指定される．これは，これらの部分におけるとのコンポーネントをどちらも0に設定することによって行う．流出条件は領域の右側で設定される．ここでは，圧力が任意に0に設定される．

境界条件の設定．

Wolfram Language code:

bcs = {DirichletCondition[u[x, y] == 4 * 0.3 * y * (height - y) / height ^ 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;

円筒の前の点と円筒の後ろの点からの圧力差を計算したいとする．さらに再循環の長さを計算する．円筒のけん引力と揚力も考慮しなければならない．は円筒の端の座標，は再循環領域の端である．よい結果を得るためには，円筒の周りの領域を梨形にして，デフォルトよりも高い確度と精度をメッシュ生成に使う必要がある．

円筒の周りの梨形に細分化された領域を設定する．

Wolfram Language code:

refinementRegion = ImplicitRegion[a ^ 4(23 (y - 1 / 5)) ^ 2 + b ^ 2(5 / 2x - 2.05) ^ 3(2a + (5 / 2x - 2.05)) < 0, {{x, 0, length}, {y, 0, height}}] /. Flatten[{rules, a -> 1, b -> 5 / 2}];

細分化領域とシミュレーション領域

を可視化する．

Wolfram Language code: Show[RegionPlot[Ω], RegionPlot[refinementRegion], AspectRatio -> Automatic]

細分化の領域に基づいて細分化関数を作成する．

Wolfram Language code:

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

速度が圧力よりも高い次数で補間される場合には，安定した解が求まる．NDSolveでは，それぞれの従属変数の補間次数が指定できる．

含める点，より高い確度と精度，細分化関数を指定しながら，方程式を解く．

Wolfram Language code:

{xVel, yVel, pressure} = 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}}];

有限要素法のオプションの補間次数の使い方については，「有限要素法の使用上のヒント」に説明がある．

円筒の前と後ろの圧力差を計算する．

Wolfram Language code: pressure[0.15, 0.2] - pressure[0.25, 0.2]

予想される圧力差の値は，0.1172と0.1176の間の値である．

円筒の後ろの再循環の長さを求める．

Wolfram Language code: (x /. FindRoot[xVel[x, 0.2] == 0, {x, 0.251, 0.4}]) - 0.25

予想される再循環の長さの値は，0.0842と0.0852の間の値である．

円筒の周りのけん引力と揚力を計算する関数を設定する．

Wolfram Language code:

force = With[{umean = 0.2, radius = 1 / 20, ρ = 1, μ = 10 ^ -3, dux = D[xVel[x, y], x], duy = D[xVel[x, y], y], dvx = D[yVel[x, y], x], dvy = D[yVel[x, y], y]}, Function[X, Block[{x, y, at, nx, ny, fx, fy, p}, 
	{x, y} = X;
	p = pressure[x, y];
	at = ArcTan[x - 1 / 5, y - 1 / 5];
	nx = -Cos[at];
	ny = -Sin[at];
	fx = nx * p + μ * (-2 * nx * dux - ny * (duy + dvx));
	fy = ny * p + μ * (-nx * (dvx + duy) - 2 * ny * dvy);
	2 * {fx, fy} / (ρ * umean ^ 2 * 2 * radius)
	]]];

円筒上のけん引力と揚力を計算する

Wolfram Language code: {fdrag, flift} = NIntegrate[force[{x, y}], {x, y}∈Circle[{0.2, 0.2}, 1 / 20], AccuracyGoal -> 3, PrecisionGoal -> 3]

予想されるけん引力の値は5.57と5.59の間の値であり，揚力の値は0.0104と0.0110の間の値である．

速度の

と

のコンポーネントを

，圧力の等高線を黒い線で，速度ベクトル場と一緒に可視化する．

Wolfram Language code:

Show[
	ContourPlot[Sqrt[xVel[x, y] ^ 2 + yVel[x, y] ^ 2], {x, y}∈Ω, ColorFunction -> "TemperatureMap", Contours -> 4], ContourPlot[pressure[x, y] == #, {x, y}∈Ω, ContourStyle -> Directive[Black, Thickness[0.002]]]& /@ {0.01, 0.02, 0.05, 0.08}, 
	VectorPlot[{xVel[x, y], yVel[x, y]}, {x, y}∈Ω, VectorPoints -> 9, VectorSizes -> .25, VectorColorFunction -> None, VectorStyle -> Gray], 
	 AspectRatio -> Automatic, ImageSize -> Large]

非定常ナビエ・ストークス（Navier–Stokes）方程式を解くこともできる．ベクトル化した方程式を以下に挙げる．

これには，定常の場合と同じ形状が使われる．まず時間依存のナビエ・ストークス方程式を設定する．

非定常のナビエ・ストークス演算子を粘性

と密度

に対して設定する．

Wolfram Language code:

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

左の流入をモデル化したいとして，右の流出の圧力条件を設定し，残りのすべての壁が滑らない境界条件を持つとする。流入境界条件が初期条件と一致するように，時間の経過とともに滑らかに流入を増やすヘルパ関数を作成する．

時間の経過とともに流入境界条件を上昇させる関数．

Wolfram Language code:

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]

境界条件を設定する．

Wolfram Language code:

bcs = {DirichletCondition[u[t, x, y] == sf[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;

初期条件は0に設定される．

Wolfram Language code: ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0};

ナビエ・ストークス方程式は高指数微分代数方程式であるので，効率的に方程式を時間積分するようにオプションが指定される．時間積分器の最大差分次数は，拒絶されるステップ数を最小化するために減らされる．時間依存境界条件は，時間積分をより効率的に行う一定の方程式系を作成するために微分される．最後に，デフォルトメッシュよりも幾分細かいメッシュが使われる．使用するハードウェアによっては，この方程式を解くのに数分かかることもある．

ナビエ・ストークス方程式を時間積分する．

Wolfram Language code:

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;)];]

速度範囲の最大値と最小値を求める．

Wolfram Language code: {minX, maxX} = MinMax[xVel["ValuesOnGrid"]]

効率的な可視化のために，解からメッシュを取り出す．

Wolfram Language code: mesh = xVel["ElementMesh"]

速度の等高線プロットを作成する．

Wolfram Language code:

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, 3, 10, 1 / 12}];]

速度のコンポーネントを可視化する．

Wolfram Language code: ListAnimate[frames, SaveDefinitions -> True]

流量と適用可能な境界条件についての詳細は，層流のモノグラフと流量モデルのガイドページに記載されている．

次のステップとして，ナビエ・ストークス方程式を熱方程式と組み合せて，エネルギー輸送をモデル化する．

まず領域といくつかのパラメータを指定する．

パラメータ，および一部が切り取られた長方形の領域を指定する．

Wolfram Language code:

sizes = {length -> 2, height -> 1, thickness -> 1 / 6, offset -> 1 / 4};
Ω = RegionDifference[RegionDifference[Rectangle[{0, 0}, {length, height}], Rectangle[{length / 2 - offset - thickness, height / 3}, {length / 2 - offset, height}]], Rectangle[{length / 2 + offset, 0}, {length / 2 + offset + thickness, 1 - height / 3}]] /. sizes;

領域を可視化する．

Wolfram Language code: RegionPlot[Ω, AspectRatio -> Automatic]

プラントル数

とレイリー数

を指定する．

Wolfram Language code: parameters = {ν -> Sqrt[Pr / Ra], ϵ -> 1 / Sqrt[Pr * Ra]} /. {Pr -> 7.1, Ra -> 2 * 10 ^ 5};

ブシネスク近似を利用する熱方程式と組み合せる粘性ナビエ・ストークス方程式を設定する．ベクトル化した方程式を以下に挙げる．

ここでは温度，は拡散係数，は動粘性率である．

指定した材料パラメータを使う．

熱方程式と組み合せられたナビエ・ストークス方程式を設定する．

Wolfram Language code:

ClearAll[ν]
op = {
	u^(1, 0, 0)[t, x, y] + 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}] + p^(0, 1, 0)[t, x, y], 
	v^(1, 0, 0)[t, x, y] + 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}] + p^(0, 0, 1)[t, x, y] - T[t, x, y], 
	u^(0, 1, 0)[t, x, y] + v^(0, 0, 1)[t, x, y], 
	T^(1, 0, 0)[t, x, y] + Inactive[Div][(-ϵ)*Inactive[Grad][T[t, x, y], {x, y}], {x, y}] + {u[t, x, y], v[t, x, y]}.Inactive[Grad][T[t, x, y], {x, y}]} /. parameters;

境界条件と初期条件を設定する必要がある．

先ほども述べたように，成分はFluidFlowPDEComponentとHeatTransferPDEComponentの関数を使って生成することも可能である．

すべての境界壁上における速度について滑りなし境界条件を設定する．

Wolfram Language code: wall = DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, True];

基準圧力点を設定する．

Wolfram Language code: reference = DirichletCondition[p[t, x, y] == 0, x == 0 && y == 0];

左と右の壁の温度差を指定する．

Wolfram Language code: temperatures = {DirichletCondition[T[t, x, y] == 1, x == 0], DirichletCondition[T[t, x, y] == 0, x == length]};

境界条件のパラメータを置き換える．

Wolfram Language code: bcs = {wall, reference, temperatures} /. sizes;

システムが静止しているような初期条件を設定する．

Wolfram Language code: ic = {u[0, x, y] == 0, v[0, x, y] == 0, p[0, x, y] == 0, T[0, x, y] == 0};

時間積分の進行と，細分化されたメッシュを使い，速度とおよび温度を二次で，圧力を一次で補間しながら，偏微分方程式を解くのにかかる総時間を管理する．

方程式の時間積分を管理する．

Wolfram Language code:

Monitor[AbsoluteTiming[
	{xVel, yVel, pressure, temperature} = NDSolveValue[{op == {0, 0, 0, 0}, bcs, ic}, {u, v, p, T}, {x, y}∈Ω, {t, 0, 300}, 
	Method -> {
	"PDEDiscretization" -> {"MethodOfLines", 
	"SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.00125}, "InterpolationOrder" -> {u -> 2, v -> 2, p -> 1, T -> 2}}}}
	, EvaluationMonitor :> (currentTime = Row[{"t = ", CForm[t]}])];
	], currentTime]

最後に後処理を行う．

領域の境界の可視化を構築する．

Wolfram Language code:

bmesh = ToBoundaryMesh[xVel["ElementMesh"]];
bmesh["Wireframe"]

圧力分布を可視化する．

Wolfram Language code: Show[ContourPlot[pressure[150, x, y], {x, y}∈pressure["ElementMesh"], ...], bmesh["Wireframe"]]

温度分布を可視化する．

Wolfram Language code: Show[ContourPlot[temperature[150, x, y], {x, y}∈temperature["ElementMesh"], ...], bmesh["Wireframe"]]

速度場を可視化する．

Wolfram Language code:

Show[bmesh["Wireframe"], VectorPlot[{xVel[150, x, y], yVel[150, x, y]}, {x, y}∈xVel["ElementMesh"], AspectRatio -> Automatic, PlotRange -> All]]

温度の変化と速度の流線をアニメーションにする．

Wolfram Language code:

numberOfFrames = 3;
frames = Table[Show[
	bmesh["Wireframe"], 
	ContourPlot[temperature[t, x, y], {x, y}∈temperature["ElementMesh"], ...], 
	StreamPlot[{xVel[t, x, y], yVel[t, x, y]}, {x, y}∈xVel["ElementMesh"], ...]
	], {t, 0, 300, 300 / numberOfFrames}];ListAnimate[frames, SaveDefinitions -> True]

興味があれば，切り取る部分を変更して作った別の速度場を調べてみるとよい．

別の領域を試す．

Wolfram Language code:

RegionPlot[RegionDifference[RegionDifference[Rectangle[{0, 0}, {length, height}], Rectangle[{length / 2 + offset, height / 3}, {length / 2 + offset + thickness, height}]], Rectangle[{length / 2 - offset - thickness, 0}, {length / 2 - offset, 1 - height / 3}]] /. sizes, AspectRatio -> Automatic]

流量，ブシネスク(Boussinesq)近似，および適用可能な境界条件についての詳細は，層流のモノグラフと流量モデルのガイドページに記載されている．