Development of a Simulation Method for Fluid-Structure Interaction Problems with Free Surface based on Finite Cover Method

(1)

修士論文要旨（2010年度）

有限被覆法に基づく自由表面を有する流体―構造連成解析手法の構築研究

Development of a Simulation Method for Fluid-Structure Interaction Problems with Free Surface based on Finite Cover Method

土木工学専攻 32号中村正人 Masato NAKAMURA

1. はじめに

実現象において，自由表面を有する流体と構造物が相互作用する問題は数多く存在する．土木工学における問題例としては，土石流や河川，海岸における浮体構造物の挙動などが挙げられる．これらの問題の評価方法として，数値シミュレーション手法が近年有効に用いられているが，それらの手法は移動メッシュを用いる手法と固定メッシュを用いる手法に大別される．移動メッシュを用いる手法は界面を直接的に扱うため解析精度の点で有効であるが，跳水や砕波を含むような複雑な解析を行う場合，解析メッシュに破綻が生じ解析が困難となる場合が多く，ロバスト性に問題がある．

本研究では，ロバスト性に優れる固定メッシュに基づく手法に着目して手法の開発を行う．しかし，固定メッシュを用いる手法では，流体と構造の界面での幾何学的境界条件及び力学的境界条件を正確に考慮することが一般に困難である．そこで，この問題を解決する手法として，計算固体力学の分野において提案された有限被覆法

1)（Finite Cover Method，以下FCM）を流体―構造連

成解析に適用することを試みた．なお，流体と構造の界面における境界条件処理にはペナルティ法を用いた．また，流体は非圧縮性粘性流体を仮定し，気体と液体の境界面の表現にはVOF関数を用い，その境界位置はVOF 関数に関する移流方程式をCIVA法²⁾により解くことで決定する．また、流体と構造の境界面の表現にはLevel Set関数を用いる．なお，連成解析手法には弱連成手法を用いる．

2. 数値解析手法

2.1 有限被覆法(FCM)

FCMは，図−1に示すように近似関数が定義される数学領域Ω^M と，支配方程式が満たされるべき物理領域Ω^P を独立して定義するという点で，有限要素法（以下，FEM）と大きく異なる．しかし，解析対象を要素で部分分割し，各要素間の未知量を節点値により補間近似するという点においてFEMと一致するため，FEM を一般化した手法とみなすことができる．以上の特性から，FCMでは，解析対象とは独立に節点を単位とする PU条件（Partition of Unity）を満たす近似関数（重み関数）を定義する．本論文では，任意形状を有する解析領域の適用性に優れている三角形要素に着目し，その中でも取り扱いの簡単な１次の三角形有限要素の形状関数

Gas

Liquid Solid

mathematical domain physical domain

図– 1 数学領域 Ω^M と物理領域Ω^P

を重み関数とする．FCMは従来のFEMとは異なり，

要素内に流体と構造の境界が存在することを許容するため，その境界位置を正確に考慮した流れ場を求めることが可能である．

2.2 流体に関する定式化 2.2.1 支配方程式

Euler記述された非圧縮性粘性流体の運動方程式及び，連続式はそれぞれ以下の(1)，(2)で表される．

ρ³∂u

∂t +u· ∇u−f´

− ∇ ·σ(u, p) = 0 (1)

∇ ·u= 0 (2)

ここで，ρは密度，uは流速ベクトル，f は物体力ベクトル，σは応力テンソルを表す．ここで，対象とする流体はNewton流体であるため，応力テンソルσは以下の式(3)で表される．

σ=−pI+µ£

∇u+ (∇u)^T¤

(3)

ここで，pは圧力，µは粘性係数である．Dirichlet型，

Neumann型境界条件は，それぞれ次式で与えられる．

u=g on Γg (4)

n·σ=h on Γ_h (5)

ここで，g，hはそれぞれ流速，トラクションの既知量を示し，nは外向き法線ベクトルを示す．また，自由表面位置は，次式(6)の移流方程式を解くことで決定される．

∂φ

∂t +u· ∇φ= 0 (6) ここで，φはVOF関数（界面関数）を表し，気体であれば0.0，液体であれば1.0，自由表面上であれば0.5となる．なお，各要素における気体，液体の密度と粘性係数は，VOF関数を用いて次式のように求められる．

ρ=ρliqφ+ρgas(1−φ) (7)

(2)

µ=µliqφ+µgas(1−φ) (8)

ここで，ρ_liq，ρ_gas，µ_liq，µ_gasはそれぞれ液体，気体の密度及び粘性係数である．

2.2.2 SUPG/PSPG法に基づく安定化手法の適用支配方程式(1)，(2)に対して，SUPG/PSPG法に基づくFCMを適用すると，次式の弱形式(9)が得られる．

ρ Z

Ω^P

µ∂u

∂t +u· ∇u−f

¶

·u^∗dΩ

