応用数値解析特論第 13 回

(1)

応用数値解析特論第 13 ^回

〜定常Navier-Stokes^{方程式に対する}Newton^法〜

かつらだ

桂田

まさし

祐史

http://nalab.mind.meiji.ac.jp/~mk/ana/

2020^年12^月21^日

(2)

本日の講義内容、連絡事項

当初は、第13^回で定常Navier-Stokes^{方程式に対する}Newton^{法と非定常}

Navier-Stokes方程式に対する特性曲線有限要素法の解説をする予定で

あったが、真面目に説明を用意したところ、思いの外 Newton^法の説明が長くなってしまった。

しかし、この手の話は意外ときちんと説明されないことが多いので(というかきちんとした説明を目にしたことがない)、時間をかけて説明することも意味があるだろうと考え直した次第。

参考に出来る資料を知らないので、意外とてこずりました。

(4)

11.4 定常 Navier-Stokes 方程式に対する Newton 法

11.4.1例題プログラムの紹介

前回、定常Stokes方程式を有限要素法で解いた。今回は定常

Navier-Stokes方程式を有限要素法で解いてみよう。

例題プログラムは FreeFem++ 4.7に付属する

/usr/local/ff++/share/FreeFEM/4.7-1/examples/examples/NSNewton.edp (2012 January, Hecht)

である。

ターミナルでこうすればコピーできる

cp -p /usr/local/ff++/share/FreeFEM/4.7-1/examples/examples/NSNewton.edp .

これはマニュアル (Hecht [1]) ^の§2.12 に掲載されているプログラムとほぼ同じである (本来完全に同じにすべきであろう)^。

かつらだ桂田

まさし

(5)

11.4.2 領域

Ω^はCe+Beb+Beo+Beu−Cc で囲まれる領域とする。ただし、R= 5,L= 15, Ce : (x,y) = (Rcost,Rsint) (t∈[π/2,3π/2]),

Cc : (x,y) = ((cost)/2,(sint)/2) (t∈[0,2π]), Beb : (x,y) = (Lt^1.2,−R) (t∈[0,1]),

−Beu: (x,y) = (Lt^1.2,R) (t∈[0,1]), Beo : (x,y) = (L,t) (t∈[−R,R]).

図1:領域Ω (∂Ω =Ce+Beb+Beo+Beu−Cc)

(6)

11.4.2 領域

Γ1:=−Cc+Ce+Beb+Beu, Γ2:=Beo

とおく。Γ1∪Γ2=∂Ωである。

Bebは(x,y) = (Lt,R) (t∈[0,1])と書いても同じことであるが、FreeFem++ではtについて等分することを想定してプログラムを書いたため、この式を採用しているのであろう。(こういう工夫が実際に効果があるのかどうかは分からない。)

Beu はプログラムでは

border beu(tt=1, 0){real t=tt^1.2; x=t*L; y=R; label=1;}

として与えられている。FreeFem++では、tt=1,0と書くことで、パラメーターが1から0まで減少する曲線を表すことが出来る、ということであろう(こういうことが出来ることは、このプログラムを見て初めて知った)。

かつらだ桂田

まさし

(7)

11.4.2 領域

Bebは(x,y) = (Lt,R) (t∈[0,1])と書いても同じことであるが、FreeFem++ではtについて等分することを想定してプログラムを書いたため、この式を採用しているのであろう。(こういう工夫が実際に効果があるのかどうかは分からない。)

Beu はプログラムでは

(8)

11.4.2 領域

Bebは(x,y) = (Lt,R) (t∈[0,1])と書いても同じことであるが、FreeFem++ではtについて等分することを想定してプログラムを書いたため、この式を採用しているのであろう。(こういう工夫が実際に効果があるのかどうかは分からない。)

Beuはプログラムでは

かつらだ桂田

まさし

(9)

11.4.3 領域の三角形分割

三角形分割をする際に

Bei: (x,y) = (L/2,t) (t∈[−R/4,R/4]) という曲線(線分)も用いている。

