渦度場の離散化（渦要素モデルと壁面渦導入法）

(1) 渦に関する諸定理^{[71], [72]}

① Kelvinの循環定理

保存力場における非粘性流体の運動を考える．流体と共に動く任意の閉曲線Cをと り，Cに沿う循環 の時間変化を調べる．

C u ds (2.3.81)

C上の点の座標をxで表し，C上の線要素ds＝dxが流体と共に運ばれることを考慮 すると， の時間変化は次式で書ける．

C C

C C

Dt d d D

Dt D Dt d

D Dt d

D Dt D

u x u x

x u

x u

(2.3.82)

ここで，上式の右辺第一項はEuler方程式から以下のようになる．

Dt p

Du 1grad

(2.3.83)

C C C C

p dp

d p Dt d

D

1 1 grad

x u x

(2.3.84)

ここに[ ]Cは閉曲線Cを一周したときの，[ ]内の値の差を表す．

また，式（2.3.82）の右辺第二項は以下のようになる．

C C C C C

u du

d Dt d D Dt

d D

2 2

2 1 2 1

u u

u x u x

(2.3.85)

ただし，u ²＝|u|²である．

したがって，Cに沿う循環 の時間変化は以下のようになる．

C

u p Dt

D ₂

2

1 (2.3.86)

速度uおよび圧力pは物理量として座標の一価関数であり，これらの量が連続であ ると仮定すれば上式の右辺はゼロとなり，以下の結果が得られる．

Dt 0

D (2.3.87)

const. (2.3.88)

すなわち，流体の連続的な運動においては，流体と共に動く任意の閉曲線に沿う循 環は時間的に不変である．これをKelvinの循環定理という．

② Helmholtzの渦定理

流れの中に一本の曲線をとり，曲線上の各点を通る渦線によって形成される曲面を 考える．これを渦面という．渦面の法線ベクトルをnとすれば，以下の式が成り立つ．

0 t, x ω

n (2.3.89)

ある時刻 t において一つの渦面 S を考え，この面が流体と共に運ばれる際の時刻 t

t

t における変化を考える．

渦面S上に任意の閉曲線Cをとれば，Cに沿う循環 C は，Stokesの定理とKelvin の循環定理を考慮して以下のようになる．

0

dS d C

C S C

n ω

s u

(2.3.90)

ここに，S (C )はCで囲まれた渦面Sの部分を表す．閉曲線Cが後の時刻t において 面S 上の閉曲線C になったとすれば，C に沿う循環 C は Kelvinの循環定理によ って以下のようになる．

0 C

dS C SC ω n

(2.3.91) 閉曲線C は元の閉曲線 C を適当に選ぶことによって面S 上で任意にとれるから，

面S 上の至るところで（ω n）となり，面S も渦面であることがわかる．

このようにして，Kelvin の循環定理からの一つの帰結として，一つの渦面は流体運 動を通じて常に一つの渦面として保たれることがわかる．渦管は一つの渦面にほかな らず，また渦管の強さは渦管にとっての不変量であることを考慮すれば，上の結果は 次の定理として述べることができる．

理想流体の連続的な運動においては，一つの渦管は常に一つの渦管として保たれ，

その渦管の強さ は時間的に一定である．

const. (2.3.92)

これを Helmholtz の渦定理という．渦糸（無限に細い渦管）については，上の定理

から，一本の渦糸は流体運動を通じて常に一本の渦糸として保たれる．渦糸の断面積 を とすれば，渦糸の強さ は時間的に一定であることが結論される．

const. (2.3.93)

③ Lagrangeの渦定理

渦糸に対するHelmholtzの渦定理から，ある時刻において渦なしω= 0であった流体 は，その後の時刻においても常に渦なしである．逆に，ある時刻において渦ありω 0 であった流体は，その後も渦ありであることが結論される．すなわち，理想流体の連 続的な運動においては，渦度は発生することも消滅することもない．これを Lagrange の渦定理という．

