講義予定 案. 9/ 数値シミュレーションの手続き テキスト第 章. 9/ 9 偏微分方程式と解析解 テキスト第 章 3. 9/6 休講 4. 9/30 差分方程式とそのスキーム テキスト第 3 章 変換 テキスト第 4 章 5. 0/ 7 計算 テキスト第 5 章 連立一次方程式の解法 テキスト第 6 章 6. 0/ 流れ関数 ポテンシャルによる解法 テキスト第 7 章 7. 0/8 流速 圧力を用いた解法 テキスト第 7 章 8. / 4 熱流体解析と多相流解析 9. / 乱流の数値解析 b 金子暁子先生 0. /8 数値解析の実際 b 渡辺正先生 JAEA. /5 予備日
圧力 流速場 を用いる解法 Los Alamos 国立研究所グループ T-3 米国 MAC 法 Make ad Cell method: Hallow ad Welch965; Welch et al.966 SMAC 法 Smlfed Make ad Cell method: Amsde ad Hallow970 S 法 ctoal Ste method または Two Ste Poecto method: Kothe et al. 99 HSMAC 法 Hghl Smlfed MAC method または SOLA 法 : Ht et al. 975 ICE 法 Hallow ad Amsde 97 ALE 法 Ht et al. 974 O 法 Ht ad Nchols 98 Imeal Collage 英国 SIMPLE 法
基礎方程式 : 連続の式と Nae-Stokes 方程式 基礎方程式 0 τ t g --------- --------- 差分式 流速と圧力のみを陰的に取り扱う 0 t τ g --------- 3 --------- 4
MAC 法 0 g t τ 差分式 4 の両辺の発散 degece をとる 連続の式 3 を満足する圧力をとすると このとき であるから g t τ 時刻 において 右辺は全て既知であるから 上式はについてのポアッソン方程式となる 一旦が求まると ] [ g t τ より が決定できる
MAC 法の解法手順初期値 : ポアソン方程式よりを求める 速度を求める g t τ ] [ g t τ
SMAC 法 g t τ δ 上式の両辺の発散 degece をとる 連続の式が満足される圧力をとすると 0 P δ 求めるべき圧力場をと分解すると 上式を 陽的部分と陰的部分に分解する g t τ ~ ~ t δ ~ t δ ~ t δ 陰的部分 陽的部分
SMAC 法の解法手順 初期値 : ~ を求める ポアソン方程式より ~ τ t δ δ ~ t を求める g 速度 δ を求める を求める ~ t δ
S 法 Two Ste Poecto 法 上式の両辺の発散 degece をとると 連続の式が満足されているとするとであるから 0 このポアッソン方程式を解くことによってが求まり 最終的に 速度場が求められる を 直接陽的部分と陰的部分に分解する g t τ ~ ~ t ~ t ~ t 陰的部分 陽的部分 g t τ
S 法の解法手順 初期値 : ~ を求める ~ t τ g ポアソン方程式より を求める t ~ 速度 を求める t ~
HSMAC 法 SOLA 法 セル における発散 D を定義する D { } { } これは 変数 が質量保存式を満足していない程度を表す指標となりうる D<0ならば このセルに向かって正味で質量の流入が生じることに対応し このセルにおける圧力を増加させる必要がある 逆に D>0ならば このセルより正味で質量の流出が生じることに対応し このセルにおける圧力を減少させる必要がある 故に 圧力を調整してDを0に調節すればよいことがわかる
圧力の補正 セル における発散 D をそのセルの圧力 の関数とみなす つまり D D である D0の根を求めるためにNewto 法を用いると δ D k k k / D ここで kは第 k 番目の反復を表す まず質量保存式をDとおき を求めてから同式に代入する 求められたD を で偏微分すれば次式のようになる k δ D t
流速の更新 この式にD k を代入してδ k が計算されると 発散を0とすべく流速が更新される t k δ / k k k t k δ / t k δ / k k tδ / k k k k ここで加速係数を取り入れて 早く収束させる δ k k ω D t ω <
対流項の差分対流項 f f f f f 4 α α f 4 α α
対流項の差分 f 4 α α f 4 α α ここで重み a は α0 のときは中心差分となり差分近似の精度が高くなり α の時には風上差分となり安定性が確保されやすくなるという効果をもつ
ドナー セル差分 風上差分 変化の激しい現象を扱う場合 0 0 0 0 / / / / < > < > またはであった場合 数値的不安定性が発生しやすい たとえば 0 0 / / / / < t t とする つまり流速の向きを考えてドナーが風上にいるようにする 上式では - がドナーで がアクセプターとなる この差分をドナー セル差分という
粘性項の差分 s f ν 粘性項 s s f f s f ν
スタガード メッシュ セル に対しして 圧力をセル中心 流速はセルエッジで定義する スタガード メッシュ 質量保存式を陰スキームで近似し 0 t - Nae-Stokes 方程式を陽スキーム近似すると g f f f s t g f f f s -
スタガード メッシュの物理的意味 右図のように物理的考察によって初めから格子点あるいはメッシュ系をずらすものがある 圧力 密度 温度 T はセルの中心 で定義し 流速 はセルエッジ / で定義する このようなメッシュ系をスタガード メッシュという - -/ / T この系では物理学的直感によく訴える 圧力 - と の圧力差によって流体運動が生じ -/ が駆動される また セル中心の密度は -/ で運び込まれる質量と / で運び出される質量によって蓄積量が定まると考えたのである
HSMAC 法の解法手順 初期値 : ~ D δ k k 速度 を求める を求める を求める を求める k D δ とする g f f f t s t g f f f s D k k k / k k k tδ / k k t k δ / k k t k δ / t k δ / k k D k
カルマン渦列の計算結果 計算条件 tme0.0 sec tme5.0 sec Ilet eloct 5.0 cm/s kematc scost 0.5 cm/s 0. 0.4 IMAX0 JMAX40 tme0.0 sec tme5.0 sec tme0.0 sec
O 関数の輸送方程式 関数 の時間空間での変化をセル中心で定義し 流速 は セルエッジで定義するスタガード メッシュを採用すると O 関数の輸送方程式 t 0 の差分近似は次式で与えられる t / / / / / / / / 上式を解くことによって 関数 の輸送が数値的に計算される
O 法 - O 関数 自由表面を記述する方法の つに O 法がある これは なるセル体積に占める流体体積の比率を O 関数または単に 関数と定義して流体の領域を記述するものである 図中のセル IJ では流体占有率は 0.7 でセル I- では 0.3 であると読む 表現の便利のために 0.9-0.3 のように記す - - 0 0. 0.5 0.3 0.9 0.7 -
界面位置 - - T R - L P R - B - -
O 関数による界面位置の決定 方向の自由表面位置 L L 方向の自由表面位置 L L 表面こう配の差分近似 / / / / P R L P これらを内分すると次式のようになる / / / / / / L P P R d d
ドナーアクセプター法 関数 を数値的に解く工夫のひとつとしてドナーアクセプター法がある ドナーとアクセプターにおける気相と液相の割合をそれぞれ α D α A とする セルエッジ / の密度を次式のようにする / l gα g A l α D 0 α A 0 < α D < α ここで 気相の密度を g 液相の密度を l とした D α D D / α A A
ドナーアクセプター法 セルエッジの流速をと表し ドナーをD アクセプターをAと表す Dセルには全体 D δ D の流体がある これを斜線を引いた長方形とする δt 時間にセルエッジの単位面積あたりを通過してDからAに輸送される体積 は D δ D A δt である これは点線で囲んだ領域内である δt 内にDからAに 輸送される変数 は AD δt AD A したがって D セルの流体は全て A セルに輸送されずに D セルに流体が残る D δt δt A
ドナーアクセプター法 右図のように A > D δ D の場合 D δ D よりも多い流体を D から A に送ることができない そこで 輸送量は D - A A AD δ t D δ D D δ D A となる
ドナーアクセプター法 右図の場合 D セルのボイド保有量 - D δ D よりも多い - A のボイド量が A セルに伝播され 物理的解釈から好ましくない したがって 次式で求められる流体の量 C A D δd - D δ D C- A -- D δ D を追加して A セルに輸送する つまり δt 時間内に D から A に伝播される変数 の量は AD δt A である C D - A A A
ドナーアクセプター法 以上の式をまとめると AD m[ C δ ] / t sg / AD δ D D C / ma δt [ δ 0] AD D D ここで 式の演算 m は D セルが保有する流体 D δ D 以上の流体が輸送されるのを防ぎ この式の演算 ma は D セルが保有するボイド以上のボイドが輸送されるのを禁止している
壊れたダムの流れ 0 K 3.8 s 0 W L W R 0 5 W S O 0
壊れたダムの流れ 非圧縮性の流体の運動方程式を速度 と圧力 で記述すると 基礎式は質量保存式ならびにNae-Stokes 方程式である 図に示すように長方形 0 5 0 3.8 のように水が蓄えられている t0において板 Kを瞬時に取り除く 水は形状を変化させながら下流の障害物 Oを乗り越えて右壁 W R に到達する 自由表面の変化に注目しつつ水の運動を求めよ この図の初期条件と境界条件を付帯させて質量保存式ならびにNae-Stokes 方程式を解け ただし壁面 WL WB WR と障害物 O 壁面はすべりなしの壁とする Pogam: BO
MAC 法 マーカーセル法 右図に示すようにスタガードメッシュ系とし セルエッジで流速が定義されるとする セルの中に質量のない粒子を配置する この粒子は流速とともに移動する このマーカー粒子の位置と動きを点や線で描くと条痕線が得られる その結果として表面形状が記述できる と考えるものである / -/ P 3 P P -/ /
MAC 法 マーカーセル法 マーカー粒子の移動の計算は セルエッジで定義された流速を内挿し 個々の粒子の速度を定めて時間積分することである 時間 ttにおける第 番目の粒子の位置 は 時 0 刻 t0の初期位置 から次式によって求められる 0 t t
MAC 法 マーカーセル法 右図を元にしてマーカー粒子の流速 の内挿を考 える 点 P の周辺に の大きさの領域を設定し その分割領域 A A A 3 A 4 を重みとして次式のように 方向の線形補完して A A A 3 3 A 4 / -/ A -/ A 4 A 4 4 3 / A のように が決定される
マーカーセル法におけるセルの定義 E E E E E E E B E E E E B E B E E B B E E B E :ll cell E:Emt cell B :Boda cell
水滴の衝突 右図に示すように点 X C Y C に半径 R の水滴が速度 0 で厚さの小さい静水に衝突する 自由表面の変化に注目しつつ過渡状態を求めよ ただし対称の条件とし 左壁 底面はすべりなし壁とする Y C W L LHT Pogam: DROP SMAC W B 0 X C