Th=buildmesh(cc(-40)+ce(20)+beb(15)+beu(15)+beo(8)+bei(10));

図2:領域の三角形分割. 中央部にB_ei (朱色の線分)が見える

(10)

11.4.4 方程式

このプログラムが解く問題の方程式は(境界条件が、マニュアルにも、プログラムの注釈にも明記されていないが…)、多分、次のものである。

(u· ∇)u−ν△u+∇p= 0 in Ω, (1a)

∇ ·u= 0 in Ω, (1b)

u=U on Γ1\Cc, (1c)

u=0 onCc, (1d)

−pn+ν∂u

∂n =0 onBeo. (1e)

(判断の根拠: 1というラベルをつけたΓ1でδu=0としていることから、そこで Dirichlet境界条件をつけているのだろう。初期値u0は

u0=

(1,0)^⊤ (x²+y²>2) (0,0)^⊤ (x²+y²≤2)

であったから、δuもCc では0, Γ1\Cc =Ce+Beb+Beu では(1,0)^⊤である。) (1e)は自然境界条件である(弱形式にはexplicitに現れない— Poisson方程式の問題における同次Neumann境界条件のようなもの)。

一様流(u=U:= (1,0)^⊤)の中に置かれた円柱の周りの流れを求めよという問題である。前回講義スライド14ページ「もう1つの弱形式」と類似点が多い。見比べてみよう。

かつらだ桂田

まさし

(11)

11.4.4 方程式

(u· ∇)u−ν△u+∇p= 0 in Ω, (1a)

∇ ·u= 0 in Ω, (1b)

u=U on Γ1\Cc, (1c)

u=0 onCc, (1d)

−pn+ν∂u

∂n =0 onBeo. (1e)

u0=

(1,0)^⊤ (x²+y²>2) (0,0)^⊤ (x²+y²≤2)

であったから、δuもCc では0, Γ1\Cc =Ce+Beb+Beu では(1,0)^⊤である。)

(1e)は自然境界条件である(弱形式にはexplicitに現れない— Poisson方程式の問題における同次Neumann境界条件のようなもの)。

(12)

11.4.4 方程式

(u· ∇)u−ν△u+∇p= 0 in Ω, (1a)

∇ ·u= 0 in Ω, (1b)

u=U on Γ1\Cc, (1c)

u=0 onCc, (1d)

−pn+ν∂u

∂n =0 onBeo. (1e)

u0=

(1,0)^⊤ (x²+y²>2) (0,0)^⊤ (x²+y²≤2)

かつらだ桂田

まさし

(13)

11.4.4 方程式

(u· ∇)u−ν△u+∇p= 0 in Ω, (1a)

∇ ·u= 0 in Ω, (1b)

u=U on Γ1\Cc, (1c)

u=0 onCc, (1d)

−pn+ν∂u

∂n =0 onBeo. (1e)

u0=

(1,0)^⊤ (x²+y²>2) (0,0)^⊤ (x²+y²≤2)

一様流(u=U:= (1,0)^⊤)の中に置かれた円柱の周りの流れを求めよという問題である。

前回講義スライド14ページ「もう1つの弱形式」と類似点が多い。見比べてみよう。

(14)

11.4.4 方程式

(u· ∇)u−ν△u+∇p= 0 in Ω, (1a)

∇ ·u= 0 in Ω, (1b)

u=U on Γ1\Cc, (1c)

u=0 onCc, (1d)

−pn+ν∂u

∂n =0 onBeo. (1e)

u0=

(1,0)^⊤ (x²+y²>2) (0,0)^⊤ (x²+y²≤2)

一様流(u=U:= (1,0)^⊤)の中に置かれた円柱の周りの流れを求めよという問題である。

前回講義スライド14ページ「もう1つの弱形式」と類似点が多い。見比べてみよう。

かつらだ桂田

まさし

(15)

部分積分公式

今回弱形式を導くのに必要な部分積分公式は以下のものである。

∫

Ω

△u·v dx =

∫

∂Ω

∂u

∂n ·v dσ−

