「移流輸送の 高次精度かつ高解像度スキームの開発」 山口大学大学院 理工学研究科 システム工学専攻 社会システム工学大講座 坪郷 浩一 山口大学大学院 坪郷 浩一
<研究の背景> 移流計算スキームの開発 数値振動を起こさない 数値拡散が少ない 従来,多くの移流項計算スキームが提案されているがいずれも一長一短あり決定版と呼べるものがない. 山口大学大学院 坪郷 浩一
<目的:移流方程式を数値的に高精度に解く計算手法の開発> 支配方程式(簡単のため一次元で考える) 対流形式,非保存形式 u=const.の場合に 両式は一致する. 保存形式 上式を用いて,従来の計算スキームの問題点を指摘する. Ex1 高次精度スキーム:非保存形式 長所:ガウス分布などの「極値の再現性」に優れている. 短所:「様々な勾配の値を含む分布」の再現性には改良の余地がある. Ex2 高解像度スキーム:保存形式 長所:「矩形分布での数値振動が抑制」され,再現性に優れている. 短所:ガウス分布などの「極値の再現性」には改良の余地がある. 山口大学大学院 坪郷 浩一
<研究の目的> 流れ解析では,基礎方程式中の移流項の計算スキームの離散化精度が計算解の精度を大きく左右する一つの要因である.本研究では,数値振動や数値拡散の影響の少ない「高次精度かつ高解像度スキームの開発」を目指す. Komatsuらが開発した6-point schemeに基づいて高次精度の6-point schemeを提案するものである.(第2章) 1.で提案した計算スキームを移動境界問題の数値解析に適用させ有効性の検証を行なう. (第3章,第4章) 山口大学大学院 坪郷 浩一
<本研究における高精度,高次精度,高解像度の違い> 高精度(high-accuracy): 数値解がより厳密解に近い場合を言う. 高次精度(high-order): スキームの打ち切り誤差のオーダーが高いことを意味する時に用いている. 高解像度( high-resolution): 数値振動が発生する箇所を滑らかな単調曲線分布として求め,それ以外では高精度な解を与えるようなスキーム いわゆるTVD条件を満足するTVD schemeである. 山口大学大学院 坪郷 浩一
高次精度かつ高解像度6-point schemeの開発(第2章) <各章の関連性> 移流計算スキームの開発(第1章) 数値振動を起こさない 数値拡散が少ない 高次精度かつ高解像度6-point schemeの開発(第2章) 気液界面の界面を捕捉する移流計算に第2章で開発した計算スキームを適用させる. 自由水表面流れ解析 密度関数法(第3章) Level Set法(第4章) 気液界面を精度よく捕捉するために,数値振動や数値拡散が発生しにくい移流項計算スキームを用いる必要がある. 山口大学大学院 坪郷 浩一
<第2章 高次精度6-point schemeの開発> 目的 Komatsuらが開発した1次精度の6-point schemeの導出の概念に従って,新たな高次精度の6-point schemeの構築を主要な目的とする.またフラックス制限関数(limiter)が導入できるように保存形式にも書き換え,さらに朝位らのdiscriminatorの改善を図り,高次精度の6-point schemeと組み合わせて様々な分布形状に対し「高次精度かつ高解像度移流計算手法を提案」する. ・特性曲線形式高次精度6-point scheme ・保存形式高次精度6-point scheme ・ 数値振動とクリッピングエラーの除去 ・ 多次元問題への適用 山口大学大学院 坪郷 浩一
<基礎方程式> 支配方程式(簡単のため一次元で考える) 対流形式,非保存形式 u=const.の場合に 両式は一致する. 保存形式 山口大学大学院 坪郷 浩一
<特性曲線法> を近隣の値を用いて内挿式で求める. これは特性曲線上でFの値が変わらないことを意味する. t n+1 n x ξ i-3 i-2 i-1 i i+1 i+2 on これは特性曲線上でFの値が変わらないことを意味する. 特性曲線表示 時刻n+1のi点の濃度を求める問題は,特性曲線によって時刻nの点ξの濃度を内挿によって求める問題に帰着される. を近隣の値を用いて内挿式で求める. 山口大学大学院 坪郷 浩一
<1次精度6-point schemeの最終形> ここで, a(=uDt/Dx)はクーラン数である. 山口大学大学院 坪郷 浩一
<Komatsuらが開発した1次精度の6-point schemeの問題点> 様々な勾配の値を含む半楕円分布や矩形分布の再現性には改良の余地がある ガウス分布などの極値 の再現性に優れている 山口大学大学院 坪郷 浩一
Fi-3,Fi-2,Fi-1,Fi,Fi+1,Fi+2 <6 point scheme (Komatsu et al,1985)> i-2 i-1 i i-3 i+1 i+2 i-4 i+3 3次曲線を内挿 Δx x 時刻nのi-3からi+2までの格子点上の濃度を Fi-3,Fi-2,Fi-1,Fi,Fi+1,Fi+2 とする. 境界条件:線形外挿 これら8個の濃度の値を使用する.ただし実質には6個の濃度の値しか使用していないことに注意する. 山口大学大学院 坪郷 浩一
<3次多項式> 4個の濃度値とそのx座標が与えられれば,その濃度値を通過する3次多項式F(x)を構築できる. i-2 i-1 i i-3 3次曲線を内挿 i-1/2 x 4個の濃度値とそのx座標が与えられれば,その濃度値を通過する3次多項式F(x)を構築できる. 山口大学大学院 坪郷 浩一
以下のような手順で6-point schemeを構築できる. i-1点 i点 i-1点の勾配 i-1/2点の勾配 i点の勾配 i-1 i i-2 i+1 i-1/2 x 以下のような手順で6-point schemeを構築できる. 6-point schemeの誘導で特徴的なのはi-1/2の位置の濃度勾配を使用していることである. 山口大学大学院 坪郷 浩一
<3次精度6-point schemeの導出> 4個の濃度値とそのx座標が与えられれば,その濃度値を通過する3次多項式F(x)を構築できる. 山口大学大学院 山口大学大学院 坪郷 浩一 坪郷 浩一
<Φi-1の勾配> 山口大学大学院 坪郷 浩一
<Φiの勾配> 山口大学大学院 坪郷 浩一
<Φi-1/2の勾配> 山口大学大学院 坪郷 浩一
a(=uDt/Dx)はクーラン数である. ここで, a(=uDt/Dx)はクーラン数である. 山口大学大学院 坪郷 浩一
Taylor級数展開による打ち切り誤差解析を行えば次式を得る.なお,打ち切り誤差項は4次までを記す. 右辺第1項は0次の誤差項である.この項は恒等的に0になる必要がある.恒等的に-2q +1=0となれば良いので重みq=1/2が得られる. 山口大学大学院 坪郷 浩一
<特性曲線形式3次精度6-point scheme> ここで, a(=uDt/Dx)はクーラン数である. 山口大学大学院 坪郷 浩一
<4次精度6-point schemeの導出> 4階微分は中央差分で近似 3次精度6-point scheme 山口大学大学院 坪郷 浩一
<特性曲線形式4次精度6-point scheme> ここで, a(=uDt/Dx)はクーラン数である. 山口大学大学院 坪郷 浩一
<モデル計算> 初期分布 計算格子幅 Δx=200m 計算時間間隔 Δt=100sec 流速 0.5m/sec クーラン数 0.25 クーラン数 0.25 初期分布 中心位置1400m,標準偏差264m,ピーク値10 中心位置2400m,標準偏差264m,ピーク値6.5 の二種類のガウス分布の重ね合わせ 幅2000mの矩形分布 中心位置5800mで最大値10の半楕円分布(空間方向の半径2000m) 山口大学大学院 坪郷 浩一
<モデル計算> Δx=200m , Δt= 100sec 1次精度6-point scheme(9600sec) CIP scheme(9600sec) 山口大学大学院 坪郷 浩一
<モデル計算> Δx=200m , Δt= 100sec 1次精度6-point scheme(192000sec) CIP scheme(192000sec) 山口大学大学院 坪郷 浩一
<モデル計算> 初期分布 クーラン数一定のまま空間解像度を倍にする. 計算格子幅 Δx=100m 計算時間間隔 Δt=100sec 流速 0.25m/sec クーラン数 0.25 初期分布 中心位置1400m,標準偏差264m,ピーク値10 中心位置2400m,標準偏差264m,ピーク値6.5 の二種類のガウス分布の重ね合わせ 幅2000mの矩形分布 中心位置5800mで最大値10の半楕円分布(空間方向の半径2000m) 山口大学大学院 坪郷 浩一
<モデル計算> Δx=100m , Δt= 100sec 1次精度6-point scheme(192000sec) CIP scheme(192000sec) 山口大学大学院 坪郷 浩一
<モデル計算> Δx=100m , Δt= 100sec 1次精度6-point scheme(192000sec) CIP scheme(192000sec) 山口大学大学院 坪郷 浩一
4次精度6-point scheme(192000sec) <モデル計算> Δx=100m , Δt= 100sec ガウス分布の最大値8.585 ガウス分布の最大値7.935 4次精度6-point scheme(192000sec) CIP scheme(192000sec)
<保存形式高次精度6-point scheme > 矩形分布で数値振動が発生している. 特性曲線形式高次精度6-point schemeにフラックス制限関数を導入する. なお,これ以後の特性曲線形式高次精度6-point schemeは,TVD化するために保存形式に直すため保存形式高次精度6-point schemeと呼ぶことにする. CIP schemeは保存形式に直すことができず,TVD化できない. 山口大学大学院 坪郷 浩一
保存形式の移流方程式を有限体積法で離散化 <保存形式高次精度6-point schemeへの拡張 > 保存形式の移流方程式を有限体積法で離散化 u :流速一定 離散化 整理すれば はコントロールボリュームの界面濃度だから5点の濃度で以下のように近似できる. まとめると a(=uDt/Dx)はクーラン数である. 山口大学大学院 坪郷 浩一
<本論文で提案する保存形式高次精度6-point scheme> a(=uDt/Dx)はクーラン数である. 山口大学大学院 坪郷 浩一
<モデル計算> 初期分布 計算格子幅 Δx=200m 計算時間間隔 Δt=100sec 流速 0.5m/sec クーラン数 0.25 クーラン数 0.25 初期分布 中心位置1400m,標準偏差264m,ピーク値10 中心位置2400m,標準偏差264m,ピーク値6.5 の二種類のガウス分布の重ね合わせ 幅2000mの矩形分布 中心位置5800mで最大値10の半楕円分布(空間方向の半径2000m) 山口大学大学院 坪郷 浩一
特性曲線形式1次精度6-point scheme <universal limiterの導入> Δx=200m , Δt= 100sec 特性曲線形式1次精度6-point scheme 保存形式1次精度6-point scheme universal limiter適用 山口大学大学院 坪郷 浩一
特性曲線形式1次精度6-point scheme <universal limiterの導入> Δx=200m , Δt= 100sec 数値振動は除去されているが,ガウスの極値の再現性の低下が問題 なお,極値の再現性の低下をクリッピングエラーと呼ぶ. 特性曲線形式1次精度6-point scheme 数値振動だけを除去でき,物理的に意味のある極値近傍ではlimiterを機能させない朝位ら(1999)のdiscriminatorを導入する. 保存形式1次精度6-point scheme universal limiter適用 山口大学大学院 坪郷 浩一
朝位らのdiscriminatorのアルゴリズム 朝位らのdiscriminatorの概念図 山口大学大学院 坪郷 浩一
universal limiterとdiscriminatorの適用 <朝位らのdiscriminatorの導入> Δx=200m , Δt= 100sec 保存形式1次精度6-point scheme universal limiter適用 保存形式1次精度6-point scheme universal limiterとdiscriminatorの適用 山口大学大学院 坪郷 浩一
universal limiterとdiscriminatorの適用 <朝位らのdiscriminatorの導入> Δx=200m , Δt= 100sec ガウスの極値の再現性は向上したが矩形分布の裾の近傍で数値振動が発生している. 保存形式1次精度6-point scheme universal limiter適用 保存形式1次精度6-point scheme universal limiterとdiscriminatorの適用 山口大学大学院 坪郷 浩一
<朝位らのdiscriminatorの誤作動> 改善 「勾配D1とD3が大きく異なる場合」においてlimiterを作用させない. 山口大学大学院 坪郷 浩一
discriminatorのアルゴリズム Maxは2個の引数の中で大きい方を選択する演算子である.dは分母が0にならないようにするために導入される非常に小さな値である.eは任意の閾値である. 修正された朝位らの discriminatorのアルゴリズム 山口大学大学院 坪郷 浩一
数値実験によりeは,3≦e ≦10が推奨される.本論文では, e=3を用いる. <修正された朝位らのdiscriminator> 移流問題を用いて提案した修正型discriminatorの機能の評価を保存形式1次精度6-point schemeを対象に閾値の検討を行った. 閾値とガウス分布のピーク値の関係 閾値と最小値の関係 数値実験によりeは,3≦e ≦10が推奨される.本論文では, e=3を用いる. 山口大学大学院 坪郷 浩一
<修正された朝位らのdiscriminator> Δx=200m , Δt= 100sec 特性曲線形式1次精度6-point scheme 保存形式1次精度6-point scheme universal limiter適用 保存形式1次精度6-point scheme universal limiterとdiscriminatorの適用 保存形式1次精度6-point scheme universal limiterと修正型discriminatorの適用 山口大学大学院 坪郷 浩一
<保存形式高次精度6-point scheme > universal limiterと修正型discriminatorの適用 山口大学大学院 坪郷 浩一
<保存性の検証> ここでCi (t)は計算終了時刻t におけるのi点の濃度,Ci (0)は初期のi 点の濃度,imax は計算点の総数である. 質量の相対誤差は,10-6%程度である,limiterやdiscriminatorを導入しても保存性は良好である. 山口大学大学院 坪郷 浩一
limiterやdiscriminatorを適用した. <多次元問題> ~剛体回転する流れ場におけるモデル計算~ 2次元平面で角速度ω=2π/12000(sec-1) 計算格子幅 Δx= Δy=100m 計算時間間隔 Δt=50sec 計算領域 4000×4000 初期分布 中心位置(1400m,1400m) 標準偏差200m,ピーク値10 のガウス型濃度分布 計算スキーム 特性曲線形式4次精度6-point scheme 保存形式4次精度6-point scheme QUICKEST scheme limiterやdiscriminatorを適用した. 山口大学大学院 坪郷 浩一
~剛体回転する流れ場における移流計算の結果~ (x=y上で) ピーク値の再現性は特性曲線6-point schemeが最も良い. 質量の相対誤差 保存形式4次精度6-point scheme:3.95×10-3% 特性曲線形式4次精度6-point scheme : 7.56×10-2% QUICKEST scheme:3.43×10-2% 質量の相対誤差は保存形式6-point schemeが最も良い. 山口大学大学院 坪郷 浩一
計算条件 計算スキーム <高次精度6-point schemeの非線形移流項> ~2次元キャビティ問題~ Ghia et al.の結果と比較できるように正方キャビティを計算対象とした.MAC法の一種であるFractional Step法で数値計算を行った. 計算条件 等間隔で50×50のセルに分割しstaggered格子を用いる. Re=1000 計算スキーム 保存形式3次精度6-point scheme 保存形式4次精度6-point scheme 5次精度風上差分法 limiterを適用した. 山口大学大学院 坪郷 浩一
y方向流速vのx方向分布 ~2次元キャビティ内の流動計算結果~ キャビティ中心を含む鉛直方向(y方向)および水平方向(x方向)の流速分布 x方向流速uのy方向分布 y方向流速vのx方向分布 山口大学大学院 坪郷 浩一
高次精度かつ高解像度6-point schemeの開発(第2章) 1. 極値の再現能力は,CIP scheme より優れている. 2.フラックス制限関数と修正されたdiscriminatorを用いれば,高精度かつ高解像度数値解を得ることが可能となった.一方,CIP schemeは保存形式に直すことができず,TVD化できない. 3. 高次精度6-point scheme はtime splitting法を用いることで 多次元問題にも容易に適用でき,精度の良い解が得られることが示された. 2.の操作が人為的な側面を持っている.今後の課題として開発した計算スキームのみで単調性や凹凸性の保存ができるような計算スキームの開発を検討する必要がある. 山口大学大学院 坪郷 浩一
高次精度かつ高解像度6-point schemeの開発(第2章) <適用事例(第3章,第4章)> 自由水表面流れ解析 気液界面を捕捉する関数の移流項に,提案した計算スキームを使用する. 界面捕捉法の有力な手法の 密度関数法 (第3章) Level Set法(第4章) を用いて提案した計算スキームの有効性を移動境界流れの検証モデルの1つであるダムブレイクモデルで検証する. 山口大学大学院 坪郷 浩一
<第3章 密度関数法を用いた自由水表面流れ解析> <第3章 密度関数法を用いた自由水表面流れ解析> 目的 本章では,密度関数の移流方程式に提案した高次精度差分スキーム(第2章)を適用させその有効性の検証を行う.また,計算進行に伴う界面のぼやけを回避し体積を補正する手法の開発について述べるものである. ・ 数値計算法 液相体積の計算 体積補正法 質量補正法 ・ 提案した計算スキーム(第2章)の検証 ・ まとめ 山口大学大学院 坪郷 浩一
・VOF法(1981,C.W.Hirtら):液相のみ <自由水表面が存在する流れの数値解析手法> 界面捕捉法 複雑で大規模な自由水表面問題に対しては界面捕捉法が有利である. 代表例 ・VOF法(1981,C.W.Hirtら):液相のみ 最近の研究:気液二相流を直接解く手法が発展している ・密度関数法(1996,Kanaiら) ・Level Set法(1994,Sussmannら) 山口大学大学院 坪郷 浩一
界面捕捉法の問題点:体積の保存性に問題がある. 既往の体積補正の手法として, Level Set法では姫野らの手法(1999年) VOF法では桜庭らの手法(2003年) が挙げられるが,密度関数法に関する体積補正の手法は,まだ開発されていない. 山口大学大学院 坪郷 浩一
< 目 的 > 本研究では,自由水表面を持つ流れ解析の有力な手法の一つである密度関数法を用いて以下のような解析・検討を行った. < 目 的 > 本研究では,自由水表面を持つ流れ解析の有力な手法の一つである密度関数法を用いて以下のような解析・検討を行った. <移流項計算スキームの検証> 本章では,密度関数の移流方程式に提案した高次精度差分スキーム(第2章)を適用させその有効性の検証を行う. <体積補正の手法の提案> ① この手法のための「オリジナルの体積補正法」を開発する. ② ①で必要な液相体積を求めるアルゴリズムを開発する. 山口大学大学院 坪郷 浩一
<基礎方程式> 連続の式 運動方程式 密度の保存式 ここでuは流速ベクトル,pは圧力, μは粘性係数,Fは体積力である. 山口大学大学院 山口大学大学院 坪郷 浩一
<密度の保存式> 密度関数Φ(0<φ<1)の判別 気相Φ=0,界面Φ=0.5,液相Φ=1 密度および粘性係数 密度の保存式を直接解く代わりに密度関数(0<φ<1)の保存則を解く 密度関数Φ(0<φ<1)の判別 気相Φ=0,界面Φ=0.5,液相Φ=1 密度および粘性係数 ここに ρLip は液相の密度 は気相の密度 は気相の粘性係数である. ρGas μLip μGas , は液相の粘性係数 山口大学大学院 坪郷 浩一
<体積補正法の概念> 界面セルの液相占有率が真のそれよりも小さいために液相体積が過小評価する 界面セルの液相占有率を高くする必要がある. 界面の総体積 真の界面 界面セルの液相占有率を高くする必要がある. 液相体積誤差量VErr 計算界面 液相占有率の補正量を算出し,それを各界面セルに加えることで体積誤差を減少させる. 山口大学大学院 坪郷 浩一
<界面セルの把握状態の一例> 気液界面セル 気相セル 液相セル Φ=0.5の等高線 <液相体積(面積)を求める> 山口大学大学院 坪郷 浩一
<体積補正の手法> 現在時刻tの液相体積誤差量VErr(t)を求める 体積補正量LErrを求める VLiq(t):任意時刻tの液相体積 Vini:初期液相体積 体積補正量LErrを求める A(t):任意事項tの計算領域内の界面セルの総体積 Φn+1を求める時点で得られた密度関数に次のような操作を行う. 山口大学大学院 坪郷 浩一
「体積補正の後に全質量の保存」を保証する目的のルーチンを作成 <質量補正の導入> 体積補正法の導入により「計算領域内の全質量保存を破綻」させてしまう可能性がある. 原因 質量誤差は界面・界面近傍セルのΦ値が適切でないために生じるものと考えられる. 改善点 「体積補正の後に全質量の保存」を保証する目的のルーチンを作成 山口大学大学院 坪郷 浩一
計算領域内のセルから質量誤差を一様に差し引く. <質量補正の手法> 現在時刻tの全質量誤差量MErr(t)を求める M(t):任意時刻tの全質量 Mini:初期全質量 質量補正量MErrを求める. MErr(t):任意時刻tの質量誤差 AM(t):界面近傍の0<Φ<1となる 領域の任意時刻tの体積 計算領域内のセルから質量誤差を一様に差し引く. 山口大学大学院 坪郷 浩一
<計算のアルゴリズム> 山口大学大学院 坪郷 浩一
<移流項計算スキームの検証> 空気 水 密度:1.25kg/㎥ 粘性係数:1.50×10-5Pa・s 密度:1000kg/㎥ 表面張力なし L=0.15m,0.6m×0.6m 計算領域40×40 計算格子間隔Δx=Δy = 0.015m 時間刻みは Δt = 0.0001sec 山口大学大学院 坪郷 浩一
<移流項計算スキームの検証> 体積補正法なしで最も数値拡散が少ない組み合わせの検討 山口大学大学院 坪郷 浩一
<全質量の相対誤差> 山口大学大学院 坪郷 浩一
<液相体積の相対誤差> 液相体積・全質量の相対誤差からTVD-3C6Pの組み合わせが最も数値拡散が少ない. 山口大学大学院 坪郷 浩一
<液相挙動 静止画> TVD-3C6P 山口大学大学院 坪郷 浩一
<液相挙動 動画> TVD-3C6P 体積補正法なし 体積補正法あり 山口大学大学院 坪郷 浩一
<水柱先端位置> TVD-TVD TVD-3C6P 山口大学大学院 坪郷 浩一
<液相体積と質量の相対誤差> TVD-3C6P 全質量の相対誤差 液相体積の相対誤差 山口大学大学院 坪郷 浩一
第2章で提案した計算スキームの有効性の検証結果を以下に示す. <第3章のまとめ> 第2章で提案した計算スキームの有効性の検証結果を以下に示す. 体積補正法を施さない場合 界面を補足する移流計算に提案した計算スキーム(3次精度保存形式6-point scheme)を適用したところ液相体積の保存能力が他の計算スキームと比べて高い,それでも十分な保存性が得られない. 体積補正法を施した場合 本研究で提案する体積補正法の適用により体積保存と質量保存の精度が向上することを確認した.しかしながら,流動が抑制される傾向も現れる.この問題は密度関数の輸送方程式に計算精度の良い提案した計算スキーム(第2章)を用いることである程度回避することができる. 山口大学大学院 坪郷 浩一
<まとめ> 1. 極値の再現能力は,CIP scheme より優れている. 2. フラックス制限関数と修正されたdiscriminatorを用いれば,高精度かつ高解像度数値解を得ることが可能となった. 3. 高次精度6-point scheme はtime splitting法を用いることで 多次元問題にも容易に適用でき,精度の良い解が得られることが示された. 4. 自由水表面流れ解析の一つである密度関数法の輸送方程式の移流項に提案した計算スキームの適用事例を以下に示す.なお,モデル計算にはダム破壊現象を用いた. (1)体積補正法を適用しない場合 従来の計算スキーム比べて界面近傍のぼやけや液相体積の減少をある程度抑制している. (2)体積補正法を適用した場合 従来の計算スキームでは,流体運動が抑制される傾向があった.しかし,提案した計算スキームを用いることで流体運動の抑制をある程度回避することができた. 山口大学大学院 坪郷 浩一