− Z

Ω^P

p∇ ·u^∗dΩ +µ Z

Ω^P

³∇u:∇u^∗+∇u: (∇u^∗)^T´ dΩ

+ Z

Ω^P

q^∗∇ ·udΩ +

nel

X

P=1

Z

Ω^P{τ_supgu· ∇u^∗+τ_pspg1 ρ∇q^∗}

·{ρ µ∂u

∂t +u· ∇u−f

¶

− ∇ ·σ}dΩ +

nel

X

P=1

Z

Ω^P

τ_cont∇ ·u^∗ρ∇ ·udΩ

| {z }

衝撃捕捉項

+ Z

Γ_g

¯

p(u−u)dΓ¯

| {z }

ペナルティ項

= 0 (9)

ここで，u^∗，q^∗はそれぞれ運動方程式と連続式に対する重み関数を表す．また，τ_supg，τ_pspg，τ_contは安定化パラメータを表す．式(9)中の衝撃捕捉項は自由表面付近の数値振動を回避するためのもので，ペナルティ項は構造壁面上での境界条件を満足させるために付加したものである．積分領域が物理領域Ω^PであることがFEMと大きく異なる点である．式(9)に対して，P1/P1（流速・

圧力一次）要素を用いて補間を行い，時間方向の離散化にはCrank-Nicolson法を用いる．また，移流速度は2 次精度Adams-Bashforth法により陽的に近似する．連立一次方程式の解法にはElement by Elementに基づく GPBi-CG法を用いる．

一方，気体，液体の界面関数であるVOF関数に関する移流方程式(6)の解法にはCIVA法を用いる．

2.2.3 FCMにおける領域積分

FCMでは，従来のFEMとは異なり，流体と構造の境界とは独立に近似領域を定義できるため，空間固定のメッシュを配置することができる．その際，物理領域で完全に満たされた要素と要素内に物理領域が部分的に存在する要素が生成される．FEMの形状関数を重み関数に採用していれば，前者の要素は通常の有限要素と同等であり，後者の要素は領域積分を工夫することで弱形式 (9)を満足させることができる．ここでは本論文で用いる領域積分の積分方法について説明する．FCMでは要素内に物理領域が部分的に存在する要素では，領域積分を行う際，面積座標の積分公式を直接利用できないため，

図−2の(a)に示すような面積座標系を用いて数値積分

（ガウス積分）を行う．なお，要素内の積分領域が四角形の場合、図のように細分化された領域ごとに積分を行

a b Gas

Liquid Solid

: Fluid-Structure Interface : Subdivision

: Degree of freedom (Fluid)

: Gaussian integration point : Intersection points Line integral

based on natural coordinate

Domain integral based on area coordinate

a

b c

(b)

(a)

図– 2 領域積分と線積分

う．流体領域Ω^Pの領域積分は次式(10)で表される．

Z

Ω^P

f(x, y)dΩ =

n_g

X

j=1

f(L_1j, L_2j, L_3j)W_j^gA_f (10)

ここで，f は被積分関数，ngは積分点の数，Lは面積座標，W^g は重み，Af は流体領域の面積を表す．本論文では，図のように3点での積分を行うため，ng = 3とする．

2.2.4 ペナルティ法による境界条件処理

図−2の(b)のように流体と構造の境界上に自由度を持つ節点が一般には存在しないため，節点値の処理による境界条件処理を行うことができない．そのため，

本論文ではペナルティ法による境界条件処理を行う．

なお，ペナルティ項の積分においては自然座標系に基づくガウス積分を適用する．ガウス積分は以下に示す Gauss-Legendreの公式(11)を用いる．

Z b a

f(x)dΓ∼= b−a 2

n_g

X

i=1

c_if µa+b

2 +b−a 2 ξ_i

¶ (11)

ここで，f は被積分関数，n_gは積分点の数，c_iは重み係数，ξは積分点を表し，領域積分と同様n_g= 3とする．

2.3 構造に関する定式化

構造体を剛体と仮定すると，運動方程式は次式で表される．

md²x

dt² =F−mg (12) ここで，mは物体の質量，xは物体の変位，tは時間，F は物体が受ける流体力，gは重力等の外力を表す．この

(3)

o

x y

Traction-free

No-slip

図– 3 解析モデル

FCM (Mesh A)

FCM (Mesh B)

FCM (Mesh C) Fluid

: Wall : Interface Penalty Method

図– 4 FCMメッシュと Mesh Aの境界近傍拡大図

図– 5 W = 3における流速分布

微分方程式の解法にはNewmark-β法を用いることで，

時刻ステップn+ 1での加速度，速度，変位を求めて，

時間進行していく．

3. 数値解析例

3.1 ポアズイユ流れ問題

