牛島 省

水工学論文集，第52巻，2008年2月

MICS と有限要素法による自由水面流と 弾性体の連成運動に対する 3 ^{次元数値計算}

3D NUMERICAL PREDICTION FOR INTERACTION BETWEEN FREE-SURFACE FLOWS AND ELASTIC BODIES WITH MICS AND FINITE ELEMENT METHOD

牛島 省

・黒田 望

・禰津 家久

Satoru USHIJIMA, Nozomu KURODA and Iehisa NEZU

1正会員 工博 京都大学大学院准教授 社会基盤工学専攻(〒615-8540京都市西京区京都大学桂Cクラスタ)

2学生員 京都大学工学部 地球工学科(〒615-8540京都市西京区京都大学桂Cクラスタ)

3フェロー会員 工博 京都大学大学院教授 社会基盤工学専攻(〒615-8540京都市西京区京都大学桂Cクラスタ) A computational method has been proposed to predict the interactions between free-surface flows and elastic bodies included in the flows. A solid model, whose deformations due to fluid forces are solved with a finite element method, is introduced into the MICS, a computational method for incompressible multiphase fields. A solid object included in the flows is divided into multiple tetrahedron elements, through which fluid-solid interactions are taken into account using a tetrahedron sub-cell method. The developed computational method was applied to the experimental results, which were obtained in an accelerated tank with elastic plates. As a result, the validity of the prediction method was confirmed.

KeyW ords : fluid-solid interaction, free-surface flow, elastic body, FEM, MICS

１．はじめに

流れの影響を受けて変形する物体の挙動を把握する ことは，各種の工学的な問題において重要な課題となっ ている．流れに対する物体の変形応答は，プラント等 における流体関連振動や，波浪に対する弾性浮体構造 物の応答，地震時に発生するスロッシングと薄肉構造 物の変形など，多くの問題に関連している．水工学分 野では，変形を伴う植生の流体抵抗の変化や，洪水時 における河畔林の倒伏の問題などに関係している．

流体と構造物の連成問題に関する数値的な検討例と しては，1自由度あるいは2自由度の一様流中の弾性 円柱の振動を扱う数値計算が行われている¹⁾．一方，

自由水面流れが関係する多自由度の物体変形計算を行 う連成解析の例は比較的少数である²⁾．

本報では，多相場の数値解法(MICS)に，物体の変 形を有限要素法で求める固体モデルを導入する．自由 水面流中に存在する物体は弾性体であるとして，これ を四面体1次要素で表現する．そして，有限要素法に より離散化された弾性体の動的問題の基礎式を，MICS による多相場の計算と連成させて解く手法を示す．複 数の弾性板の一端が底面に固定された水槽に水を満た し，水槽を加振することにより，水面変動と弾性板の 変形を発生させる実験を行い，これを対象とする数値 計算を行って解法の適用性を検討する．

２．数値解析手法

自由水面流れの中に存在する弾性体が，流体力によっ て有意な変形を示す問題を扱う数値解法では，流れと 弾性体の動的応答を適切に扱う必要がある．本研究で は，このような問題に対する解法を構築するため，多 相場の解法であるMICS³⁾に，変形や動的特性を有限 要素法によって計算する固体モデルを導入する．

(1) 3次元自由水面流れの計算法

MICSにおける多相場の解法の概要を以下に示す．

基礎式は，以下のEuler表記による質量保存則，非圧 縮条件，保存形表示された運動方程式の3式である．

∂ρ

∂t + ∂

∂xj(ρuj) = 0 (1)

∂uj

∂xj = 0 (2)

∂ui

∂t + ∂

∂xj(uiuj) =fi−1 ρ

∂p

∂xi

+1 ρ

∂

∂xj

∂

∂xj(µui) + ∂

∂xi(µuj)

(3)

tとxiは時間と3次元直交座標系の座標成分である．

ρ，µ，pは順に計算セル内の体積平均操作によって求 められる密度，粘性率，圧力である．また，uiはセル 内の質量平均により算出される流速成分である．fiは 外力の加速度成分を表す．

水工学論文集,第52巻,2008年2月