例えば，流体の運動がある時刻に静止状態uから出発したとすると，静止状態では いたるところでω= rot u=0であるから，その後の運動もいたるところで渦なしでなけ ればならない．同様に流体流れが無限上流で一様流u = const.であったとすれば，そこ

ではω= rot u = 0であるから，下流においても至るところで渦なしでなければならない．

(2) 渦要素のモデル化

① 二次元モデル

一般化されたBiot-Savartの式において，渦度の空間積分項を二次元で離散化すると，

以下のように表せる．

M i M i

i S dSi R

R

i 1 2

1 2 2

1 2

1 ω R k R

(2.3.94) このように，速度は自由渦分布からの誘起速度として表されることが分かる．

初期における離散渦法では，比較的高いレイノルズ数の流れを取り扱ってきたので，

粘性の影響は物体近傍およびはく離せん断層に限られてきた．そのため，粘性による 渦度の拡散は比較的小さく，流れ場の物体表面およびはく離せん断層を除く領域は非 粘性のポテンシャル流れと近似的にみなすことにより解析が行われてきた．この際に 導入される渦要素には，渦度が一点にデルタ関数的に存在する自由渦が用いられてき た．しかし，Biot-Savart の式からも明らかなように，このモデルでは微小渦要素近傍 や，渦要素の中にある位置での誘起速度の計算には不適切であり，計算精度の低下を 招く原因とも成り得た．

そこで，各渦要素に粘性核構造（渦コア半径）を設けて，その渦度分布を精度良く 積分するか，もしくは合理的な速度分布を仮定するかによって，その不適切さを取り 除く必要がある．以下に，渦モデルを二つ示す．

(a) コア付き自由渦点

各渦要素に半径 の粘性核構造を持たせ，その中では渦度 が一様に分布すると し， ²を循環量とするコア付き自由渦点モデルがある．渦要素からの誘起 速度は，次式に示すように渦中心からコア半径 までは強制渦型の速度分布を持ち，

それより外部では自由渦型の速度分布を持つものとする．

2 2

2

2 2

i i

i c

i i

v R

R u k

：R (2.3.95)

0 2

2

2 2

i

i c

i i

R R v

R u k

：R (2.3.96)

(b) Blobモデル

コア付き自由渦点モデルでは，コア半径 において誘起速度および渦度が不連続 となるため，実際の流れとの相違が生じる．そこで，誘起速度および渦度を連続的 に分布させるモデルとして，Leonard^[37]の提案した Gauss 関数型に渦度の分布する Blobモデルが通常よく用いられている．

渦要素からの距離をrとするとき，渦度 _zの分布は次式のように表される．

f A f

i i z i

i i

i i

z 2 2

1 r r

r (2.3.97)

ここに， Aiは渦要素の微小面積であり， _iは微小面積における渦度を示す．Γ_i はΓ_i= i Aiで定義される循環量を示す．f ( )は点対称な渦度の分布関数であり，

Gauss関数分布型の核構造を採用すれば，以下のようになる．

exp 2

f 1 (2.3.98)

任意の位置rにおける渦要素からの誘起速度uiは，Biot-Savart の法則より以下の ように算出される．

R g

i

i 2 2

1 R

u (2.3.99)

2 0

i i

exp 1 2 t f t dt g

R

i

R r r

r r R

上式から明らかなように，渦要素の持つ渦度分布は渦要素の中心に対して点対称 な粘性核半径 を分散とする Gauss 型分布を用いるため，渦要素は自分自身に対す る誘起速度を持たない．

本研究では，ここで示した渦Blobモデルを用いる．また，渦Blobの他に，渦度 が一様に分布する矩形渦要素（Vortex sheet）を用いる．

② 三次元モデル

従来研究において提案された三次元渦モデルは，渦糸モデル，Stick モデル，Vorton モデル，Blobモデルの４種類に大きく分類できる．本研究では，三次元渦モデルとし て，二次元同様のBlobモデルを用いる．

流れ場に存在する渦度からの誘起速度は，Biot-Savartの式より次式で表せる．

R dV

V 3

4

1 ω R

r

u (2.3.100)

このBiot-Savartの式の空間積分を離散化し，N個の渦要素からの誘起速度の和とし

