複素関数と流体力学
桂田 祐史
2015 年 6 月 17 日, 2015 年 7 月 29 日
目次
1
はじめに
2
2
流体力学の方程式
2.1 連続の方程式 (質量保存) . . . . . . . . . . . . .
2.2 応力テンソルと圧力 . . . . . . . . . . . . . . .
2.3 運動方程式 (運動量の保存) . . . . . . . . . . . .
2.4 数学的チャレンジ: Euler 方程式、Navier-Stokes
2.5 状態方程式 . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
2
2
4
5
6
6
.
.
.
.
7
7
7
8
11
.
.
.
.
.
.
.
.
12
12
13
15
17
17
17
18
18
3
4
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
方程式は解けるか？
. . . . . . . . . . . .
非圧縮流体の渦なしの流れ
3.1 渦度と渦なしの流れ . . . . . . . . . . . . . . . . . . . . .
3.2 非圧縮流体の渦なしの流れとポテンシャル (3 次元の場合)
3.3 非圧縮流体の渦なしの流れとポテンシャル (2 次元の場合)
3.4 流線と等ポテンシャル線 . . . . . . . . . . . . . . . . . . .
Poisson 方程式・Laplace 方程式に対する数値計算法
4.1 はじめに . . . . . . . . . . . . . . . . . . . . . . .
4.2 差分法 . . . . . . . . . . . . . . . . . . . . . . . .
4.3 有限要素法 . . . . . . . . . . . . . . . . . . . . .
4.4 基本解の利用 . . . . . . . . . . . . . . . . . . . .
4.4.1 基本解とは . . . . . . . . . . . . . . . . .
4.4.2 Poisson 方程式の特解 . . . . . . . . . . . .
4.4.3 Green の積分公式 (続き) . . . . . . . . . .
4.4.4 基本解の方法 . . . . . . . . . . . . . . . .
A ベクトル解析駆け足の復習
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
20
1
はじめに
1
何か一つだけ覚えるとすると、次の 1 行になるだろう。
2 次元の縮まない (非圧縮) 流体の渦なしの流れ = 正則関数
複素関数論を用いて流体力学の問題を取り扱うことも可能であるし、その逆に複素関数の問題
を流体力学のイメージで考えることも出来る。
流体 (fluid) とは大まかに言って、気体や液体のように “流れるもの” である。圧縮性と粘
性という 2 つの性質の有無によって大きく分類される。
複素関数と流体力学がオーバラップするところを学ぶ際の参考書としては、まずは今井 [1]
をあげておく (実はあまり読みやすくはない…)。著者の今井功先生 (1914–2004) は、流体力
学の権威で、著書の [2] は流体力学の基本的な文献であるとされている。航空機の飛行に関す
る流体力学の研究で著名であり、「物理学の散歩道」を書いたロゲルギストのメンバーでもあ
る。佐藤の超関数の流体力学的解釈を述べた [3] も非常に興味深い著作である。
流体力学のテキストとしては、[2] 以外に、巽 [4] をあげておく。話題が豊富であることに
加え、基本的な部分の説明は特に良く行き届いている。
なお、複素関数とのからみでは、最近の研究を背景として書かれた新井 [5] もとても参考に
なる (読みやすい教科書的な記述で、しかも研究ネタが見つかるかもしれない)。
非圧縮流体の方程式の数学的取り扱いについては、岡本 [6] を見よ。
複素関数以外にベクトル解析を用いる。ベクトル解析は、歴史的には電磁気学や流体力学
を記述するために発達したもので、現在でもそれらへの応用を意識しつつ学ぶことは有益で
ある。
流体力学の方程式
2
いきなり連立偏微分方程式の話になる。
流体の状態は、ふつう次の 3 つ (2 つ？4 つ？) のものを求めて定まる。
• 速度 v(x, t) = (u(x, t), v(x, t), w(x, t)) = (v1 (x, t), v2 (x, t), v3 (x, t))
• 圧力 p(x, t)
• 密度 ρ(x, t)
• 温度 ... 今回は考えない。
2.1
連続の方程式 (質量保存)
連続の方程式と呼ばれる次式が成り立つ。
(1)
∂ρ
+ div(ρv) = 0.
∂t
2
証明
流体中の任意の領域 V について
∫
∫
d
ρ dx = −
ρv · n dσ
dt V
∂V
が成り立つので、(左辺に微分と積分の順序交換、右辺に Gauss の発散定理を用いて)
∫
∫
∂ρ
dx =
div (ρv) dσ.
V ∂t
V
これが任意の V について成り立つことから、(1) を得る。
(1) は積の微分法により、次のようにも書ける。
∂ρ
+ (v · ∇) ρ + ρ div v = 0.
∂t
D
∂
:=
+v·∇
Dt
∂t
(2)
で定義される Lagrange 微分 (物質微分) を用いると、次のようにも書ける。
Dρ
+ ρ div v = 0.
Dt
(3)
物質微分について
流体の流れに沿って運動する粒子の、時刻 t での位置を x(t) とするとき、
∑ ∂f
∂f
d
∂f
∂f
Df
ẋj (t) +
f (x(t), t) =
=
+ ẋ(t) · ∇f =
+ v · ∇f =
.
dt
∂xj
∂t
∂t
∂t
Dt
j=1
3
密度の Lagrange 微分が 0 であること、すなわち条件
(4)
Dρ
=0
Dt
が成り立つとき、流体は非圧縮である (縮まない) という。このとき、連続の方程式は
(5)
div v = 0
と同値である (逆にこの条件が成り立つとき、連続の方程式を仮定すると、Dρ/Dt = 0 が導
かれる)。この条件を非圧縮条件と呼ぶ。
時々、ρ が定数になることを非圧縮条件と考えている人もいるが、定義は上のようになる (ρ
が定数ならば非圧縮であるが、逆は必ずしも成り立たない)。場所ごとに濃度の異なる食塩水
の流れを考えると、食塩水は縮まないが、ρ は定数にはならない。非圧縮かつ (ある瞬間に) 空
間的に一様であれば、密度は定数となる。
3
2.2
応力テンソルと圧力
重力のような力は体積に比例するため、体積力と呼ばれるが、流体が接触することによって
及びす力などは、接触面の面積に比例するため、面積力と呼ばれる。面積力の場合は単位面積
当たりの力を応力 (stress) と呼ぶ。
応力は、もちろん場所 a に依存するが、それだけでなく面の方向と向きに依存する。しば
らく a は固定しておいて、応力が面の方向と向きにどのように依存するか考察しよう。外向
き単位法線ベクトルが n である面についての応力を p(n) と書くことにする。
(点 a において
座標平面に平行な平面 xi = ai を通して正の側が負の側に及ぼす単位面積
)
pi1
当たりの力を pi2 として、P := (pij ) とおく。P を応力テンソルと呼ぶ。定義から
pi3
pi1
pi2 = p(ei ) (i = 1, 2, 3)
pi3
であるが、実は任意の n に対して
(6)
p(n) = P n
である。
授業では、(6) が成り立つと言っただけで、その根拠を説明しなかった。こういうことが気
になる人には (気にして欲しい) 時間をかけて学ぶことを勧める。
問 1.
このことを示せ。
問 2.
応力テンソルは一般に対称テンソルであることを示せ (角運動量の保存による)。
流体があるとき、
(7)
E := (eij ),
1
eij :=
2
(
∂vj
∂vi
+
∂xj ∂xi
)
で定まる E をひずみ速度テンソル (変形速度テンソル) と呼ぶ。
(反省点: 図でも描いて意味を説明すれば良かった。)
等方的な流体の多くで、
(8)
P = −pI + 2µE,
が成り立つ (そうである)。ここで µ は粘性率 (粘性係数, viscosity) と呼ばれる非負定数であ
り、p = p(x, t) は圧力と呼ばれる。
以下この文書では、(8) が成り立つような流体について論じることにする。(µ が定数である
場合を Newton 流体、定数でない場合を非 Newton 流体と呼ぶそうで、我々は以下では Newton
流体だけを考察する、ということでもある。)
4
もしも µ = 0 あるいは v = 0 (静止流体) ならば、力は −pn であり、方向は面に垂直、大
きさは面の方向に依らず一定、となる。
µ = 0 である流体を完全流体 (perfect fluid) または非粘性流体 (inviscid fluid), µ ̸= 0 であ
る流体を粘性流体 (viscous fluid) と呼ぶ。
この項の内容には天下りの部分が多いが、より詳しい説明が読みたい場合は、巽 [4] を見よ。
2.3
運動方程式 (運動量の保存)
最初に流体を構成する粒子の加速度は ∂v/∂t ではなく、Dv/Dt であることを注意しておく。
流体の運動方程式は
(9)
ただし
Dv
1
= div P.
Dt
ρ
∂p
11
∂x211
div P = ∂p
∂x1
∂p31
∂x1
+
+
+
∂p12
∂x2
∂p22
∂x2
∂p32
∂x2
+
+
+
∂p13
∂x3
∂p23
∂x3
∂p33
∂x3
である。
証明 任意の領域 V において
∫
∫
∫
Dv
ρ
dx =
P n dσ =
div P dx
Dt
V
∂V
V
が成り立つことによる。成分ごとに Gauss の発散定理を用いている。
(等方的な) 完全流体では div P = −∇p であるから、運動方程式は
1
Dv
= − ∇p.
Dt
ρ
すなわち
(10)
∂v
1
+ (v · ∇) v = − ∇p.
∂t
ρ
これが完全流体の運動方程式として有名な Euler 方程式である。
(等方的な) 粘性流体では P = −pI + 2µE となることから、
div P = −∇p + µ (△ v + grad div v) .
さらに非圧縮の場合は、右辺の最後の項が 0 であることに注意すると
(11)
1
∂v
+ (v · ∇) v = − ∇p + ν △ v.
∂t
ρ
これが非圧縮粘性流体の運動方程式として有名な Navier-Stokes 方程式である。ただし
ν :=
5
µ
ρ
とおいた。この ν を動粘性率 (kinematic viscousity) と呼ぶ。
Euler 方程式も、Navier-Stokes 方程式も非線形微分方程式である (速度の物質微分 Dv/Dt =
∂v/∂t + (v · ∇)v が非線形であることに注意する)。Navier-Stokes 方程式の粘性項を落とした
ものが Euler 方程式になる。
一方、Navier-Stokes 方程式で非線形項 (v · ∇)v を落として線形近似した
∂v
1
= − ∇p + ν △ v.
∂t
ρ
(12)
を Stokes 方程式と呼ぶ (粘性項 ν △ v は線形である)。線形方程式であるので、数学的には
大分簡単化したことになる。乱暴な近似であるようにも感じられるが、流速が小さい場合は、
Navier-Stokes 方程式の良い近似になっていると言われている。
2.4
数学的チャレンジ: Euler 方程式、Navier-Stokes 方程式は解けるか？
一般に、未知関数の個数と方程式の個数が一致することが方程式が解けるための必要十分条
件というわけではないが (例外はいくつもあげられる)、以下のように考えると分かりやすい
と思われる。
ベクトルを成分ごとに分けて考えると、連続の方程式は 1 つの方程式で、運動方程式は 3 つ
の方程式である。速度と圧力で 4 つの未知関数であるから、密度が既知であるとすれば、連続
の方程式と運動方程式を連立すると、未知関数の個数と方程式の個数が 4 でつりあって、速
度と圧力が求まる可能性があると考えられる。
例えば 密度が定数の場合、連続の方程式は非圧縮条件となり、これと Euler の方程式、ま
たは Navier-Stokes 方程式を連立させ、適当な境界条件と初期条件を与えて、解を求めよ、と
いう数学的な問題 (初期値境界値問題) が得られる。
古典的な問題ではあるが、非線形偏微分方程式であるため難しく、現時点でも完全な解決に
は至っていない。これについては、岡本 [6] を見よ。
2.5
状態方程式
ρ が未知関数のときは、連続の方程式と運動方程式だけでは、方程式が不足していて解けな
い。そこで、以下に説明する状態方程式を補うことが多い。
p が ρ のみの関数である、というバロトロピー流体の仮定
p = f (ρ)
をおくことで解ける場合が多い。具体的には、例えば断熱変化を仮定すると、∃γ > 1 s.t.
ρ ∝ p1/γ .
等温変化を仮定すると、
ρ ∝ p.
6
より一般に
1
∇p = ∇w
ρ
を満たす w が存在するとき、isentropic であるといい、w をエンタルピー (enthalpy) と呼ぶ。
バロトロピー流体は isentropic である。
非圧縮流体の渦なしの流れ
3
(粘性はあっても良い。言い換えると完全流体でも、粘性流体でも、以下の話は通用する。以
下の議論を見ると、Laplace 方程式や Poisson 方程式の重要性が感じられると思われる。これ
については別の解説を用意する。)
3.1
渦度と渦なしの流れ
ω := rot v
うずど
を渦度 (vorticity) と呼ぶ。これは流体粒子の自転の角速度の 2 倍に等しい (そうである)。
ω = 0 であるとき、流体は渦なしであるという。
流体が渦なしであれば局所的に v のポテンシャル ϕ が存在する:
v = ∇ϕ.
3.2
非圧縮流体の渦なしの流れとポテンシャル (3 次元の場合)
流体が非圧縮で渦なしであり、領域全体で v のポテンシャル ϕ が存在すると仮定する (例
えば領域が単連結であれば、渦なしであることから、ポテンシャルの存在が導かれる)。ϕ を
速度ポテンシャルと呼び、流れはポテンシャル流であるという。
v = ∇ϕ,
i.e. u =
∂ϕ
,
∂x
v=
∂ϕ
,
∂y
w=
∂ϕ
.
∂z
一般に div ∇ = div grad = △ であるから、
△ ϕ = div grad ϕ = div v.
非圧縮条件 div v = 0 を用いると
(13)
△ ϕ = 0.
すなわち、非圧縮流体のポテンシャル流において、速度ポテンシャルは調和関数である。
まとめると
渦なしならば (局所的には) 速度ポテンシャルが存在し、
さらに非圧縮ならば速度ポテンシャルは調和関数である。
7
例えば、流体の占める領域 Ω の境界で v が分かっていれば、(∂ϕ/∂n = ∇ϕ · n より)
∂ϕ
= v · n (on ∂Ω)
∂n
(14)
であるから、Laplace 方程式の Neumann 境界値問題 (13), (14) が得られ、それを解けば ϕ が
求まる (v も求まる)。
ここでは、連続の方程式から導かれる非圧縮条件だけから (運動方程式を用いずに) v が求
まった。それが出来た理由は、ポテンシャルが存在する場合は、非圧縮条件は 1 つの未知関数
に関する 1 つの方程式 (13) に帰着されるから、と言えるであろう。
上のようにして、速度が求まったとして、残る圧力は運動方程式から求められると期待され
る。それについて、一つの定理を紹介しよう。
完全渦なし流体のポテンシャル流においては、(さらに適当な仮定を追加でおいて) 「一般
化された Bernoulli の定理」が成り立つ。ここでは簡単のため、
ρ が定数 (つまり一様で非圧
( )
縮) と仮定しよう。∇
p
ρ
= ρ1 ∇p であるから、Euler 方程式に代入して、
∂v
+ (v · ∇) v = −∇
∂t
ベクトル解析の公式
(
v × rot v = ∇
より
∂v
= −∇
∂t
(
( )
p
.
ρ
)
1 2
|v| − (v · ∇) v
2
1 2 p
|v| +
2
ρ
)
+ v × rot v.
v = ∇ϕ と、rot v = ω = 0 を代入すると
(
)
p
∂ϕ 1
2
+ |∇ϕ| +
∇
= 0.
∂t
2
ρ
これから t のみの関数 g = g(t) が存在して
[
]
∂ϕ 1
2
(15)
p = −ρ
+ |∇ϕ| + g(t) .
∂t
2
これを一般化された Bernoulli の定理と呼ぶ。(以上の議論は完全渦なしのバロトロピー流体
に拡張できる)
3.3
非圧縮流体の渦なしの流れとポテンシャル (2 次元の場合)
流体の流れが 2 次元的であるとは、v(x, t) = (u, v, w) が次の条件を満たすことを言う。
(a) z 成分である w は 0.
(b) z によらない。
8
すなわち
u(x, y, t)
v(x, y, z, t) = v(x, y, t) .
0
このとき渦度は
∂ u e wy − vz
0
1
∂x
∂
v e2 = uz − wx = 0 .
ω = rot v = ∂y
∂
vx − uy
vx − uy
w e3 ∂z
ゆえに
0
ω = 0 ,
ω
ω := vx − uy .
以下では、2 次元の速度場 v = (u(x, y), v(x, y)) について考えることにして、rot v := vx − uy
と定める (割と良くやる定義)。
三次元の場合と同様に、
• v がポテンシャルを持てば (v = ∇ϕ = (ϕx , ϕy ) を満たす関数 ϕ が存在すれば)、ω = 0.
• ω = 0 であれば、v は局所的にはポテンシャルを持つ。
• ω = 0 であり、考えている領域 Ω が単連結であれば、v は Ω 全体でポテンシャルを持つ。
ϕ のことを流体の速度ポテンシャルと呼ぶのは 3 次元と同様である。
△ ϕ = ϕxx + ϕyy = ux + vy = div v.
ゆえに流体が非圧縮であれば △ ϕ = 0.
ω = 0 であることを渦なしの定義とする本が多いが、実際には領域全体でポテンシャルが存
在することを仮定しているので、ポテンシャルが存在することと渦なしの定義と考えるのが良
いかもしれない。あるいはポテンシャルが多価関数になることを許して扱うことになる (そう
している場合もある)。
ところで
非圧縮 ⇔ div v = 0
⇔ ux + vy = 0
⇔ ux = (−v)y
⇔ 局所的に ∃ψ
s.t. ψx = −v,
一般に
(16)
ψx = −v,
9
ψy = u
ψy = u.
を満たす ψ を v = (u, v) の流れ関数 (stream function) と呼ぶ。このとき、流れ関数は次の
Poisson 方程式を満たす:
− △ ψ = − (ψxx + ψyy ) = − (−vx + uy ) = ω.
ゆえに渦なしの非圧縮流ならば、流れ関数は調和関数である:
△ ψ = 0.
さらに、このとき
(17)
ϕx = u = ψy ,
ϕy = v = −ψx .
が成り立つ。
これは
(18)
f (z) := ϕ(x, y) + iψ(x, y)
についての Cauchy-Riemann 方程式である。ゆえに f は正則関数である。
f を v の複素速度ポテンシャルという。
f = ϕ + iψ であるとき、
f ′ = ϕx + iψx = u − iv
であるから、
u = Re f ′ , v = − Im f ′ .
(
)
q
cos
θ
f ′ の極形式を f ′ = qe−iθ とすると、v =
. ゆえに q が速さ、θ が速度の方向となる。
q sin θ
このことを背景に、f ′ のことを複素速度と呼ぶ (共役複素速度と呼ぶ人もいるらしく、その
気持もわかるけれど、複素速度という方が普通である)。
まとめると、2 次元の流れについて
• 渦なしならば (局所的には) 速度ポテンシャルが存在する。
(逆に速度ポテンシャルが存在するならば渦なしである。)
• 非圧縮ならば (局所的には) 流れ関数が存在する。
(逆に流れ関数が存在するならば非圧縮である。)
• 非圧縮かつ渦なしならば (局所的には) 速度ポテンシャル ϕ, 流れ関数 ψ が存在して、と
もに調和関数である。ψ は ϕ の共役調和関数であり、f := ϕ + iψ は複素速度ポテンシャ
ルと呼ばれる正則関数となる。
• v = (u, v) の複素速度ポテンシャル f について、f ′ = u − iv.
注意 3.1 Ω が単連結でなくても、Ω 全体で一価関数の流れ関数が存在する場合もある。実際、
一価関数となるためには、Ω 内の任意の閉曲線 C に対して、
∫
∫
u dy − v dx =
ψx dx + ψy dy = 0
C
C
が成り立つことが必要十分である。(穴があると、その周りを一周積分すると 0 にならないか
もしれないが、いつも ̸= 0 というわけではない。)
10
問 3. 流れ関数が存在するとき、Euler 方程式は次のように書ける (p が消去できることに
注目)。
∂
△ ψ − J(ψ, △ ψ) = 0.
∂t
ただし
J(f, g) := fx gy − fy gx .
(流れ関数が存在するならば非圧縮なので、この方程式は、Euler 方程式と非圧縮方程式を連
立させた方程式と同値であることになる。)
3.4
流線と等ポテンシャル線
速度場や圧力が時刻によらない場合を定常解、定常流と呼ぶ (流体が止まっているのではな
い)。定常流を可視化するための最も適当な方法は流線を描くことである。
定義 3.2 時刻 t における流れの流線 (streamline) とは、流れの領域の中の曲線で、その
各点における接線ベクトルと、その点における流れの速度ベクトルの方向が一致するもの
をいう。式で表すと、曲線 s 7→ x(s) で、
x′ (s) ∥ v(x(s), t)
を満たすものが流線である。
• 流線は時刻ごとに異なる。
• 流線は粒子の軌道とは異なる概念であるが、定常流の場合は流線は粒子の軌道と一致
する。
• 他に流脈線というものもあるが省略する。(流れの実験で良く可視化されるが、これは一
般には流線とも、粒子の軌道とも一致しない。)
2 次元の場合、流れ関数の任意の等高線は、流線である。
流体内にある仮想的な曲線 (障害物にはならないという意味) の微小部分 dr = (dx, dy) を
通過する流量 (単位時間あたりの体積) は、v · n ds であり、
(
)( ) (
)
0 1
dx
dy
n ds =
=
.
−1 0
dy
−dx
ゆえに
v · n ds = −v dx + u dy = ψx dx + ψy dy = dψ.
従って曲線 C 全体を、左から右に通過する流量は、C の始点を P0 , 終点を P1 として
∫
∫
v · n ds =
dψ = ψ(P1 ) − ψ(P0 ).
C
C
11
ゆえに、2 つの流線ではさまれた領域を流れる流体の流量は、流線上の流れ関数の値の差で
ある。
単純閉曲線 C で囲まれる範囲内に連続の方程式が成り立たない領域あるいは点が存在する
場合、湧き出しかあるいは吸い込みがあることになる。このとき C から外に湧き出る流量は、
いわゆる流束積分
∫
v · n ds
C
に等しい。
領域 D の境界 ∂D に対して、
∫
∫
v · n ds =
∂D
div v dx
D
非圧縮流体であればこれは 0 である。
流体の占める領域内の任意の閉曲線 C に対して、
∫
Γ :=
v · dr
C
を閉曲線 C に沿う循環と呼ぶ。
Poisson 方程式・Laplace 方程式に対する数値計算法
4
4.1
はじめに
Poisson 方程式・Laplace 方程式の境界値問題は、素性の良い問題なので、色々な数値解法
が適用できる。
モデル問題として
(19)
− △u = f
(20)
u = g1 (on ΓD ),
∂u
= g2 (on ΓN ).
∂n
(21)
(in Ω),
簡単のための 1 次元版として、Ω = (0, 1), ΓD = {0}, ΓN = {1}, g1 (0) = α, g2 (1) = β の場
合は、
(P)
(22)
− u′′ (x) = f (x) (x ∈ (0, 1)),
(23)
u(0) = α,
(24)
u′ (1) = β
ポピュラーな差分法 (FDM), 有限要素法 (FEM) を手短に紹介してから、Laplace 方程式
(f = 0) の場合に基本解を利用する方法の解説を行う。
12
4.2
差分法
常微分方程式の初期値問題に対する Euler 法, Runge-Kutta 法などを学んだことがあれば、
理解しやすいであろう。
差分法 (finite difference method, FDM) は次の 2 つの考え方を用いる。
• 微分方程式に含まれる導関数を差分商で置き換えた差分方程式の解を近似解に採用する。
• 領域を格子に区切って、格子点上の値を求めることを目標にする。
(25)
f ′ (x) =
f (x + h) − f (x)
+ O(h) (h → 0).
h
(26)
f ′ (x) =
f (x) − f (x − h)
+ O(h) (h → 0).
h
f ′ (x) =
(27)
f ′′ (x) =
(28)
(29)
f ′′′ (x) =
(30) f ′′′′ (x) =
f (x + h) − f (x − h)
+ O(h2 ) (h → 0).
2h
f (x + h) − 2f (x) + f (x − h)
+ O(h2 ) (h → 0).
h2
f (x + 2h) − 2f (x + h) + 2f (x − h) − f (x − 2h)
+ O(h2 ) (h → 0).
2h3
f (x + 2h) − 4f (x + h) + 6f (x, y) − 4f (x − h) + f (x − 2h)
+O(h2 ) (h → 0).
h4
多変数関数の偏導関数はこれらを適当に組み合わせて近似する。例えば
△ u(x, y) =
u(x + hx , y) − 2u(x, y) + u(x − hx , y) u(x, y + hy ) − 2u(x, y) + u(x, y − hy )
+
+O(h2x +h2y ).
2
2
hx
hy
例 4.1 1 次元 Poisson 方程式
−u′′ (x) = f (x) (x ∈ (0, 1)),
u(0) = α,
N ∈ N, h =
u′ (1) = β.
1
,
N
xi = ih とおいて、途中略で
α
f (x1 )
U1
2 −1
f (x2 ) 0
U2
−1 2 −1
..
... ... ...
+ ...
... = h2
.
−1 2 −1 UN −1
f (xN −1 ) 0
2βh
f (xN )
UN
−2 2
13
という差分方程式が得られる。境界条件 u′ (1) = β は仮想格子点 xN +1 を導入して
UN +1 − UN −1
=β
2h
と近似した。
例 4.2 (2 次元 Poisson 方程式の同次境界値問題) Ω := (0, W ) × (0, H) における Poisson 方
程式の
−(uxx (x, y) + uyy (x, y)) = f (x, y) ((x, y) ∈ Ω),
u(x, y) = 0 ((x, y) ∈ ∂Ω).
Nx , Ny ∈ N として、
hx :=
W
,
Nx
hy :=
H
,
Ny
xi = ihx ,
yj = jhy
として u(xi , yj ) の近似値 Uij を求めることを目標とする。
f (x1 , y1 )
U11
U21
f (x2 , y1 )
..
.
..
.
f (xNx −1 , y1 )
UNx −1,1
f (x1 , y2 )
U12
f (x , y )
U
22
2 2
..
..
U =
.
.
, f =
f (x
U
Nx −1 , y2 )
Nx −1,2
.
..
..
.
f (x1 , yNy −1 )
U1,Ny −1
f (x2 , yNy −1 )
U2,Ny −1
..
..
.
.
f (xNx −1 , yNy −1 )
UNx −1,Ny −1
とおくと
(
)
1
1
INy −1 ⊗ 2 (2INx −1 − JNx −1 ) + 2 (2INy −1 − JNy −1 ) ⊗ INx −1 U = f .
hx
hy
例 4.3 例えば熱方程式
ut (x, t) = uxx (x, t)
は (x, t) = (xi , tn ) において、左辺を前進差分、右辺を中心差分で置き換えると
u(xi , tn+1 ) − u(xi , tn )
u(xi+1 , tn ) − 2u(xi , tn ) + u(xi−1 , tn )
+ O(∆t) =
+ O(∆x2 )
2
∆t
∆x
14
を得る。そこで
un − 2uni − unu−1
un+1
− uni
i
= i+1
∆t
∆x2
という差分方程式の解 uni を u(xi , tn ) の近似解とすることが考えられる。
n
n
).
+ Ui−1
Uin+1 = (1 − 2λ)Uin + λ(Ui+1
差分解 (差分方程式の解) が偏微分方程式の解に収束することを示すには、偏微分方程式の
問題ごとの分析が必要である。Poisson 方程式や熱方程式については、最大値原理の離散版が
成り立つことを利用した証明がある。
the Lax equivalence theorem (Lax の同等定理)
適切な (well-posed) 線形初期値問題の適合 (consistent) な差分スキームが収束するにはそ
れが (Lax-Richtmeyer の意味で) 安定であることが必要十分である。
普通の差分方程式を作ると適合性は明らかに成り立つので、後は安定性をチェックするだけ
で収束が証明できることになる。
• 差分方程式の導出はとりあえず簡単だが、差分解の収束や安定性の証明は、もとの偏微
分方程式の性質の証明と関連深い。(結局勉強をサボるのは難しい。)
• 多次元の場合、境界が曲がっている領域を扱うには工夫が必要になる。(四角い領域に写
像してから解く, 不等間隔差分を用いる)
有限要素法
4.3
Laplace 方程式の境界値問題に対する Dirichlet の原理は、Poisson 方程式に対しても拡張
できる。
これも簡単のため、1 次元版で論じる (多次元でも本質的な違いはない)。
(P)
(31)
− u′′ (x) = f (x) (x ∈ (0, 1)),
(32)
u(0) = α,
(33)
u′ (1) = β
H (0, 1) を 1 階の Sobolev 空間として
1
{
}
X := v ∈ H 1 (0, 1) | v(0) = 0 ,
{
}
Xg1 := v ∈ H 1 (0, 1) | v(0) = α
とおく。
上の境界値問題の解は、次の 2 つの問題 (W), (V) の解でもある。
まず弱形式 (weak from)、あるいは弱定式化 (weak formulation) した問題。
15
(W)
Find u ∈ Xg1 s.t.
∫
1
′
∫
′
1
f (x)v(x)dx + βv(1) (v ∈ X).
u (x)v (x)dx =
0
0
変分問題にしたもの。
(V)
Find u ∈ Xg1 s.t.
J[u] = min J[w].
w∈Xg1
ただし
1
J[u] :=
2
∫
1
∫
′
2
1
|u (x)| dx −
0
f (x)v(x) dx − βv(1).
0
(W) と (V) は同値な問題であり、常に一意的な解を持つことが分かる。
逆に f がある程度滑らかであれば、(W), (V) の解は (P) の解であることが示される。
そこで (P) を解く代わりに、(W) あるいは (V) を解くことを目指す (変分法の直接法と呼
ばれる解き方)。
{xi }N
i=0 を
0 = x0 < x 1 < · · · < x N = 1
を満たす数列として、
e := {v ∈ C([0, 1]) | v は小区間 [xi − 1, xi ] では 1 次多項式と一致 } ,
X
{
}
e
X̂ := v ∈ X | v(0) = 0 ,
{
}
e | v(0) = α
X̂g1 := v ∈ X
e の要素を区分 1 次多項
とおくとき、X を X̂ で、Xg1 を X̂g1 で置き換えた問題を考える。X
式と呼ぶ。
次の 2 つの問題は同値であり、常に一意的な解 û を持つ。それを近似解として採用する。
c
(W)
Find û ∈ X̂g1 s.t.
∫
1
′
′
∫
f (x)v(x)dx + βv(1) (v ∈ X̂).
û (x)v (x)dx =
b
(V)
0
1
0
Find û ∈ X̂g1 s.t.
J[û] = min J[w].
w∈X̂g1
16
実は {xi } が [0, 1] の N 等分点であるとき、有限要素解 û の xi での値は、差分解 Ui と一
致する。もちろん、いつもそうなるわけではない (もしそうならば、2 つの方法を考える意味
がない)。
有限要素法には以下の利点がある。
• 弱形式の議論を済ませてあれば、有限要素解の厳密解への収束の議論は簡単になる。
• 多次元問題の場合に、長方形領域以外でも、それほど苦労なく解析が可能である。
• プログラムの自動生成がしやすい。
4.4
4.4.1
基本解の利用
基本解とは
1 1
1
, 2 次元の場合 E(x) = − log |x| を − △ の基本解 (the
4π |x|
2π
fundamental solutionn) と呼ぶ。
E は原点以外では調和関数であることは簡単な計算で分かる:
3 次元の場合 E(x) =
△ E(x) = 0 (x ∈ Rn \ {0}).
さらに、超関数の言葉で言うと (原点まで込めて)
−△E = δ
(34)
が成り立つ。ここで δ は Dirac のデルタ関数である。
物理的な解釈: 原点に置かれた単位点電荷の作る電場のポテンシャルが E(x) である。
4.4.2
Poisson 方程式の特解
∫
E(x − y)f (y)dy
U (x) :=
Ω
とおくと、
−△U = f
が成り立つ (証明は結構難しい)。つまり U は Poisson 方程式の特解である。
v := u − U とおくと、△ v = 0 が成り立つので、f = 0 の場合の問題が一般に解ければ良
いことになる。(実際には U を計算することは難しいことが多く、数値計算向きではないかも
しれない。)
17
4.4.3
Green の積分公式 (続き)
Green の第 2 積分公式
∫
∂Ω
(
∂u
∂v
v
−u
∂n
∂n
)
∫
(v △ u − u △ v) dx
dσ =
Ω
を利用して、次の Green の第 3 積分公式を得る (証明の詳細は略するが、関数論の Cauchy の
積分公式の証明のように、x を中心とする球を除いた領域で Green の公式を適用してから、球
の半径を 0 に近づける)。
∫
∫
∫
u(x) (x ∈ Ω)
∂
∂u
1
(y)dσy −
u(y)
E(x−y)dσy =
− E(x−y) △ u(y)dy+
E(x−y)
u(x) (x ∈ ∂Ω)
2
∂n
∂n
y
y
∂Ω
Ω
∂Ω
0
(x ∈ Rm \ Ω).
(x ∈ ∂Ω の場合は左辺第 3 項は主値積分である。)
(a) u が △ u = 0 を満たすならば、x ∈ Ω に対して、
∫
∫
∂
∂u
(y)dσy −
u(y)
E(x − y)dσy .
u(x) =
E(x − y)
∂ny
∂ny
∂Ω
∂Ω
すなわち、u が調和関数であるとき、u, ∂u/∂n の ∂Ω での値が分かれば、u(x) の値がこの
式で求まることになる (正則関数の Cauchy の積分公式に似ていて、使いでのある公式)。
境界条件から半分は分かっているので、もう半分求めれば良いことになる。
以下は細かい話になるが: 例えば ΓN = ∅ のとき、
∫
∫
∂
∂u
1
(y)dσy −
g1 (y)
E(x − y)dσy .
g1 (x) =
E(x − y)
2
∂ny
∂ny
∂Ω
∂Ω
これから ∂Ω 上で ∂u/∂n を求めることが出来る。
(b) u が ∂Ω の近傍で 0 ならば、x ∈ Ω に対して、
∫
− E(x − y) △ u(y)dy = u(x).
Ω
この事実を超関数解釈すると − △ E(x − ·) = δ(· − x) となる。
4.4.4
基本解の方法
簡単のため、ΓN = ∅ の場合の Laplace 方程式の境界値問題
△ u(x) = 0 (x ∈ Ω),
u(x) = g1 (x) (x ∈ ΓD = ∂Ω)
を取り上げる。
18
Ω が滑らかな境界を持つ有界領域の場合に、この境界値問題が一意的な解を持つことは知
られている。
Ω が (円盤とか、長方形とか) 特別な形をしている場合に、解 u を求める公式はいくつか知
られているが、ここでは多くの場合に使える数値解法を紹介する。一見素朴であるが、多くの
場合に良好な近似解を得ることが出来る。
Ω の外部に Ω を「囲むように」点 {yj }N
j=1 を取り、
(N )
u
(x) =
N
∑
Qj E(x − yj )
j=1
とおく。ここで Q1 , Q2 , . . . , QN は未定係数である。これらが何であっても
△ u(N ) (x) = 0 (x ∈ Rn \ {y1 , y2 , . . . , yN })
が成り立つ。もし u(N ) (x) = g1 (x) (x ∈ ∂Ω) が成り立てば u(N ) は解である。さすがにそんな
都合の良いことはめったにおこらないが、多くの場合、境界 ∂Ω 上で選んだ点 x1 , . . . , xN に
対して
u(N ) (xi ) = g1 (xi ) (i = 1, 2, · · · , N )
を満たすように Qj を決めることが出来る。
u(N ) ≒ u
が成り立つことは期待できる。この近似解法を the method of fundamental solutions (基本解
の方法, fundamental solution method), あるいは代用電荷法 (charge simulation method) と
呼ぶ。
次のような利点がある。
• しばしば誤差の指数関数的減少
(∃C)(∃ρ ∈ (0, 1))(∀N ∈ N)
u − u(N ) ≤ CρN
が成り立つ (この場合、差分法や有限要素法と比較して効率が高くなりうる)。
• u(N ) 自身が調和であり、特に
N
x − yj
1 ∑
Qj
2 (2 次元の場合)
− 2π
|x
−
y
|
j
j=1
grad u(N ) (x) =
N
1 ∑
x − yj
−
Qj
4π
3 (3 次元の場合)
|x
−
y
|
j
j=1
のように grad が直接計算できる。
基本解の方法以外に、基本解を利用する方法として、境界要素法 (boundary element method,
BEM) がある。
19
付録
A
ベクトル解析駆け足の復習
さすがに線積分、面積分の定義の復習は省略する。
Ω を Rn の領域とする。Ω を定義域とする n 次元ベクトル値関数 f : Ω → Rn のことを Ω
上のベクトル場と呼ぶ。
f1
f2
f = . とするとき、
..
fn
n
∑
∂fj
div f :=
∂xj
j=1
とおき、div f を f の発散と呼ぶ。div f のことを ∇ · f とも書く。
“良い” 境界を持つ領域 Ω の閉近傍で C 1 級のベクトル場 f に対して、
∫
∫
div f dx =
f · n dσ
Ω
∂Ω
が成り立つ (Gauss の発散定理)。ただし n は、∂Ω 上の点における、外向き単位法線ベクト
ルを表すとする。
C 1 級の関数 F : Ω → R に対して、
∂F
∂x1
∂F
∂x2
grad F := .
..
∂F
∂xn
とおき、F の勾配 (gradient) と呼び、∇F とも書く。これはベクトル場である。
Ω 上のベクトル場 f に対して、
grad F = f
(in Ω)
を満たす関数 F : Ω → R が存在するとき、F を f のポテンシャルと呼ぶ。
f の C 2 級のポテンシャル F が存在するならば、
(35)
∂fi
∂fj
=
∂xj
∂xi
(in Ω; i, j = 1, 2, . . . , n)
が成り立つ。
逆に C 1 級のベクトル場 f が (35) を満たしていて、Ω が単連結ならば、f のポテンシャル
が存在する。実際、Ω 内に任意の選んだ定点 a を始点として、Ω 内の任意の点 x を終点とす
る、Ω 内の区分的に滑らかな曲線 Cx に対して
∫
f · dr
Cx
20
の値は、F : Ω → Rn が Cx の取り方によらずに定まり、
∫
F (x) :=
f · dr
Cx
とおくことで得られる関数 F : Ω → Rn が f の一つのポテンシャルになる。
Ω 自身が単連結でない場合も、任意の開球は単連結領域であるから、f が (35) を満たして
いれば、任意の点の十分小さな近傍でポテンシャルが存在することになる。
3 次元の C 1 級ベクトル場 f に対して、
∂f3
∂f2
∂ f e −
1
1
∂x
∂x3
2
∂x
∂
∂f1
∂f3
f2 e2 = ∂x
rot f := ∂y
−
∂x
3
1
∂
∂f
∂f
2
1
− ∂x2
f 3 e3 ∂x1
∂z
とおき、f の回転と呼び、∇ × f , curl f とも表す。rot f は 3 次元のベクトル場である。
3 次元の C 1 級ベクトル場 f に対して、(35) は
rot f = 0
と同値である。ゆえに 3 次元のベクトル場 f がポテンシャルを持てば、rot f = 0 が成り立
ち、逆に rot f = 0 が成り立つならば、単連結領域においてポテンシャルが存在する。
( )
u
2 次元のベクトル場 f =
に対しては、(35) は
v
vx − uy = 0
と同値であり、同様のことが成り立つ。
参考文献
[1] 今井功：複素解析と流体力学, 日本評論社 (1989).
[2] 今井功：流体力学 前編, 裳華房 (1973), 後編は書かれなかった。
[3] 今井功：応用超関数論 I, II, サイエンス社 (1981).
たつみともまさ
[4] 巽 友正：流体力学, 培風館 (1982).
[5] 新井紀夫：複素流体力学ノート, コロナ社 (2012).
[6] 岡本久：ナヴィエ・ストークス方程式の数理, 東京大学出版会 (2009).
21
© Copyright 2024 ExpyDoc