計算方法の概要は以下のとおりである．最初に，四 面体要素のサブセル法⁴⁾により，計算セルに含まれ る物体体積を算出し，体積平均された物性値等を求め る．次に，コロケート格子を用いる非圧縮性流体計算 法に従い，まずセル中心で流速の推定値を求め(予測 段階)，これをセル境界に空間内挿して圧力勾配を考 慮し，C-HSMAC法による圧力計算を行う．予測段階 の計算では，陰的解法であるC-ISMAC法⁵⁾ を使う．

自由水面形状は，式(1)を数値拡散を抑制する保存形 スキームで解いて求める．

(2) 有限要素法による物体計算

本報では，流体中に存在する物体は，弾性体である とする．物体の動的挙動を数値的に扱う場合には，流 体力が外力として作用したときの物体の変形量や各点 の加速度および速度を正確に計算することが重要とな る．既報⁶⁾では，物体を図–1に示すような四面体要 素により表現する固体モデルを提案した．この固体モ デル(T型質点バネモデル)では，四面体要素の各節点 (node)に質点が配置され，それらがダンパ付きのバネ で接続される．この固体モデルは数値的に比較的安定 で取り扱いやすく，水面変動に対する弾性円柱群の動 的な挙動が再現されることが示されている⁶⁾．

node

element

図–1 四面体要素と節点

一方，この質点バネモデルを利用する場合には，バ ネ定数や減衰係数などのパラメータを，対象物を用い た変形実験と試行計算との比較を通じて定める手続き が必要となる．そこで，本報では物体を弾性体と仮定 し，ヤング率などの一般的なパラメータを利用する有 限要素法により，物体の変形や動的挙動を扱う解法を 新たに作成した．以下では，この解法をT型FEMモ デルと表記する．

(3) 物体の動的応答の基礎式

T型FEMモデルでは，図–1のように物体を四面体 要素に分割し，その各節点上に変数を定義する．本報 では，形状関数が座標の1次関数で表される，簡単な 四面体1次要素を利用することとした．有限要素法に よる連続的な弾性体の動的計算では，質量マトリック スM と各節点の加速度成分から構成される加速度ベ クトルとの積により慣性力が定義され，これに弾性力，

減衰力および加振力(外力)を考慮した方程式が動的応 答の基礎式となる．

各節点の3次元変位を成分とするベクトルをdとす れば，物体の動的挙動に関する支配方程式は，次式で 与えられる．

Md¨+Cd˙+Kd=f_e (4)

ここで，上付のドットは時間微分(2つのドットは2階 微分)を表し，Cは減衰マトリックス，Kは剛性マト リックス，f_eは流体力などの外力ベクトルである．質 量マトリックスM の要素は，一般に物体の密度ρbに 形状関数を乗じて，これを四面体要素内で積分して得 られた行列から構成される．本報では，集中質量を対 角要素とする集中質量マトリックスを用いる．

式(4)の減衰マトリックスCは，質量マトリックス と同様に，減衰係数と形状関数の積を要素内で積分し て得られた行列から構成される．本報では，集中質量 マトリックスと同様に，対角行列として表される減衰 マトリックスを利用する．

式(4)の剛性マトリックスKは，各要素で得られる 要素剛性マトリックスKeから構成される．四面体1 次要素を用いる場合には，要素剛性マトリックスKe

は，歪みと変位の関係，応力と歪みの関係，そしてエ ネルギー原理により導かれ，次のように表される．

Ke=V B^TDB (5)

ここに，V は四面体要素の体積であり，B^T は行列B の転置行列を表す．行列Bは次の歪みベクトルと変 位ベクトルdの関係を表す6×12行列である．

=Bd (6)

行列B の要素は，形状関数を座標成分で偏微分した 値から定められる．また，式(5)のDは，次の応力ベ クトルσと歪みベクトルの関係を表す6×6行列で ある．

σ=D (7)

行列Dの要素は，ヤング率およびせん断弾性係数，ポ アソン比から定められる．

(4) 動的応答の基礎式の解法

動的応答の基礎式である式(4)を時間積分すること により，節点の速度と変位の3次元成分が得られる．