∫

Ω

∇u :∇v dx, (2)

∫

Ω

∇p·v dx =

∫

∂Ω

pv ·n dσ−

∫

Ω

p∇ ·v dx.

(3)

(自力で確認できると思われるが、桂田[2]の付録Bに一応書いてある。)

(16)

11.4.5 問題点と対処の方針

問題になることと、それへの対処について、次の2^{点を指摘しておく。}

(a) Navier-Stokes方程式は非線形方程式であるので、一度に解を求めるのが難しい。

Newton法を用いて解くことにする。

(b) 動粘性係数ν=₂₀₀¹ は素朴な数値計算にとっては、かなり小さい(“円柱”の直径1 であるから、これを代表的な長さ、∥U∥= 1を代表的な速さに選ぶと、νの逆数が Reynolds数Re に相当するので、Re = 200ということになり、これが実はかなり大きい、ということである)。Newton法の初期値は真の解に十分近いものを選ばないと、反復が収束しないかもしれない。そのため、大きめのν^{についての解を求} め、それを小さいν の問題のNewton法の初期値とすることで、解きたいνに近づいていく。

かつらだ桂田

まさし

(17)

11.4.5 問題点と対処の方針

(b) 動粘性係数ν=₂₀₀¹ は素朴な数値計算にとっては、かなり小さい(“円柱”の直径1 であるから、これを代表的な長さ、∥U∥= 1を代表的な速さに選ぶと、νの逆数が Reynolds数Re に相当するので、Re = 200ということになり、これが実はかなり大きい、ということである)。Newton法の初期値は真の解に十分近いものを選ばないと、反復が収束しないかもしれない。そのため、大きめのν^{についての解を求} め、それを小さいν の問題のNewton法の初期値とすることで、解きたいνに近づいていく。

(18)

11.4.5 問題点と対処の方針

(b) 動粘性係数ν=₂₀₀¹ は素朴な数値計算にとっては、かなり小さい(“円柱”の直径1 であるから、これを代表的な長さ、∥U∥= 1を代表的な速さに選ぶと、νの逆数が Reynolds数Re に相当するので、Re = 200ということになり、これが実はかなり大きい、ということである)。Newton法の初期値は真の解に十分近いものを選ばないと、反復が収束しないかもしれない。そのため、大きめのνについての解を求め、それを小さいν の問題のNewton法の初期値とすることで、解きたいνに近づいていく。

かつらだ桂田

まさし

(19)

11.4.6 Newton 法の復習

方程式f(u) = 0の近似解を求めるためのNewton法とは、まず適当な初期値u0 を選び、漸化式

(4) ui+1=ui−f^′(ui)⁻¹f(ui)

により{ui}^{を定め、十分大きな}i に対して、ui を近似解に採用する、というものである。

(4)は、f(u) = 0をu=ui で線形近似した

f^′(ui)(u−ui) +f(ui) = 0 の解をui+1 とする、ということである。

wi :=f^′(ui)⁻¹f(ui)とおくと、wi はf^′(ui)wi =f(ui)の解であることに注意すると、次のアリゴリズムを得る(FreeFem++のマニュアルでは、F の微分はDF と表している)。

Findu∈V such thatf(u) = 0 wheref:V →V.

1. chooseu0∈V

2. for (i=0; i<niter; i=i+1)

1. Solvef^′(ui)wi =f(ui)

2. ui+1:=ui−wi

break∥wi∥< ε (修正量が小さくなったら反復を停止する)

(20)

11.4.6 Newton 法の復習

(4) ui+1=ui−f^′(ui)⁻¹f(ui)

1. chooseu0∈V

1. Solvef^′(ui)wi =f(ui)

2. ui+1:=ui−wi

かつらだ桂田

まさし

(21)

11.4.6 Newton 法の復習

(4) ui+1=ui−f^′(ui)⁻¹f(ui)

1. chooseu0∈V

1. Solvef^′(ui)wi=f(ui)

2. ui+1 :=ui−wi

(22)

11.4.7 我々の問題での f と f