FCMにおける境界条件処理の妥当性について検討するために，固定メッシュ（メッシュと構造物の幾何形状が一致しない）を用いたポアズイユ流れ問題を行う．解析モデルを図−3に示す．レイノルズ数はRe = 10，微小時間増分量は∆t = 10⁻³を用いる．本数値解析例では，厳密解，FEMによる解析結果との比較を行い，

本手法の妥当性を示す．解析に用いたFCMメッシュを図−4に示す．Mesh A，B，Cにおける領域内部の線は流体と構造の境界を表す．それぞれ要素内における境界の位置は異なるが流体領域は同じである．図−5に W = 3における流速分布を示す．図より，本手法の解析結果は厳密解，FEMの結果とほぼ一致していることから，本手法の妥当性が確認できる．

FEM (reference) FCM

: Cylinder : Interface

図– 7 円柱近傍のメッシュ図（左： FEM，右：FCM）

3.2 円柱周り流れ問題

FCMにおける構造物に作用する流体力の評価の妥当性について検討するために円柱周り流れ問題を行う．解析モデルを図−6に示す．レイノルズ数はRe= 100，微小時間増分量は∆t= 10⁻³を用いる．境界条件として，円柱上でNo-slip条件，側壁でSlip条件，下流側流出境界ではトラクション・フリー（traction-free）条件とした．本数値解析例では，FEMの解析結果との比較を行い，本手法の妥当性を示す．解析に用いたメッシュの円柱近傍拡大図を図−7に示す．左側は参照解である FEMの解析に用いるメッシュ，右側はFCMの解析に用いるメッシュである．FCMの解析に用いるメッシュは要素内に流体と構造の界面が存在していることがわかる．図−8に抗力係数の時刻歴と平均値を示す．図より，カルマン渦の発生時刻が若干異なるが，本手法の結果（抗力係数1.4022）は参照解であるFEMの結果（抗力係数1.3892）とほぼ一致している．このことから，界面位置とメッシュの要素辺の位置が一致しない場合においても構造物に作用する流体力を定量的に評価できているといえる．

3.3 楔型構造物の落下問題

最後に応用例として，楔型構造物の水中への落下問題

3)を行う．図−9に解析モデル図を示す．構造物境界上

でNo-Slip条件，四角形領域の周りの4辺でSlip条件とした．また，青色の領域は液体部分を表す．ここで，

構造物は実験においてスライドレールに沿って落下するため，水平方向，回転方向に拘束されている．また，

構造物の重力加速度は，スライドレールとの摩擦により

(4)

Average

図– 8 抗力係数の時刻歴と平均値

8.2015 m/s²となっている．構造物の単位長さ当たりの質量は22.537 kg/m，先端角度β=π/4である．また，

構造物先端が水面に接した状態を0 sとし，その際の構造物の初期速度は1.5797 m/sとする．

図−10, 11に解析結果を示す．図−10に構造物の着水後の自由表面形状と圧力分布を示す．図より構造物と自由表面を有する流体との相互作用をロバストに解析できていることがわかる．図−11に構造物の加速度の時刻歴を示す．若干の振動はみられるものの実験値と概ね良い一致を示していることが確認できる．以上の結果より，本手法の妥当性を示すことができる．

1.65 m

1.1 m 035 m

0.2 m β

x y

4. おわりに

本論文では，複雑な自由表面を有する流体―構造連成解析を行うことを目的として，FCMに基づく流体―構造連成解析手法の構築を行い，以下の結論を得た．

• ポアズイユ流れ問題において，要素内に流体と構造の境界が存在する場合においても，厳密解，参照解であるFEMによる結果とほぼ同等の流速分布が得られ，本手法の妥当性が確認できた．

0.02 [s]

0.04 [s]

4000

0 [Pa]

0.02 [s]

0.04 [s]

図– 10 各時刻における自由表面形状と圧力分布

図– 11 構造物加速度の時刻歴

• 円柱周り流れ問題において，要素内に流体と構造の境界が存在する場合においても，参照解である FEMによる結果とほぼ同等の抗力係数が得られ，

構造物に作用する流体力を定量的に評価できていることが確認できた．

• 楔型構造物の落下問題において，自由表面を有する流体と構造物との相互作用現象をロバストに解析でき，実験値とも概ね良い一致を示した．

今後の課題として，構造物の変形等の考慮を行う予定である．

参考文献

1) Kenjiro Terada，Mitsuteru Asai，Michihiro Yamag- ishi：Finite cover method for linear and non-linear analyses of heterogeneous solids：Int. J. Numer. Meth.

Engng, Vol.58, pp.1321-1346，2003.

2) 田中伸厚：数値流体力学のための高精度メッシュフリー手法の開発，日本機械学会論文集(B編)，64巻620号，

1998.

3) G.X.Wu, H.Sun, Y.S.He：Numerical simulation and experimental study of water entry of a wedge in free fall motion：Journal of Fluids and Structures，Vol.19, pp.227-289，2004.