時間方向には，差分法を用いて離散化を行う．本報で は，次式のように，オイラー陽解法を用いてn+ 1ス テップの節点の速度ベクトルd˙を求める．

d˙ⁿ⁺¹=d˙ⁿ+M⁻¹

fⁿ_e −Cd˙ⁿ−Kdⁿ

∆t (8)

剛性マトリックスKは，式(5)の要素剛性マトリックス Keを重ね合わせることにより得られる．しかし，本報 では，この計算処理を省くため，element-by-element 法⁷⁾における行列ベクトル積の算出法と同様に，Kを 構成せずにベクトルKdを直接計算することとした．

すなわち，Kdを次のような演算処理により求める．

v(:, :) = 0.0d0

do ne = 1, nttr ! nttr^{は全要素数} (A) do nt1 = 1, 4 ! nt1^{は局所節点番号} (B)

nt1に対するグローバル節点番号nd1^を取得 do nt2 = 1, 4 ! nt2^{は局所節点番号} (C)

nt2に対するグローバル節点番号nd2^を取得 do j = 1, 3 ! (D)

do i = 1, 3 ! (E)

v(i, nd1) = v(i, nd1) &

+ ke(i, j, nt1, nt2, ne) &

* d(j, nd2) ! (F) enddo

enddo enddo enddo enddo

上記の(A)は全要素，また(B)と(C)は四面体要素 の節点に関するループ演算である．(D)と(E)は3次 元成分に関するループであり，(F)の演算で要素剛性 マトリックスKeの要素と節点変位ベクトルdの成分 との積を配列vに加算する．上記の演算では全体の剛 性マトリックスKを構成する必要はなく，配列vに 行列ベクトル積Kdの結果が格納される．

節点の変位は，次式のように，式(8)に完全陰解法 を用いて求められる．

dⁿ⁺¹ =dⁿ+d˙ⁿ⁺¹∆t (9) (5) 流体と物体の動的挙動の相互作用

物体に作用する流体力は，多相場の計算により得ら れた圧力勾配項と粘性拡散項から計算される．MICS との連成計算では，体積積分を通じて力と運動量の交 換が考慮される点に特徴がある．

T型FEMモデルでは，四面体要素の節点上の流体 力を求める必要があるが，最初に物体を構成する四面 体領域に作用する流体力を以下のようにして計算する．

図–2に，流体計算セルと物体を構成する四面体要素

の関係を示す．実際には，3次元であるが，ここでは セルを直方体，四面体要素を三角形で概略的に表示し ている．図–2に示すように，流体計算セルC内の多 相流体が，セル内に含まれる物体kの四面体要素Tkm

あるいはその一部分の体積∆TCkmに及ぼす流体力を FCkmとし，そのxi方向成分をF_Ckmⁱ と表す．

T

object-k C

∆ T

F

F

F

F

図–2 物体に作用する流体力の評価方法

F_Ckmⁱ は，四面体サブセル法⁴⁾ により求められた

∆TCkmと物体kの密度ρbkを用いて，次式から求め られる．

F_Ckmⁱ =ρbk∆TCkm

−1 ρ

∂p

∂xi

+1 ρ

∂

∂xj

∂

∂xj(µui) + ∂

∂xi(µuj)

(10)

次に，T型FEMモデルでは，上記のようにして得 られたFCkm を各節点上の流体力に変換する．すなわ ち，図–2に示すように，式(10)から得られるFCkm

に対して，セル中心からの距離の逆数に相当する重み 付けを行い，これを四面体節点上の流体力FCkmj と する(j = 1,· · · ,4)．ある節点に対して，その節点を 含む全ての四面体要素と，その要素を含む流体計算セ ルに対してFCkmjの総和を求め，その結果を式(4)右 辺の外力f_eとする．

一方，物体の動的応答計算の結果は，多相流場に反 映される．図–3にその概要を示す．ある流体計算セル Cに対して，それに含まれる物体kの四面体要素Tkm

