渦度場の時間変化の計算方法

(1) 渦度の移流による渦伸縮の計算方法

三次元流れでは，二次元とは異なり，渦度は粘性による拡散と共に，移流によって渦 度ベクトル方向の伸縮（渦度変化）を伴う．本研究では，球形の三次元渦モデル（渦Blob） によって，渦度移流による渦度ベクトル方向への伸縮を考慮する方法を用いる．

渦要素iを代表する渦度をω_i，位置ベクトルをr_i，渦要素の体積をdviとする．渦要素 は半径 _iの球する．この渦要素iによる任意位置rにおける渦度は次式によって表せる．

i i i

i i

dv f r r ω r

ω 1₃

(2.3.118)

ここに，fは球対称な渦度の分布関数であり，次式に示すようなWinckelmans＆Leonard^[56]

により提案された粘性核構造を用いる．

2 2 17

1 8

f 15 (2.3.119)

ここで，渦度の時間的変化は，次式に示す渦度輸送方程式に従う．

ω u

ω ω ₂

dt

d (2.3.120)

渦要素を代表する渦度の移流による時間的変化量は，Biot-Savart の式の勾配を取り，

渦要素及びポテンシャル流れそれぞれからの重ね合わせにより求められた速度勾配を上 式の右辺第一項に代入することにより与える．

移流による渦度の時間変化率は，上式の渦度輸送方程式の粘性項を除いた非粘性の渦

度輸送方程式に従う．

u ω ω

dt

d (2.3.121)

上式の速度uを，渦度による誘起速度uvorと，ポテンシャル流れによる速度 upotに分 けて，渦度の時間変化率を計算する．

u (r) = uvor + upot (2.3.122)

uvor：渦度による誘起速度

upot：ポテンシャル流れによる速度

初めに，流れ場に存在する渦要素からの移流による渦度の時間変化のみを考える．任 意の位置rにおける渦要素からの誘起速度u(r)は，渦要素の総数をNとすると，Biot-Savart の式より以下のように算出される．

i N

i V

i

vor g dv

R³ 4

1 ω R

r

u (2.3.123)

ここで，R，R， ，g は以下の通りである．

ri

r R R r ri

i

ri

r

2 5 2

3 5 0

1 5 2 2 4 f t dt 1

g (2.3.124)

この渦要素からの誘起速度を式（2.3.121）に代入して微分演算を行えば，渦要素から の移流による渦度の時間変化率を次式のように求めることができる．

i i i

i i

i

i i

i i i

i i V

dV dt

d

i

r r ω r r ω r

r r r

ω ω r

r r ω r

2 7 2 2

2 2

2 5 2 2 2 2

2 3 7

2 5 4

1

(2.3.125)

次に，境界面に分布させたポテンシャル値 およびポテンシャル流束値 nからの 移流による渦度の時間変化のみを考える．境界面をSとすると，解析領域の内部点rに おける速度ポテンシャルは，以下のように定式化される．

ds G n * n ds

* G

S

S r r

r (2.3.126)

ただし，r*は境界面を離散化したパネル中心の参照点（境界上の点）であり，n は境界 面での内向きの法線方向単位ベクトルである．また，R r r* とすれば，三次元流れ

におけるグリーン関数Gは以下のようになる．

G R 4

1 (2.3.127)

式（2.3.121）および式（2.3.122）から渦度の時間変化を計算するために，ポテンシャ

ル流れ場の速度および速度勾配を以下のように求める．任意の位置 r における および nからの誘起速度upot (r)は，成分表示すると以下のように算出できる．

ds , G n * ds , F

* u x

S j

S j

j j pot

r r

r r

j =1, 2, 3 (2.3.128)

i ij

j G, n

, F

ij j i

ij R R

G R 3 , ,

4

, 1 ₃

4 2

, , R G _i Rⁱ

R R,_i Rⁱ

i iR R R

i*

i

i r r

R i =1, 2, 3

ここで， _ijはクロネッカーのデルタ関数である．

また，速度 upotの速度勾配は，式（2.3.128）の勾配を取ることによって，成分表示す ると次式のように得られる．