て表せば次式ように表せる．

N

i R³

4

1 α R

r

u (2.3.101)

ここで，RおよびRは，riを渦要素の位置とし，次式で与える．

ri

r

R ，R r r_i (2.3.102)

また，次の関係を持つαを渦強度とする．

N

i i

VωdV α (2.3.103)

したがって，任意の渦要素は，渦強度ベクトルα_iと位置ベクトルr_iによって表され

る．式（2.3.101）による離散化では，渦要素の近傍において速度場が特異性を持つこ

とがわかる．渦法では，求められた速度場から渦要素の移動や渦度の変化を計算する

ことから，速度場を精度良く求めることが重要である．そこで，滑らかな速度分布を 得るために渦要素の粘性核構造として渦度の分布関数を導入する．

渦要素の持つ渦度の分布ω_i r は，分布関数 f ( )及び渦度の空間的広がりを表す粘 性核半径 _iを用いて次式で与える．

f f

i i i

i i

i

i 3 3

r α α r

r

ω (2.3.104)

i

ri

r (2.3.105)

したがって，N 個の渦要素があるとき，渦度の分布は個々の渦要素の持つ渦度の重 ね合わせにより次式で表せる．

N

i i

N i i

i α₃ f

r ω r

ω (2.3.106)

このように定義される渦要素によって，任意位置 r に誘起される速度ur は，

Biot-Savartの式により次式より得られる．

1 lim

4

4 1 4

1

0 2

3 3

g

dt t f t g

R g R

x

N i N i

i

i r R α R

r ω u

(2.3.107)

このとき，渦度の空間積分値と渦強度の間には，次に示すような関係式が成り立つ．

N

i i

N

i i i

i

V R dR

f R

dV α α

ω 0 3

4 2 (2.3.108)

本研究では，分布関数f ( )として次式に示すWincklemans&Leonard^[56]によって提案 された分布関数を用いる．

2 2 5

5 2 0

2 2 2 7

1 2

2 4 5

1 1 8 15

dt t f t g

f

(2.3.109)

(3) 渦要素の導入方法

本研究では，非定常なはく離や逆流を伴う内部流れの解析を行うため，壁面上に発生 する渦度の分布を全て離散渦要素で置き換え，それらの挙動をラグランジュ的に追跡す ることで，境界層はく離や非定常な逆流等のダイナミックな流れ現象を捉えられるよう にした．

① 二次元モデル

渦要素導入法の概略を図2.3.4に示す．壁面極く近傍（境界層排除厚さh）の速度を

Biot-Savart の式から算出し，その速度勾配から決定される一定渦強度の渦度を，壁面

上の薄い渦層内の領域（a-b-c-d）に分布させる．

渦層内の速度分布を渦層外端参照点（c, d）での速度および壁面上でのNo-slip条件 に基づき図2.3.5のように直線勾配で近似すると，渦層内の渦度ω_izは次式により算出 できる．

h y

u y u x

ω_iz v ^z

2

2 1

1 u u

n

(2.3.110) ここで，u1 , u2は渦層外端参照点での速度ベクトル，nは壁面に対する法線方向単位ベ クトルである．

流れ場への渦度の流れ出しは，この薄い渦層内の渦度自身の移流と拡散を逐次計算 することによって，渦要素を流れ場へ流出させることで表現する．

壁面から高さhの渦層外面を通過する渦度の移流速度Vcは，連続の条件より次式か ら得ることができる．

2 2

1 1 2

1 2

1 1

f f

c h

V S u n u n (2.3.111)

ここで， _Sは速度参照点間長さ，n_{f 1} , n_f2は渦層両端面での渦層内を向く法線方向単 位ベクトルである．

図2.3.4. 壁面からの渦要素導入法（二次元）

thin vorticity layer ω_iz Vc： convective velocity

Vd： diffusion velocity

vortex blob

wall surface

c

a b

d

x, u y, v

h vortex sheet

Vc +Vd

ドキュメント内 渦法および粒子追跡法の内部流れ混相流解析への適用 (ページ 42-51)