を選択する．式(8)より得られた要素Tkmの節点の速 度ベクトルd˙をvkmjとするとき(j = 1,· · · ,4)，こ れらの和を1/4倍した結果を四面体要素の速度ベクト ルvkmと近似する．着目した流体計算セルに含まれる 全ての四面体要素に対してこの処理を行い，次式より セル内の質量平均流速uを定める．

u= 1 mC

mfuf+

k

m

ρbk∆TCkmvkm (11)

T km

C

∆ T Ckm

v _km v _km3 v _km1

object-k v _km2

図–3 物体の動的挙動を多相場に考慮する方法

ここで，mCとmfは，それぞれ着目する流体計算セ ル内の全質量および気相と液相の質量である．また，

ufは気相と液相の混合体の流速ベクトルである．な お，体積∆TCkmは上述のように，四面体サブセル法 により求められる．

３．水理実験と解法の検証

(1) 弾性板を用いる水理実験

上記で述べた数値解法の適用性を確認するため，弾 性板を内部に備えた水槽に水平加速度を与えて，水面 変動と流れを発生させ，弾性板の変形を調べる実験を 行った．この実験結果を計算結果と比較する．

実験に用いた水槽を図–4に示す．水槽は長さLが 190 mm，奥行きbが60 mmの矩形容器で，内部に4 枚の弾性板が設置されている．弾性板は，下面が水槽 底面に接着され，片持ち梁の状態となっている．これ らの弾性板は，イノアックコーポレーション製のポロ ンスポンジで，高さhbは100 mm，幅b1は30 mm であり，板厚は10 mmである．弾性板の比重は0.255 であり，荷重を加えて変形させたところ，ほぼ弾性体 と見なせる特性があることが確認された．図–4の弾性 板の配置は，d1とd2がそれぞれ45および35 mmで あり，y方向には水槽中央に固定されていて，b2は15 mmである．

実験では，水槽内に初期水深がh0となるように水 を入れ，静水状態とした後，図–4のx方向に電動ス ライダを用いて水平加速度を加えた．約0.2秒間，正 の一定の加速度axを加えた後，同じ時間にわたり加 速度−a_xを加えて水槽を停止させた．この加振により

x z x

y ^b

b

1

b

b d

d

d

d

d

h

L

h

図–4 実験水槽と座標系(上は側面図，下は平面図)

水面変動と流れが生じ，それに伴って弾性板が変形す る．この状況を側面からビデオで撮影し，ビデオ画像 を解析して弾性板の変形量などを求めた．

実験では，axは1.0 m/s²とした．また，初期水深 h0は，75，100，125 mmの3条件とした．以下では，

これらを順にケースh75，h100，h125と表す．

(2) 計算結果との比較

計算では，水と上部の空気部分および弾性板を含む 3次元直方体領域に対して，各方向5 mmの大きさの 流体計算セルを設定した．水と空気の動粘性係数は，

それぞれ1.0×10⁻⁶および1.0×10⁻⁵m²/sとし，密 度はそれぞれ1.0×10³および1.0 kg/m³ とした．弾 性板の密度は実験値と同様にした．

弾性板は，節点数358，要素数1,032の四面体要素 により表現した．後述するように，有限要素法の計算 負荷は大きく，計算時間の関係から，節点数は少なく 設定されている．四面体1次要素を用いる物体変形の 計算では，より多くの節点数が必要であると考えられ るが，この点は今後の検討課題である．計算では，上 記の節点数で弾性板の変形をほぼ適切に表現できるよ うに，弾性板のヤング率を1.5×10⁵Paとし，単位体 積あたりの減衰係数は2.0 ×10⁴ N s/m/m³とした．

図–5と図–6にそれぞれケースh75およびh100の計 算結果を示す．ケースh75では，弾性板の上部は常時 水面上に露出しているが，ケースh100では水面は弾 性板上端の上下に変動する．なお，ケースh125では 弾性板は常時水面下にある状態となった．水面変動に 伴い，容器内には流れが発生する．図–5と図–6では，

弾性板周辺の渦度の水平分布を合わせて示している．

t= 0.2 (s)

t= 0.4 (s)

t= 0.6 (s)

t= 0.8 (s)