ds n G

ds x F

u

S jl

S jl

l potj

,

* ,

* r

r r

( j, l ) =1, 2, 3 (2.3.129)

i ijl

jl G, n

, F

i jl j il l ij l j i

ijl R R R R R R

G R 5 , , , , , ,

4

, 3 ₄

以上から，渦要素を代表する渦度の時間的変化率をdω dt ωt で表すと，式（2.3.121） は，渦要素およびポテンシャルからの重ね合わせとして，以下のようになる．

t t

t ωvor ωpot

ω (2.3.130)

ωvor t ：渦要素による寄与

ωpot t ：ポテンシャル流れによる寄与

式（2.3.130）を用いれば，時刻t tにおける移流による渦度の変化は次の差分式で 近似的に計算できる．

t t t t

t ω ω

ω (2.3.131)

(2) 粘性拡散の計算方法

① 二次元モデル

本研究では，粘性効果の考慮法として渦要素の粘性核半径の拡散（増大）によって 考慮する方法を用いる．導入された渦要素について，粘性による渦度の減衰を考慮す るために粘性核半径 を各時間間隔ごとに変化させるCore-spreading法を用いる．

二次元流れにおいては，各渦要素は軸方向に無限長さをもつ渦糸と考えることがで きる．そこで二次元における渦度輸送方程式を円筒座標系に変形すると，以下のよう になる．

r r r t

ω ω

ω 1

2 2

(2.3.132)

上式より，時刻t において，無限長さを持つ循環量 ₀の渦糸が存在するとき，粘 性により渦運動が減衰する際の速度分布は，以下の式で与えられる．

t r t r

, r

u 1 exp 4

2

2

0 (2.3.133)

上式の速度分布において，時刻 tにおける速度の極大値を与える半径 を粘性核半 径とする．ここで，

t x r

4

2

(2.3.134) とおくと，次のように表される．

0

2 exp exp 1

1 1 2

4 exp exp

2 1

exp 2 1

exp 2 1

2 0

0 0

0 0

t x r x

t x r x

r x x x r

r r r u

(2.3.135)

上式より次のようになる．

x x exp 2

1 (2.3.136)

これを数値的に解くと x= -1.2564 を得る．したがって，次式が得られる．

t C t x

4 (2.3.137)

この の時間変化を考えると，次式となる．

2 2

C2

t C Dt

D (2.3.138)

ここに，C =2.2418である．

したがって，上式の微分方程式を，

t0

t のとき， _t0 t

t

t ₀ のとき， _{t0 t}

と置き，時間積分すると次のような関係式を得る．

C t

t t t t

0 0

0 2

2

(2.3.139) 上式から，各渦要素の持つ粘性核半径 を時間的に変化させることにより，粘性 拡散の効果を考慮する．

② 三次元モデル

本研究では，粘性効果の考慮法として，渦要素の粘性核半径を渦度移流による渦度 方向への伸縮と，粘性による渦度ベクトルと垂直な方向への拡散を適切に考慮する方 法を用いる．粘性拡散は，Core-spreading法による渦度変化を考える．

Blobモデルは，球対称な分布関数であるから時刻tにおいて渦度は，図2.3.7に示 すように半径 t の球体となっている．これが時刻t tでは，移流による渦度方向 への伸縮，及び渦度ベクトルと垂直な方向へ粘性拡散を受けて変形し，図2.3.7に示 すような2 _r，2 _lを軸とする楕円の回転体になると考えられる．本研究では，この 楕円の回転体の体積と等しい球の半径をt tのBlobモデルの粘性核半径とする．

3 1 2

l

t r

t (2.3.140)

または，近似的に次式となる．

3

2 _{r l}

t

t (2.3.141)

渦度の移流による渦度方向への伸縮 _l，及び伸縮によるコア半径の変化量 _sは，

時刻t及びt tの渦度に対し，それぞれ等価な渦管に置き換える．Helmholtzの定理 より体積一定であるので，次式により求まる．

t t t t

l ω

ω (2.3.142)

t t t

t t

s

2 1

ω

ω (2.3.143)