^′

は？

f = 0が(1a) –(1e)と一致するようにするには

f(u,p) :=







(u· ∇)u−ν△u+∇p

∇ ·u u|Γ₁−ub

−pn+ν^∂u_∂n

B_eo







とすれば良いだろう。ここでub は次式で定まるΓ1上の関数である。

ub:=

(1,0)^⊤ (on Γ1\Cc) (0,0)^⊤ (on Cc).

かつらだ桂田

まさし

(23)

11.4.7 我々の問題での f と f

^′

は？ ( 続き )

f(u+δu,p+δp)−f(u,p)

=







((u+δu)· ∇)u−ν△(u+δu) +∇(p+δp)−((u· ∇)u−ν△u+∇p)

∇ ·(u+δu)− ∇ ·u (u+δu)|Γ₁−ub−

u|Γ₁−ub

−(p+δp)n+ν^∂(u+δu)_∂n

Beo−

−pn+ν^∂u_∂n

Beo







=







(δu· ∇)u+ (u· ∇)δu−ν△(δu) +∇(δp) +δu· ∇(δu)

∇ ·(δu) δu|Γ1

−(δp)n+ν^∂(δu)_∂n

Beo





.

であるから

f^′(u,p) δu

δp

=







(δu· ∇)u+ (u· ∇)δu−ν△(δu) +∇(δp)

∇ ·(δu) δu|_Γ₁

−(δp)n+ν^∂(δu)_∂n

B_eo





.

(書くと手間だが、要するに2次式であるから、慣れると結果がすぐに分かる。)

(24)

11.4.8 Newton 反復の弱形式

f^′(ui,pi) δu

δp

=f(ui,pi)を具体的に書くと

(δu· ∇)ui+ (ui· ∇)δu−ν△(δu) +∇(δp) = (ui· ∇)ui−ν△ui+∇pi in Ω (5a)

∇ ·δu=∇ ·ui in Ω (5b)

δu= 0 on Γ1, (5c)

−(δp)n+ν∂(δu)

∂n =0 onBeo. (5d)

当たり前のことであるが、δu,δpについての線形方程式である。

最初の2つに試験関数(それぞれv とq)をかけて積分すると Z

Ω

[((δu· ∇)ui)·v+ ((ui· ∇)δu)·v−ν△(δu)·v+∇(δp)·v]dx

= Z

Ω

((ui· ∇)ui−ν△ui+∇pi)·v dx Z

Ω

q(∇ ·(δu))dx = 0.

かつらだ桂田

まさし

(25)

11.4.8 Newton 反復の弱形式

f^′(ui,pi) δu

δp

=f(ui,pi)を具体的に書くと

(δu· ∇)ui+ (ui· ∇)δu−ν△(δu) +∇(δp) = (ui· ∇)ui−ν△ui+∇pi in Ω (5a)

∇ ·δu=∇ ·ui in Ω (5b)

δu= 0 on Γ1, (5c)

−(δp)n+ν∂(δu)

∂n =0 onBeo. (5d)

当たり前のことであるが、δu,δpについての線形方程式である。

最初の2つに試験関数(それぞれv とq)をかけて積分すると Z

Ω

[((δu· ∇)ui)·v+ ((ui· ∇)δu)·v−ν△(δu)·v+∇(δp)·v]dx

= Z

Ω

q(∇ ·(δu))dx = 0.

(26)

11.4.8 Newton 反復の弱形式 ( 続き )

それぞれ部分積分して(境界積分を消去するのに、(5d)を用いる) Z

Ω

[((δu· ∇)ui)·v+ ((ui· ∇)δu)·v−ν∇(δu) :∇v−δp(∇ ·v)]dx (6a)

= Z

Ω

q(∇ ·(δu))dx= 0.

(6b)

これが弱形式である。

かつらだ桂田

まさし

(27)

11.4.9 プログラム中の弱形式解読

次のようなマクロが定義されている。

Grad(u1,u2) [dx(u1),dy(u1),dx(u2),dy(u2)]