図–5 自由水面流れと弾性体の変形の計算結果 (ケースh75，等高線は渦度を表す)

t= 0.2 (s)

t= 0.4 (s)

t= 0.6 (s)

t= 0.8 (s)

図–6 自由水面流れと弾性体の変形の計算結果 (ケースh100，等高線は渦度を表す)

図–7に容器両端(x= 0, L)における水位変動量η0

およびηLと，図–4の左から2番目の弾性板先端のx 方向の変位量dtの時系列を示す．実験で得られた最大 水位変動量ηmは，計算結果とほぼ一致している．ま た，弾性板先端は水位変動と概ね1/4波長ずれた変化 を示している．これより，弾性板は流速が最大となる ときの流体力により変形していると考えられる．

-0.03 -0.02 -0.01 0.00 0.01 0.02 0.03

0.0 0.2 0.4 0.6 0.8 1.0

exp ( ) cal ( ) cal ( ) cal (dt)

η,

t (s) η

η

dt (m)

0

L

η_L η η₀^m

図–7 水位変動と弾性板変位の時系列(ケースh100)

図–8に弾性板先端の変位dtの時系列を示す．実験お よび計算結果では，両端2枚および中央2枚の弾性板 の変位は，概ね同様の変動を示した．そのため，図–8 の実験結果では，それらの平均をout，inと表してい る．前者より，後者の変位が一般に大きくなった．

実験および計算では，図–5から図–7の結果に示さ れるように，水槽のx方向の長さLを半波長とする1 次モードに近い水面変動が見られた．この水面変動の 固有周期は，水深75および100 mm に対して，それ ぞれ約0.53および0.51秒と求められる．一方，弾性 板の1次モードの固有周期は，数値計算で用いたヤン グ率を使うと，約0.26秒と計算される．図–8に示さ れた弾性板の変動周期は，水面変動の周期に近いので，

その振動は主として水面変動により引き起こされたも のと考えられる．

なお，1ケースの計算時間は，単一CPU(Pentium4, 2.66GHz)のPCにより，約55分であり，そのうち有 限要素法の計算が占める計算時間は約51%であった．

今後は，有限要素法による計算精度の向上と，計算時 間の短縮が課題であると考えられる．

５．おわりに

多相場の解法であるMICSに，有限要素法による固 体変形計算を用いる解法を示した．固体計算にいくつ かの検討課題が残されるが，実験結果との比較による 解法の適用性はほぼ良好であったと考えられる．

(a) ケースh75

(b)ケースh125 図–8 弾性板先端の変位の時系列

参考文献

1) 泉元,谷口伸行,川田裕,小林敏雄. 円柱周りの3次元流 動解析(第3報，弾性支持円柱の場合). 日本機械学会論 文集(B編), Vol. 66, No. 644, pp. 1013–1020, 2000.

2) 田中嘉宏,森西晃嗣,松野謙一. 弾性変形を伴う物体周り の二次元非定常流れの数値計算.日本機械学会論文集(B 編), Vol. 72, No. 718, pp. 42–49, 2006.

3) 牛島省,山田修三,藤岡奨,禰津家久. ３次元自由水面流 れによる物体輸送の数値解法(3D MICS)の提案と適用 性の検討. 土木学会論文集, Vol. 810/II-74, pp. 79–89, 2006.

4) 牛島省,牧野統師,禰津家久. 四面体サブセル法を用いる 市街地に流入する氾濫流の3次元数値計算. 水工学論文 集, Vol. 51, pp. 787–792, 2007.

5) 牛島省,禰津家久. 陰解法を用いたコロケート格子によ る高次精度の流体解析手法の提案. 土木学会論文集, No.

719/II-61, pp. 21–30, 2002.

6) 牛島省,福谷彰,牧野統師,禰津家久. 3次元流体中を運 動する接触と変形を考慮した任意形状固体モデルの数値 解法. 応用力学論文集, Vol. 10, pp. 139–146, 2007.

7) E. Barragy and G. F. Carey. A parallel element-by- element solution scheme. Int. J. Num. Meth. Engrg., Vol. 26, pp. 2367–2382, 1988.

(2007.9.30受付)