講義予定. 第 回目 阿部数値シミュレーションの手続き. 第 回目 9 阿部差分方程式と差分スキーム. 第 回目 6 田中方程式の代数化 連立 次方程式の解法. 第 回目 田中並列計算法 5. 第 5 回目 阿部 MC 法など差分の計算方法 6. 第 6 回目 田中有限要素法 7. 第 7 回目 田中有限要素法 8. 第 8 回目 田中有限要素法 N-S プログラムによる実習 9. 第 9 回目 阿部乱流 熱流体 多相流. 第 回目 8 田中津波解析など事例紹介. 第 回目 5 試験予定日 流体運動の基礎方程式質量保存式 次元非圧縮性流れの Ner-Soes 方程式 密度動粘性係数圧力 : : : 方向流速 方向流速 : : 渦度と流れ関数渦度 渦度輸送方程式 流れ関数 コーシーリーマンの関係式 流れ関数 についての Posso 方程式渦度輸送方程式の導出非保存系表示の Ner-Soes 方程式 : 渦度輸送方程式上式を で微分した式から 下式を で微分した式を差し引く : 渦度 ここで
流れ関数についての Posso 方程式の導出を 以下の渦度の定義式に代入する 流れ関数の定義式 コーシーリーマンの関係式 流れ関数 についての Posso 方程式より 渦度移動方程式を 次元とし として書き換えると渦度移動方程式 Re Re 対流拡散方程式となる 次元渦度移動方程式圧力方程式保存系表示の Ner-Soes 方程式 : 圧力 についての Posso 方程式上式を で 下式を で微分して足し合わせる S S : 発散圧力についてのポアソン方程式ここで S S S 圧力 についての Posso 方程式を流れ関数を用いて表すと
計算フローチャート速度場コーシーリーマンの関係式圧力圧力に関するポアソン方程式流れ関数表示のソース項流れ関数流れ関数に関するポアソン方程式渦度渦度移動方程式流速分布 S S S 計算手順. == での初期値から初めて = まで計算が完了しているとする したがって 変数値 は既知である. 渦度移動方程式にしたがって 内点の新しい変数値 + を計算する. 内点での新しい + を用い 流れ関数方程式にしたがって内点の新しい + を反復収束過程から計算する. 新しい値 + を用いて新しい流速 をから計算する. 内点における新しい + + を用い 境界条件から新しい境界値を計算する 下図に示したように境界 6 がある 領域 は流れに置かれた障害物である 障害物後方の流れを計算しなさい 問題設定 5 6 基礎方程式と基本関係式バックステップ流れを渦度 と流れ関数 で記述すると基礎方程式は である
境界条件 は中心線とし と 5 はそれぞれ障害物の上面と底面とする は上方境界 は上流境界 6 は流出境界と呼ばれる 境界 - 5 - は流線であるから=Cosである 記述の便利のために 5 と書く 上流境界 では一様な流速 で流れが体系に入る つまり次の条件である 6 5 境界条件 障害物の上面 ではすべりなしの壁とする つまり のうえでは 同じように障害物の底面 5 では である したがって 境界 の上では流線の条件式より 6 5 境界条件 境界 は中心線であるから における速度の条件より 上方境界 もすべりなしの壁とした この境界も流線であるから =Cos である における流線の条件より また における条件より 6 5 したがって 計算フローチャート コーシーリーマンの関係式 速度場 流速分布 渦度移動方程式 渦度 流れ関数に関するポアソン方程式 流れ関数 流れ関数表示のソース項 圧力に関するポアソン方程式 圧力 S S S
渦度移動方程式の差分表現 渦度移動方程式のTCS 近似 f f s s f f s s SOR 法による Posso 方程式の数値解析 流れ関数 についての Posso 方程式 この式は連立一次方程式となる これを SOR 法で解くこととする つまり 第 回目の反復値を とすると Posso 方程式の差分スキーム 6 Posso 方程式の差分スキーム 6 Posso 方程式 空間に中心差分近似を用いた差分方程式 Nrl order: を 9 と番号を付け直す 5 7 6 7 8 9 5 6 9 = とする 5 8
Posso 方程式の差分スキーム 6 Nrl order によって差分式は以下のように表すことができる 9 8 7 6 5 9 8 7 7 9 6 5 6 8 5 9 8 7 6 5 Posso 方程式の差分スキーム 6 以上の置き換えによって 行列式は 以下の形に書き換わる T T T T T T 9 8 7 6 5 9 8 7 6 5 : ここで疎な行列 Srse Mr 変数の数 : = 格子分割の数 : 分割とすると係数行列の全要素数 : = 非零の要素数の割合 : 5- 非零の要素数 : +- + - =5- = % = 75 % =5 6.8 % =.6 % =. % =5.968 % =.96 % SOR 法 Gss-Sedel 法を書き その左辺を仮の反復値と見る すなわち と定義する 反復値は重み として あるいは変形して は緩和係数と呼ばれる
SOR 法 > を過緩和といい < を不足緩和という 上式からを消去すると この表現を SOR 法という 行列 を分離すると この式のベクトル形式は である SOR 法 と定義するとこの式は L I L となる この式を整理して L I L L L L が反復法の行列形式による表現である は SOR 行列と呼ばれ は Gss-Sedel 行列を意味している L L SOR 法 よって SOR 法は残差修正法である この式中において R とし 右辺第 項の [ ] の中の 項は とかける こうすると上式は R Prorm SOR 境界条件の差分表現 すべりなしの壁を考える つまり 次の条件 が指定される壁である + の流れ関数 + を点 で Tlor 展開すると である 壁面 はひとつの流線を形成し =cos = である - 5 - も流線であって この場合 -5- = とした
境界条件の差分表現 すべりなしの壁では となる しかし 方向には流速 は分布をもつ つまり である ところで 壁面での渦度 は であるから となる すべりのない壁では = = であるが 渦度の生成 はある 境界条件の差分表現 これらの式を流れ関数を Tlor 展開した式に代入すると となる 一般に壁面を点 w とすると この点の境界条件は である ところで 壁面での渦度 は であるから あるいは w w w となる すべりのない壁では = = であるが 渦度の生成 はある ここで は壁面に垂直にはなれた次の格子点までの距離である 境界条件の差分表現 境界条件の差分表現 5 鋭い角 について考える 鋭い角の流れ関数 の評価は心配ない 鋭い角の渦度 の指定には以下の式が用いられる 不連続近似 + 対称近似 壁面平均値近似 + 右図のように格子状に壁を設定しないが点 + は壁として表現可能である この境界をすべりのない壁とする つまり + = である 図中の + に対して境界点 について とおけば点 + では + = となってすべりのない壁の近似となる + + - - + +
時間ステップ の制御法 Cor 条件 c 拡散数の条件 : 拡散係数 変係数の熱伝導方程式 TCS スキームの安定条件 の決定法 計算フローチャート コーシーリーマンの関係式 速度場 流速分布 渦度移動方程式 渦度 流れ関数に関するポアソン方程式 流れ関数 流れ関数表示のソース項 圧力に関するポアソン方程式 圧力 S S S 圧力 流速場 を用いる解法 Los lmos 国立研究所グループ T- 米国 MC 法 Mrer d Cell mehod: Hllow d Welch965; Welch e l.966 SMC 法 Smlfed Mrer d Cell mehod: msde d Hllow97 S 法 rcol Se mehod または Two Se Proeco mehod: Kohe e l. 99 HSMC 法 Hhl Smlfed MC mehod または SOL 法 : Hr e l. 975 IC 法 Hllow d msde 97 L 法 Hr e l. 97 O 法 Hr d Nchols 98 Imerl Colle 英国 SIMPL 法 基礎方程式 : 連続の式と Ner-Soes 方程式 基礎方程式 --------- --------- 差分式 流速と圧力のみを陰的に取り扱う --------- ---------
MC 法 差分式 の両辺の発散 derece をとる 連続の式 を満足する圧力をとすると このとき であるから 時刻 において 右辺は全て既知であるから 上式はについてのポアッソン方程式となる 一旦が求まると ] [ より が決定できる MC 法の解法手順初期値 : ポアソン方程式よりを求める 速度を求める ] [ SMC 法 上式の両辺の発散 derece をとる 連続の式が満足される圧力をとすると P 求めるべき圧力場をと分解すると 上式を 陽的部分と陰的部分に分解する 陰的部分 陽的部分 SMC 法の解法手順を求める 速度を求める 初期値 : を求める ポアソン方程式よりを求める
S 法 Two Se Proeco 法 上式の両辺の発散 derece をとると 連続の式が満足されているとするとであるから このポアッソン方程式を解くことによってが求まり 最終的に 速度場が求められる を 直接陽的部分と陰的部分に分解する 陰的部分 陽的部分 S 法の解法手順速度を求める 初期値 : を求める ポアソン方程式よりを求める これは 変数が質量保存式を満足していない程度を表す指標となりうる < ならば このセルに向かって正味で質量の流入が生じることに対応し このセルにおける圧力を低下させる必要がある 逆に > ならば このセルより正味で質量の流出が生じることに対応し このセルにおける圧力を増加させる必要がある このように 圧力を調整して を に調節すれば 最終的な速度場が得られることがわかる HSMC 法 SOL 法 セル における発散 を定義する 差分式 流速と圧力のみを陰的に取り扱う HSMC 法 SOL 法 アルゴリズム圧力変化分が求められれば 速度場が求められる
流速の更新 が計算されると 速度場は よりとすると 流速の更新 速度場の評価値 を発散の式に代入すると 圧力の補正 いま セル における発散 をそのセルの圧力 の関数とみなすことにする つまり = である = の根を求めるために Newo 法を用いると ここで は第 番目の反復を表す まず質量保存式を とおき を求めてから同式に代入する 求められた を で偏微分すれば次式のようになる 流速の更新ここで加速係数を取り入れて 早く収束させる この式に を代入して が計算されると 発散を とすべく流速が更新される
対流項の差分対流項 f f f f f f 対流項の差分 f f ここで重み は = のときは中心差分となり差分近似の精度が高くなり = の時には風上差分となり安定性が確保されやすくなるという効果をもつ ドナー セル差分 風上差分 変化の激しい現象を扱う場合 またはであった場合 数値的不安定性が発生しやすい たとえば とする つまり流速の向きを考えてドナーが風上にいるようにする 上式では - がドナーで がアクセプターとなる この差分をドナー セル差分という 粘性項の差分 s f 粘性項 s s f f s f
スタガード メッシュ スタガード メッシュの物理的意味 セル に対しして 圧力をセル中心 流速はセルエッジで定義する スタガード メッシュ 質量保存式を陰スキームで近似し - Ner-Soes 方程式を陽スキーム近似すると f f f s f f f s - 右図のように物理的考察によって初めから格子点あるいはメッシュ系をずらすものがある 圧力 密度 温度 T はセルの中心 で定義し 流速 はセルエッジ + で定義する このようなメッシュ系をスタガード メッシュという - - + + T この系では物理学的直感によく訴える 圧力 - と の圧力差によって流体運動が生じ - が駆動される また セル中心の密度は - で運び込まれる質量と + で運び出される質量によって蓄積量が定まると考えたのである HSMC 法の解法手順 初期値 : を求める f f f s f f f s カルマン渦列の計算結果 計算条件 : Ile eloc = 5. cms emc scos =.5 cm s Δ=. Δ=. IMX= JMX= me=. sec me=5. sec me=. sec me=5. sec me=. sec を求める を求める 速度 を求める とする
O 関数の輸送方程式 O 法 - O 関数 関数 の時間空間での変化をセル中心で定義し 流速 は セルエッジで定義するスタガード メッシュを採用すると O 関数の輸送方程式 の差分近似は次式で与えられる 上式を解くことによって 関数 の輸送が数値的に計算される 自由表面を記述する方法の つに O 法がある これは なるセル体積に占める流体体積の比率を O 関数または単に 関数と定義して流体の領域を記述するものである 図中のセル IJ では流体占有率は.7 でセル I- では. であると読む 表現の便利のために =.9 - =. のように記す - -..5..9.7 - + 界面位置 O 関数による界面位置の決定 方向の自由表面位置 - - T 方向の自由表面位置 - L P R - + - R - + 表面こう配の差分近似 P L R P これらを内分すると次式のようになる d d P L R P
ドナーアクセプター法 関数 を数値的に解く工夫のひとつとしてドナーアクセプター法がある ドナーとアクセプターにおける気相と液相の割合をそれぞれ とする セルエッジ + の密度を次式のようにする l l ここで 気相の密度を 液相の密度を l とした + + ドナーアクセプター法 セルエッジの流速をと表し ドナーを アクセプターをと表す セルには全体 の流体がある これを斜線を引いた長方形とする 時間にセルエッジの単位面積あたりを通過してからに輸送される体積 は である これは点線で囲んだ領域内である 内にからに 輸送される変数 は したがって セルの流体は全て セルに輸送されずに セルに流体が残る ドナーアクセプター法 右図のように > の場合 よりも多い流体を から に送ることができない そこで 輸送量は となる - ドナーアクセプター法 右図の場合 セルのボイド保有量 - よりも多い - のボイド量が セルに伝播され 物理的解釈から好ましくない したがって 次式で求められる流体の量 C - C=- -- - を追加して セルに輸送する つまり 時間内に から に伝播される変数 の量は C である
ドナーアクセプター法 壊れたダムの流れ 以上の式をまとめると m C s K C m.8 W L s = W R ここで 式の演算 m は セルが保有する流体 以上の流体が輸送されるのを防ぎ この式の演算 m は セルが保有するボイド以上のボイドが輸送されるのを禁止している 5 W S O 壊れたダムの流れ 非圧縮性の流体の運動方程式を速度 と圧力 で記述すると 基礎式は質量保存式ならびにNer-Soes 方程式である 図に示すように長方形 5.8のように水が蓄えられている =において板 Kを瞬時に取り除く 水は形状を変化させながら下流の障害物 Oを乗り越えて右壁 W R に到達する 自由表面の変化に注目しつつ水の運動を求めよ この図の初期条件と境界条件を付帯させて質量保存式ならびにNer-Soes 方程式を解け ただし壁面 WL W WR と障害物 O 壁面はすべりなしの壁とする MC 法 マーカーセル法 右図に示すようにスタガードメッシュ系とし セルエッジで流速が定義されるとする セルの中に質量のない粒子を配置する この粒子は流速とともに移動する このマーカー粒子の位置と動きを点や線で描くと条痕線が得られる その結果として表面形状が記述できる と考えるものである + - - P P P + Prorm: O
MC 法 マーカーセル法 マーカー粒子の移動の計算は セルエッジで定義された流速を内挿し 個々の粒子の速度を定めて時間積分することである 時間 =における第 番目の粒子の位置 は 時 刻 =の初期位置 から次式によって求められる MC 法 マーカーセル法 右図を元にしてマーカー粒子の流速 の内挿を考 える 点 P の周辺に の大きさの領域を設定し その分割領域 を重みとして次式のように 方向の線形補完して + - - + のように が決定される マーカーセル法におけるセルの定義 水滴の衝突 :ll cell :m cell :odr cell 右図に示すように点 X C Y C に半径 R の水滴が速度 で厚さの小さい静水に衝突する 自由表面の変化に注目しつつ過渡状態を求めよ ただし対称の条件とし 左壁 底面はすべりなし壁とする Prorm: ROP SMC Y C W L LHT W X C