これは∇u を4次元ベクトルにしたものである。これから

Grad(du1,du2)’*Grad(v1,v2)

は

dx(du1)*dx(v1)+dy(du1)*dy(v1)+dx(du2)*dx(v2)+dy(du2)*dy(v2)

となる。これは∇(δu) :∇v である。

注 FreeFem++では、’^{は転置を表す}(MATLABでも使われる、古い記法の利用)。ベクトルの内積

real[int] a=[1,2,3,4]; real[int] b=[-4,-3,-2,-1]; cout << a’*b << endl;

これで−20という値が表示される。

(28)

11.4.9 プログラム中の弱形式解読

これは∇u を4次元ベクトルにしたものである。

これから

は

かつらだ桂田

まさし

(29)

11.4.9 プログラム中の弱形式解読

は

(30)

11.4.9 プログラム中の弱形式解読

は

real[int] a=[1,2,3,4];

real[int] b=[-4,-3,-2,-1];

cout << a’*b << endl;

かつらだ桂田

まさし

(31)

11.4.9 プログラム中の弱形式解読 ( 続き )

一方

UgradV(u1,u2,v1,v2) [[u1,u2]’*[dx(v1),dy(v1)],[u1,u2]’*[dx(v2),dy(v2)]]

は

u1*dx(v1)+u2*dy(v1) u1*dx(v2)+u2*dy(v2)

,

すなわち

(u· ∇)v =

u1

∂

∂x1

+u2

∂

∂x2

v1

v2

= u1∂v₁

∂x₁ +u2∂v₁

∂x₂

u1∂v₂

∂x1 +u2∂v₂

∂x2

!

を意味する。ゆえに

UgradV(du1,du2,u1,u2)’*[v1,v2]

は((δu· ∇)u)·v を、

UgradV(u1,u2,du1,du2)’*[v1,v2]

は((u· ∇) (δu))·v を計算する。

(32)

11.4.9 プログラム中の弱形式解読 ( 続き )

一方

は

,

すなわち

(u· ∇)v =

u1

∂

∂x1

+u2

∂

∂x2

v1

v2

= u1∂v₁

∂x₁ +u2∂v₁

∂x₂

u1∂v₂

∂x1 +u2∂v₂

∂x2

!

を意味する。

ゆえに

は((δu· ∇)u)·v を、

かつらだ桂田

まさし

(33)

11.4.9 プログラム中の弱形式解読 ( 続き )

一方

は

,

すなわち

(u· ∇)v =

u1

∂

∂x1

+u2

∂

∂x2

v1

v2

= u1∂v₁

∂x₁ +u2∂v₁

∂x₂

u1∂v₂

∂x1 +u2∂v₂

∂x2

!

を意味する。ゆえに

は((δu· ∇)u)·v を、

(34)

11.4.9 プログラム中の弱形式解読 ( 続き )

int2d(Th)(

nu*(Grad(du1,du2)’*Grad(v1,v2)) + UgradV(du1,du2, u1, u2)’*[v1,v2]

+ UgradV( u1, u2,du1,du2)’*[v1,v2]

- div(du1,du2)*q - div(v1,v2)*dp - 1e-8*dp*q // stabilization term

)

は次式を意味する。

Z

Ω

ν∇(δu) :∇v+ ((δu· ∇)u)·v+ ((u· ∇) (δu))·v−q∇ ·(δu)−δp∇ ·v−10⁻⁸(δp)q dx

- int2d(Th)(

nu*(Grad(u1,u2)’*Grad(v1,v2)) + UgradV(u1,u2, u1, u2)’*[v1,v2]

- div(u1,u2)*q - div(v1,v2)*p - 1e-8*p*q

)

は次式を意味する。

−Z

Ω

ν∇u_i:∇v+ ((u_i· ∇)u_i)·v−q∇ ·u_i−pi∇ ·v−10⁻⁸piq dx

かつらだ桂田

まさし

(35)

11.4.10 初期値の選択

最初にν= 1

50 ^{について、}

u0(x,y) =

