数値解析手法

3.2.1 基礎方程式

本研究の基礎方程式は，非圧縮性粘性流体に対する連続方程式(3.1)，運動方程式として，

Navier-Stokes方程式(3.2)および(3.3)，VOF関数Fの移流方程式(3.4)によって構成される．

∂𝑢

d𝑥

+

^𝜕𝑤_𝜕𝑧

= 0

(3.1)

𝜕𝑢

𝜕𝑡

+ 𝑢

^𝜕𝑢_𝜕𝑥

+ 𝑤

^𝜕𝑢_𝜕𝑦

= 𝑋 −

¹_𝜌^𝜕𝑝_𝜕𝑥

+ 𝜈 (

^𝜕_𝜕𝑥^2𝑢₂

+

^𝜕_𝜕𝑧^2𝑢₂

)

(3.2)

𝜕𝑤

𝜕𝑡

+ 𝑢

^𝜕𝑤_𝜕𝑥

+ 𝑤

^𝜕𝑤_𝜕𝑧

= 𝑍 −

¹_𝜌^𝜕𝑝_𝜕𝑧

+ 𝜈 (

^𝜕_𝜕𝑥^2𝑤₂

+

^𝜕_𝜕𝑧^2𝑤₂

)

(3.3)

∂F

dt

+

^{𝜕(𝐹𝑢)}_𝜕𝑥

+

^{𝜕(𝐹𝑤)}_𝜕𝑧

= 0

(3.4)

ここで，u，wはそれぞれx，z方向の流速成分，X，Zはそれぞれx，z方向の外力，pは圧 力，νは動粘性係数，Fは流体の体積率である．

VOF法は元々Newton流体を対象とした解析手法であるが，本研究ではこれを，非Newton 流体を対象とした解析を行えるように拡張した．具体的には動粘性係数を，さまざまな非

Newton流体のひずみ速度とせん断応力の関係式から算出し，計算の過程で時々刻々変化す

るひずみ速度に応じて動粘性係数が時空間的に変化するよう改良を行った．

3.2.2 計算アルゴリズム

図2.1に計算アルゴリズムを示すように，解析領域，水深等の初期条件を入力する．そし て，流速が境界条件を満足するように Navier-Stokes 方程式を解くことにより，次の時間ス テップの流速を計算する．しかし，運動方程式から得られた流速は必ずしも連続方程式を 満足していないので，流速と圧力を調整しながら連続方程式を繰り返し計算する必要があ る。その後、連続方程式を満たした流速値を用いて，VOF関数Fの移流方程式を計算し，

自由表面形状を決定する．このように，上記の計算フローを各時間ステップで繰り返すこ とにより時系列解析を行うことができる．Newton 流体と非 Newton 流体の計算フローの違 いは，上記でも述べたが，動粘性係数の取り扱いである．Newton流体は定数値を与えるが，

非Newton流体では，時空間的に動粘性係数が変化するので，その都度求めなければならな

い．

41

図3.1 本計算手法のチャートフロー

Yes 反復 No

時 間 ステ ッ プ Newton流体

VOF関数F値を計算 Start

境界条件の設定 初期条件の入力

出力 t>t_end

End 連続式=0(発散)

自由表面の決定 N-S式の計算 動粘性係数(ν)の計算

ν=定数値 ν=変数値

ひずみ速度の計算

動粘性係数を算出 非Newton流体

Yes

No

42

上記で述べたように，Navier-Stokes式中に時間項が存在しない圧力は次時間ステップの値 を求めることができない．そこで，式(3.5)のように，連続式を満たすための発散Di,kを定義 して，Di,k = 0となるように修正する．

∂𝑢

∂𝑥

+

^∂𝑤_∂𝑧

=

^{𝑢𝑖+1/2,𝑘}_∆𝑥^{−𝑢𝑖−1/2,𝑘}

+

^{𝑤𝑖,𝑘+1/2}_∆𝑥^{−𝑤𝑖,𝑘−1/2}

= D

_ik

(

^Dij>0 圧力を減少

D_ij<0 圧力を増加

)

(3.5)

Di,kを解くにはNewton-Raphson法を用い，次式を得る．

δp_i,k^(m)= p_i,k^(m+1)− p_i,k^(m)= −D_i,k^(m)/ (∂D_i,k^(m)/ ∂p_i,k) (3.6)