粘性拡散による変化 _vは，二次元の場合と同様にして次式により近似する．

t t C

v 2

2

（C=2.242） (2.3.144) そして，渦度と垂直方向への変化 _rは，渦の伸縮による変化 _sと粘性拡散による 変化 _vの重ね合わせで表現できるものとすれば，次式で与えられる．

v s

r t (2.3.145)

ここで，粘性拡散による渦度変化を除外し，渦度の移流のみ考慮した渦度輸送方程 式から得られた渦度ωt t を，粘性拡散を考慮して得られた粘性核半径 t t を 用いて修正する．時刻t及びt tにおける渦要素をそれぞれ等価な渦管要素と考え，

Kelvin の循環定理より，粘性拡散を考慮した渦度ω t t を次式により求め，これ

を時刻t tにおける渦度とする．

t t t

t t

t ^r ω

ω

2

(2.3.146)

s

r t (2.3.147)

図2.3.7. 渦要素の伸縮と拡散のモデル

(3) 渦要素の再配置

① 渦要素の急激な伸縮に対するモデル化

中西・亀本ら^[39]によって提案された三次元 Core-spreading 法は，渦度ベクトル方向 の伸縮を伴う移流による渦度変化，および粘性拡散による渦度変化を考慮する手法で ある．三次元流れへの適用性および拡張性が高く，これまで様々な解析に用いられて きた．ここで，三次元Core-spreading法では，毎ステップごとに渦要素を一つの球要素 に置き換える．そのため，渦要素の急激な変形を伴うような速度変動が大きな領域で は，渦構造の変形を解像度良く再現することができない．また，乱流場のように非定 常性が強く，局所的に速度変動の大きさが異なる流れ場においては，空間解像度を維 持するために要求される渦要素スケールが，空間および時間において同一ではない．

このため，流れ場全体で空間解像度を維持するためには，非常に微細な渦要素を多数 導入することが必要となり，計算負荷が増大する．

一方，Biot-Savart 法を用いた渦法においては，流れ場に格子を設けないため，渦要

素のスケールを局所的に変化させることが可能である．そこで，福田ら^[48]は，渦要素 の変形量を見積もることで，変形が大きな渦要素をより小さな渦要素に分割して再配 置するラグランジュ型渦要素再配置モデルを提案した．この手法において，再配置さ れる渦要素のスケール・要素数は，渦要素の変形量により決定される．このような局 所渦要素再配置モデルを用いることで，速度変動が大きな領域においては，小スケー ルの渦要素を用いて渦度分布を離散化することになり，空間解像度を高めることが可 能となる．また，渦要素の変形に応じて渦要素スケールが選択されるため，速度変動 の大きな領域にのみ小スケールの渦要素が導入され，効率の良い解析が可能となる．

本研究では，福田らが提案したラグランジュ型渦要素再配置モデルを，混相流内部 流れの解析へ適用する．ここでは，福田らが提案した渦要素再配置モデルの詳細につ いて記す．

② 渦要素の再配置

渦要素が渦度ベクトル方向へ急激な引き伸ばしを受ける場合には，その渦度は増加 する．従来の三次元Core-spreading法では変形した要素を一つの球要素に再び置き換え るため，渦度ベクトル方向の伸張効果が大きな領域では渦度ベクトル方向の解像度が 低下することがあり，急激な渦度の伸張に対する再配置モデルが必要となる．福田ら が提案したラグランジュ型渦要素再配置モデルでは，要素の持つ軸方向長さおよび断 面半径を用いて，渦度ベクトル方向の伸張効果に基づき要素を分割・再配置すること で空間分解能を維持する．モデルの詳細を以下に記す．

時 刻 t に お け る 渦 要 素 の 長 さl_t お よ び 断 面 半 径 _t の 時 間 変 化 率 を 三 次 元

Core-spreading 法を用いて見積もることで， t時間後の渦要素の長さl_{t t}および断面

半径 _{t t}を求めることができる．図2.3.8に示すように，変形した渦要素をn個の渦要 素に再配置する．再配置される渦要素数nは，以下の式で決定する．

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