(1,0)^⊤ (x²+y²>2) (0,0)^⊤ (x²+y²≤2)

をNewton法の初期値としてNavier-Stokes方程式の解を求め、それをν を小さくした Navier-Stokes方程式に対するNewton法の初期値としている。

最終的には、ν= 1

200 ^の場合のNavier-Stokes方程式の解を求めている。

Newton法の収束が早ければ、νをより速く小さくしたり、Newton法が収束しなければ、

ν を少し(戻して)大きくしたり、細かい制御を行っている。余裕があればプログラムを読んでみよう。

それが筋が通ったものであるか分からないが「そこまでやるか」と感心した。こういうプログラムを公開してもらえるのは、大変ありがたいことである。

(36)

11.4.10 初期値の選択

最初にν= 1

50 ^{について、}

u0(x,y) =

(1,0)^⊤ (x²+y²>2) (0,0)^⊤ (x²+y²≤2)

かつらだ桂田

まさし

(37)

11.4.10 初期値の選択

最初にν= 1

50 ^{について、}

u0(x,y) =

(1,0)^⊤ (x²+y²>2) (0,0)^⊤ (x²+y²≤2)

(38)

11.4.11 余談 : Oseen 方程式

NSNewton.edpの弱形式を定義しているところでsolve Oseen()^{と書いてある。}

オゼーン

Oseen方程式は、Navier-Stokes方程式を一様な定常流の周りで線形近似したものである。無限遠で小さいが0ではない、一様な速度Uex = (U,0,0)^⊤を持つ流れを扱いたい。

v :=u−Uex

とおく(v はδu と書くのが良いかもしれない)。このとき (u· ∇)u=

(U+v1) ∂

∂x1

+v2

∂

∂x2

+v3

∂

∂x3

U+v1

v2

v3



=U∂v

∂x1

+ (v· ∇)v≒U∂v

∂x1

.

そこで近似方程式として

(7) ∂v

∂t +U∂v

∂x1

=−1

ρ∇p+ν△v

の採用が考えられる。この方程式が普通Oseen方程式と呼ばれる(と私は習った)。

Wikipedia等を見ると、もう少し意味を広げて解釈されることもあるようである。

しかし、Hechtのプログラムで、Oseenと言っているのは、かなりの拡張解釈であると思

われる(うるさい人がツッコミを入れそうだ)。

“Oseen近似”という検索語でヒットした玉田[3]は、Oseen方程式が研究された経緯が分かる解説である(と私は感じた)。

かつらだ桂田

まさし

(39)

11.4.11 余談 : Oseen 方程式

オゼーン

Oseen方程式は、Navier-Stokes方程式を一様な定常流の周りで線形近似したものである。

無限遠で小さいが0ではない、一様な速度Uex = (U,0,0)^⊤を持つ流れを扱いたい。 v :=u−Uex

とおく(v はδu と書くのが良いかもしれない)。このとき (u· ∇)u=

(U+v1) ∂

∂x1

+v2

∂

∂x2

+v3

∂

∂x3

U+v1

v2

v3



=U∂v

∂x1

+ (v· ∇)v≒U∂v

∂x1

.

(7) ∂v

∂t +U∂v

∂x1

=−1

ρ∇p+ν△v

(40)

11.4.11 余談 : Oseen 方程式

オゼーン

Oseen方程式は、Navier-Stokes方程式を一様な定常流の周りで線形近似したものである。

無限遠で小さいが0ではない、一様な速度Uex = (U,0,0)^⊤を持つ流れを扱いたい。

v :=u−Uex

とおく(v はδu と書くのが良いかもしれない)。

このとき

(u· ∇)u=

(U+v1) ∂

∂x1

+v2

∂

∂x2

+v3

∂

∂x3

U+v1

v2

v3



=U∂v

∂x1

+ (v· ∇)v≒U∂v

∂x1

.

(7) ∂v

∂t +U∂v

∂x1

=−1

ρ∇p+ν△v

かつらだ桂田

まさし

応用数値解析特論 第 13 回