(m)はm回目の繰り返し計算を表している．この方程式(3.6)に運動方程式(3.2)，(3.3)の𝑢^𝑛+1， 𝑤^𝑛+1を𝑢^(𝑚)，𝑤^(𝑚)に置き換え，(3.6)に当てはめて整理すると次式を得る．

δp_i,k^(m) = −ωD_i,k^(m)/ [∆t {_∆x¹

i(_∆x¹

i−1/2+_∆x¹

i+1/2) +_∆z¹

k(_∆z¹

k−1/2+_∆z¹

k+1/2)}] (3.7)

ここでωは反復計算の収束を早めるための緩和係数である．以上より，発散Dik = 0を満た

す流速𝑢^(𝑚+1)，𝑤^(𝑚+1)を以下の式で決定することができる．

p_i,k^(m+1)= p_i,k^(m)+ δp_i,k^(m) (3.8) u_i+1/2,k^(m+1) = u_i+1/2,k^(m) + ∆tδp_i,k^(m)/∆x_i+1/2 (3.9)

u_i−1/2,k^(m+1) = u_i−1/2,k^(m) − ∆tδp_i,k^(m)/∆x_i−1/2 (3.10)

w_i,k+1/2^(m+1) = w_i,k+1/2^(m) + ∆tδp_i,k^(m)/∆z_k+1/2 (3.11)

w_i,k−1/2^(m+1) = w_i,k−1/2^(m) − ∆tδp_i,k^(m)/∆z_k−1/2 (3.12)

本計算では，Di,kの収束判定基準はεp = 1.0×10^-3(1/s)とし，これを満足した流速と圧力が次の 時間ステップにおける値となる．

43

3.2.3 VOF法による自由表面の計算法

流体中に存在する物理量FをLagrange的に捉える方程式は式(3.13)となる．F = 0のとき

気体，F = 1のとき流体と仮定すると，式(3.13)のFは流体の体積率を表現しているようにみ

えるが，これはF = 0とF = 1という物理量が流速u，wで移動されることを表している．従

って，F = 0とF = 1に間に自由表面が存在することのみが分かり，Fは流体の体積率として

考えることができない．そこで，Fを体積率とみなすために，式(3.13)を上述した式(3.4)に 変形する．これにより，式(3.4)は流体領域だけでなく，全解析領域に用いることが可能とな る．つまり，F = 0のとき気体セル，F = 1のとき流体セル，0 < F < 1のとき表面セルとして 表現することができる．しかし，表面セルをVOF関数の値のみで判断すると，境界セルを 計算する際に不都合になる場合があるので，表面セルはVOF関数の値だけで判断せず，必 ず気体セルに隣接するという条件をつけなければならない．セルの分類については以下に 説明する．

∂F

dt

+ 𝑢

^𝜕𝐹_𝜕𝑥

+ 𝑤

^𝜕𝐹_𝜕𝑧

= 0

(3.13)

a) 解析セルの分類

次に図3.2に示すように，Staggered Meshを用いて計算領域をセル群に分割する．図に示す ように，各セルに流速u，wをセル境界上に，圧力PおよびVOF関数Fをセル中央に設定 する．図中のΔx(i)とΔy(k)はそれぞれx，z方向のメッシュの長さである．セル全体に流体が ある場合を流体セル(F：Fluid cell)，セル全体に気体がある場合を気体セル(E：Empty cell) および流体と気体が混在している場合は表面セル(S：Surface cell)と設定する．

図3.2 セルの分類とStaggered Meshの設定

44

b) 自由表面の向きの決定

先述したように，VOF関数であるF値によりセルの分類を行うことができるが，表面セ ルの自由表面の向きはF値のみでは決めることはできない．そこで，自由表面の向きを決 定するために，表面セルのどちらの方向に流体セルが存在しているのかみつける必要があ る．その方法は，自由表面の勾配を計算することにより，どの座標軸に垂直になるかを判 断することで，自由表面の向きを決定することができる．図3.3に示すように，セル(i,k)は

F値の0 < F < 1であり，気体セルと接しているため表面セルと判断できる．すると表面セル

の向きは①，②のどちらかになる．次に，セルにおける自由表面勾配を式(3.14)，式(3.17) で求め，比較することにより自由表面の向きを決定することができる．仮に，自由表面の 勾配が緩やかな方がより実際の自由表面に近いとすると, (dHZ/dx)＜(dHX/dz)のときの自由 表面形状はx軸に平行, 反対に(dHZ/dx)＞(dHX/dz)のときの自由表面形状はz軸に平行とな る．

(

^dHZ_dx

)

