修士論文発表要旨
(2011
年度)特性有限要素法を用いた
2
次元界面問題解析Analysis of Two Dimensional Interface Problem Using Characteristic Finite Element Method
土木工学専攻
32
号 戸来 友哉Tomoya HERAI
1
はじめにNavier-Stokes
方程式のような流体を扱う支配方程式には,移流項と拡散項が含まれており,これらの項によって 流れは変化を見せる.また,それらに対する解き方によっ ても数値解析結果は大きく変化している.近年,数値流体 解析において,移流項を高精度に計算するための
Cubic Interpolated Propagation / Constrained Interpolation Pro- file (CIP)
法をはじめとする, Multi Moment Schemeが提 案され数多くの分野で使用されている.また,高Reynolds
数問題で
Galerkin
法が不安定になることに対して,風上差分法, SUPG(streamline upwind / Petrov-Galerkin),特性 曲線などの手法が提案されている. 本研究では,特性曲 線に基づいた特性法に着目する.
CIP
法とは,双曲型微分方程式を解くための高次精度 差分法の1
つである.この手法の補間に対して1
次や2
次での関数で補間を行うと,数値拡散による数値振動等 が発生してしまい,解が不安定になってしまう. この問 題を解決するために,本研究では導関数値を考慮した, 3次の
Hermite
補間によって近似させる. さらに,非移流計算についても同様の要素を適用し, Galerkin法によっ て離散化を行う.本手法では非移流計算においても移流 計算で用いる
3
次要素を適用するため,精度が損なわれ ない.本研究では,非圧縮性
Navier-Stokes
方程式を用いて, 表面張力を考慮した2
次元界面問題に特性有限要素法 を適用する. 補間関数として,流速場には移流計算と同 様の3
次三角形要素,圧力場には1
次三角形要素を適用する.また,界面を
3
次のB-スプライン関数を用いて曲
線で表し,曲率の計算を行う.
2
支配方程式境界
Γ = ∂Ω
を有する2
空間領域Ω ∈ R 2
を定義し,領 域Ω
は2
つの部分領域Ω α = Ω α (t)(α = 1, 2)
によっ て構成され,非圧縮性流体で満たされているものとする.これらの
Ω α
は密度ρ = ρ α
と粘性係数µ = µ α
を有 している.また, 2つの領域の流体は混ざらなく,密度と 粘性係数は各流体で一定であると仮定し,界面をΣ
と する. 支配方程式には,以下の無次元化された非圧縮性Navier-Stokes
方程式を用い,各流速と圧力はu (α) , p (α)
と表す.ρ (α) ( ˙ u (α) i + u (α) j u (α) i,j ) + p (α) ,i
− µ (α) (u (α) i,j + u (α) j,i ) ,j = 0 in Ω (α) (1) u (α) i,i = 0 in Ω (α) (2)
また,界面Σ
において以下の条件を与える.{ ρ(u (1) n − u (2) n ) = 0 σ ij n i = σκn i
(3)
ここで,
σ
は表面張力係数,κ
は界面においての曲率を 表す.上式の第1
式は流速の連続性,第2
式は式(1)
を部 分積分した際に現れる,境界での応力ベクトルを表面張 力として考えており,以下のように表される.f i = σκn i (4)
3
特性法3.1
特性曲線に基づく近似u = u(x, t)
を既知の流速とする.X = X (t)
を時刻t
における粒子の位置とする.X
は以下の常微分方程式を 満たしている.dX
dt (t) = u(X (t), t) (5)
このとき,実質微分項はX
を使って,以下のように表す ことができる.( ∂u
∂t + uu ,i )(X(t), t) = d
dt u(X (t), t) (6)
∆t
を時間増分量として,差分近似したものは実質微分 項の近似となり,以下のようになる.u(X(t), t) − u(X(t − ∆t), t − ∆t)
∆t (7)
すなわち,これは起点と上流点との差分を表している.こ れが特性曲線に基づく近似の基本的な考え方である.
3.2
上流点の計算上流点
X
を求め,粒子の起点と差分をとることで実 質微分項を近似している. 図1
に示すように,節点l
の 位置をx l ,
要素e
の重心位置をx e ,
またそれぞれの上流 点の位置をX l n , X e n
とすると,上流点の位置は, 2次精度 の予測子・修正子法である次式によって求められる.X 1 (x) = x − u(x)∆t (8) X 2 (x) = x − (u(x)X 1 + u(x)) ∆t
2 (9)
式
(9)
の右辺第2
項は1
回目の上流点を求めた点にお いての流速である.図
1:
移流計算4
特性有限要素法4.1
移流計算本章で用いている, Navier-Stokes方程式に含まれてい る移流項に対して特性法を適用すると,以下のような近 似式を得ることができる.
ρ( ˙ u i + u j u i,j ) ≈ ρ u n+1 i (x) − u i (x)X n
∆t (10)
ここで,
∆t, n
はそれぞれ時間増分量,時間ステップを表 している.特性法による移流項の近似式
(10)
の右辺をゼロとし, さらに移流計算による解の更新をu ˜
とすると,移流計算 は以下のように行われる.˜
u(x) = u(x)X n (11)
また,節点l
での1
階導関数値の更新は,以下のように 行われる.∂˜ u i
∂x j = (δ jk − ∆t ∂u k
∂x j ) ∂u i
∂x k (X l n ) (12)
移流計算による解の更新をU ˜
にまとめると,以下のよ うになる.U ˜ = [
˜
u l , ∂˜ ∂x u | l , ∂ ∂y u ˜ | l , u ˜ e ] T
4.2
時間方向の離散化本手法は移流計算と非移流計算に分離し,あらかじめ 移流計算をした後,非移流計算を行う. 非移流計算にお いての,時間方向の離散化には,安定性に優れ,時間増分 を大きく取れる
Crank-Nicolson
法を適用する.ρ u n+1 i − u ˜ i
∆t + p n+1 ,i − 1
2 µ(˜ u i,j + u n+1 i,j + ˜ u j,i + u n+1 j,i ) ,j
= 0 (13) u n+1 i,i = 0 (14)
4.3
空間方向の離散化基礎方程式の重み付き残差方程式は以下のように表 される.
ρ
∫
Ω
ω i
u n+1 i − u ˜ i
∆t dΩ
+ 1 2 µ
∫
Ω
ω i,j (u n+1 i,j + ˜ u i,j + u n+1 j,i + ˜ u j,i )dΩ
−
∫
Ω
ω i,i p n+1 dΩ = −
∫
s
ω i pn i ds+µ
∫
s
ω i (u i,j +u j,i )n j ds
∫ (15)
Ω
qu i,i dΩ = 0 (16)
ここで,ω i , q
は重み関数である. 式(15)
の右辺項は, 部分積分した際に現れる,境界積分項である. これらは 式(3)
の第2
式により,表面張力として考える.補間関数として,流速場には
Hermite
型3
次三角形要 素,圧力場には1
次三角形要素を用いた,混合補間を行う.領域
Ω
における要素ごとの三角形領域Ω e
での,流速 場の有限要素近似をu h | Ω e
とすると,以下のように表さ れる.u h | Ω e =
∑ 3 i=1
(H 0i u i +H xi ∂u
∂x | i +H yi ∂u
∂y | i )+H 0e u e (17)
ここで,H 0i , H xi , H yi , H 0e
は補間関数であり,面積座 標L i
の関数によって以下のように表される.
H 0i = L 2 i (3 − 2L i ) − 7L 1 L 2 L 3
H xi = L 2 i (x ji L j − x ik L k ) − (x ji − x ik )L 1 L 2 L 3
H yi = L 2 i (y ji L j − y ik L k ) − (y ji − y ik )L 1 L 2 L 3
H 0e = 27L 1 L 2 L 3
また, 1次三角形要素による圧力場の有限要素近似
p h | Ω e
は以下のように表される.
p h | Ω e =
∑ 3 i=1
N i p i , N i = L i (18)
図
2: Hermite
型要素 図3: 1
次三角形要素有限要素方程式は以下のように示される.
M αβ (u n+1 βi − u ˜ βi ) + 1
2 D α,jβ,i (u n+1 βj + ˜ u βj ) (19) + 1
2 D α,jβ,j (u n+1 βi + ˜ u βi ) − B α,iλ p n+1 λ = f i (20) B λ,αi u n+1 αi = 0 (21)
5
曲率の計算法5.1 3
次B-スプライン関数の定義
界面の曲率を求めるために, 3次
B-スプライン関数を
用いて曲線を表現する. B-スプラインは,いくつかの制 御点を定義することで,滑らかな一本の曲線を得ること ができる. 本研究では,パラメトリック表現という手法 を紹介する.x
とy
の直接的な関係式ではなく,別の変数
(媒介変数)
を定義し,その変数とx,
その変数とy
の関係を独立して考える. 3次
B-スプライン関数を数式で
表すと,以下のようになる.x(t) =
n+2 ∑
j= − 2
N 3 (t − j)p x (j) ( − 1 ≤ t ≤ n + 1) (22)
y(t) =
n+2 ∑
j= − 2
N 3 (t − j)p y (j) ( − 1 ≤ t ≤ n + 1) (23)
B-スプライン曲線はいくつかの制御点の座標値に重み
づけの係数をかける,補間に似た考え方で曲線を近似し ている. この係数を重み関数で得るが,この関数は次数 によって一定である.N 3 (t) =
(3 | t | 3 − 6 | t | 2 +4)
6 ( − 1 < t < 1)
− ( | t |− 6 2) 3 ( − 2 < t < 2) 0 (t ≤ − 2, t ≥ 2)
(24)
以上の関数を定義し,曲線を近似することで界面を表現 し,曲率を求めることができる.
5.2
曲率の計算曲率
κ
はB-スプライン関数によって表現される,
パラメータ曲線を用いて求める.パラメータ曲線を用いた場 合の曲率は以下の式で求めることができる.
κ = x¨ ˙ y − y ˙ x ¨
( ˙ x 2 + ˙ y 2 ) 3 2 (25)
˙
x, x ¨
はそれぞれパラメータt
の微分,x ˙ = dx dt , x ¨ = d dt 2 x 2
を 示している.6
移動境界問題6.1
界面の移動方法界面に与えられる表面張力を用いて有限要素方程式 を解き,界面においての流速を得る. その流速を界面の 移動量として節点を移動させ更新し,リメッシングを行 い,次の解析へと移る.
6.2
リメッシングの方法初期の解析を行った後,図のように解析領域を
2
つに 分ける. その2
つは2
種類の流体によって分けられて いる.図
4:
解析領域(Ω 2 )
図5:
解析領域(Ω 1 )
解析の条件は表面張力のみを考慮して計算をするた め,界面周辺を中心に流速が発生する. その求められた 流速をメッシュの移動量として計算し, 2つの解析領域 の界面に適用し,メッシュを移動させていく. 本研究で は,要素分割の手法として, Delaunay分割法を適用して いる. 移動させた後,境界部分のメッシュをつなぎ合わ せてリメッシングが完了し,次のステップの解析対象と する.リメッシングに必要なデータは, 2つの解析領域の 節点,要素データとそれぞれの境界の節点番号があれば よい.
6.3
流速の補間リメッシング後の各節点の流速については,リメッシ ング後の節点がリメッシング前のどの要素内に含まれる かを探索し,探索された要素を構成する節点
(リメッシ
ング前の節点)からHermite
型3
次要素を用いて補間す る. そのため,リメッシングを行う前の解析領域を保存 しておく必要がある.U i =
∑ 3 i=1
(H 0i u i + H xi ∂u
∂x | i + H yi ∂u
∂y | i ) + H 0e u e (26)
U i
は更新された流速であり,u i
とH α
は更新前の要素 が持つ流速と補間関数である.7
数値解析7.1
数値解析例異なる
2
種類の流体を取り扱う問題として,正方形領 域においての表面張力による変形シミュレーションを 解析する.図7
で示す解析領域において,流速およびそ の1
階導関数値の境界条件を与える.図6
は有限要素分 割図であり,節点数は4,249,
要素数は8,352
である. ま た,界面まわりの節点数は128
である. それぞれの流体 の密度と粘性係数の値は,ρ 1 = 1.5, µ 1 = 1.0, ρ 2 = 1.0, µ 2 = 1.0
とする.図
6:
有限要素分割図 図7:
解析領域と境界条件7.2
数値解析結果X
Y
-1 -0.5 0 0.5 1
-1 -0.5 0 0.5 1
図
8: T=0.0
X
Y
-1 -0.5 0 0.5 1
-1 -0.5 0 0.5 1
図
9: T=0.5
X
Y
-1 -0.5 0 0.5 1
-1 -0.5 0 0.5 1
図
10: T=1.0
X
Y
-1 -0.5 0 0.5 1
-1 -0.5 0 0.5 1
図
11: T=1.5
解析結果より,表面張力の影響がはたらき,楕円型か ら徐々に流体が移動し,円形となっている. これは曲率 の大きい箇所の表面張力が大きくなり,円形になったと ころで一定になっていると見られる.
8
結論本研究では,特性有限要素法を用いて表面張力を考慮 した流体挙動解析を行った.運動方程式の境界積分項を 界面においての表面張力として扱った.界面の関数を
3
次
B-スプライン関数を用いて表現することができ,
曲率の計算を行うことができた. また,移動境界問題を解く ためのリメッシング手法についても示した. 以上より, これらの手法を用いて界面の挙動解析を行うことがで きた.
参考文献