i,k

=

^{(dHZ/dx)i+1/2,k∆𝑥}_∆x^{𝑖−1/2+(dHZ/dx)i−1/2,k∆xi+1/2}

i−1/2+∆x_i+1/2 (3.14)

(

^dHZ_dx

)

i+1/2,k

=

^HZi+1,k+1_∆x ^−HZi,k

i+1/2 (3.15)

(

^dHZ_dx

)

i−1/2,k

=

^HZi,k_∆x^{−HZi−1,k−1}

i−1/2 (3.16)

(

^dHX_dz

)

i,k

=

^{(dHX/dz)i,k+1/2∆𝑧}_∆z^{𝑘−1/2+(dHX/dz)i,k−1/2∆zk+1/2}

i−1/2+∆z_i+1/2

(3.17)

① 自由表面がz軸に垂直な場合 ②自由表面がx軸に垂直な場合

図3.3 自由表面形状の決定

45 c) donor-acceptor法

表面セル，自由表面の向きが決定したら，donor-acceptor法により，気体，流体を移流さ せる必要がある．donor-acceptor法は，移流により移動するVOF関数Fがdonorセル(風上 側セル)とacceptor側セル(風下側セル)のF値によって決定される方法である．

例えば，図3.4に示すように，ある時間における移流量が赤枠の範囲だとすると，自由表 面が移流面に対して垂直な場合，移流面におけるVOF関数F値をdonor側のF値に一致さ せる．反対に，自由表面が移流面に対して並行の場合は，①，②のように，移流面におけ

るVOF関数F値をacceptorのF値と一致させる．しかし，決定される移流量によっては，

③，④のように，donorセルに十分な気体，流体がない場合がある．まず，③のように，donor セルに移流させるための気体が十分に存在しない場合は，そのときは不足分として流体を 移流させなければならない．一方，④のように，donorセルに移流させるための流体が十分 に存在しない場合は，不足分として気体を移流させる．上記のように，自由表面の形状や 気体・流体の存在状況を考慮すると式(3.18)となる．また式(3.18)中の式(3.19)，式(3.20)の

minはdonorセルが保有する流体以上の流体が移動するのを防ぎ，式(3.21)，式(3.22)のmax

はdonorセルが保有する気体以上の気体が移動するのを防いでいる．式中の添字については，

Dはdonorセル，ADはdonorセルを表わしている．

図3.4 VOF関数Fの移流方法の取り扱い

46

𝐹

_𝑖,𝑘^𝑛+1

= 𝐹

_𝑖,𝑘^𝑛

− [

^{𝑅𝑋𝑖+1/2,𝑘}_∆𝑥^{−𝑅𝑋𝑖−1/2,𝑘}

𝑖

+

^{𝑅𝑍𝑖+1/2,𝑘}_∆𝑧^{−𝑅𝑍𝑖−1/2,𝑘}

𝑘

]

(3.18)

R𝑋

_𝑖,𝑘

= 𝑠𝑖𝑔𝑛(𝑢

_𝑖,𝑘^𝑛+1

) ∙ 𝑚𝑖𝑛{𝐹

_𝐴𝐷

|𝑢

_𝑖,𝑘^𝑛+1

∆𝑡| + 𝐶𝐹𝑋, 𝐹

_𝐷

∆𝑥

_𝐷

}

(3.19)

R𝑍

_𝑖,𝑘

= 𝑠𝑖𝑔𝑛(𝑤

_𝑖,𝑘^𝑛+1

) ∙ 𝑚𝑖𝑛{𝐹

_𝐴𝐷

|𝑤

_𝑖,𝑘^𝑛+1

∆𝑡| + 𝐶𝐹𝑍, 𝐹

_𝐷

∆𝑧

_𝐷

}

(3.20)

CFX = 𝑚𝑎𝑥 {(1 − 𝐹

_𝐴𝐷

)|(𝑢

_𝑖,𝑘^𝑛+1

∆𝑡)| − (1 − 𝐹

_𝐷

)∆𝑥

_𝐷，

0}

(3.21)

CFZ = 𝑚𝑎𝑥 {(1 − 𝐹

_𝐴𝐷

)|(𝑤

_𝑖,𝑘^𝑛+1

∆𝑡)| − (1 − 𝐹

_𝐷

)∆𝑧

_𝐷，

0}

(3.22)

ドキュメント内 河川堤防の浸透性破壊を対象とした段階安全診断に関する研究 (ページ 42-48)