〜前処理法による圧縮性流れの数値計算手法の紹介

Presentation on theme: "〜前処理法による圧縮性流れの数値計算手法の紹介"— Presentation transcript:

1 〜前処理法による圧縮性流れの数値計算手法の紹介
密度に基づく圧縮性流れの数値解法による 非静力学モデルについて 〜前処理法による圧縮性流れの数値計算手法の紹介 神戸大学 地球および惑星大気科学研究室 修士 2年 河合 佑太 2012/06/14 大気セミナー

2 はじめに イントロダクション 本研究について 本日の発表について 圧縮性低速流れに対する数値手法
前処理法について (Weiss and Smith, 1995) 前処理法による密度ベース圧縮性流れの数値解 法による数値計算 GFD シミュレーションへの適用に向けた課題 まとめ

3 イントロダクション

4 本研究について 系の全質量・全エネルギー保存性. 長時間積分を行う際に系の保存量が保存すること が, 重要となる. 理想気体近似をはずす.
木星型惑星大気の数値シミュレーションを想定した堅牢な 完全圧縮系非静力学モデルの開発 非静力学モデルの定式化にあたって次の点を重視したい. 系の全質量・全エネルギー保存性. 長時間積分を行う際に系の保存量が保存すること が, 重要となる. 理想気体近似をはずす. 木星型惑星の数値シミュレーションでは, 実在気体の効果(real gas effect)を真面目に考える必要がでてくるだろう(eg Cloutman, 1984). コンピュータ・クラスタ上での並列計算のための 計算手法の選択と設計.

5 本日の発表について 一方, 非静力学コアで使われるのは, 主に圧力ベー スの手法である.
系の全質量・全エネルギー保存性, 理想気体近似なしを満 たす数値モデル開発のために, 圧縮性低速(低マッハ数)流れ の数値解法を, 気象分野および CFD 分野に渡り調査した. 本研究では, 圧縮性低速流れを数値計算するために, CFD の伝統的な密度ベースの圧縮性流れの数値解法に前処理法を 施した数値手法を採用することにした. 一方, 非静力学コアで使われるのは, 主に圧力ベー スの手法である. 本発表では, 始めに圧縮性低速度流れのための圧力ベース と密度ベースの数値解法を簡単に比較した後, 密度ベースの 数値解法に前処理を施す方法を紹介する.

6 圧縮性低速流れに対する数値解法 1 従来の完全圧縮系の非静力学モデル
CFD 的には非圧縮流れの解法を圧縮性流れに拡張し た圧力ベースの数値解法が取られる. 圧力方程式が解かれる. 密度は診断的に決定. 音波に伴う CFL の制約は音波関連項を陰的に取り扱 う(HE-VI, HI-VI 法)ことで回避.

7 圧縮性低速流れに対する数値解法 2 CFD における伝統的な密度ベース圧縮性流れの数値解法 を大気現象のような低速(低マッハ数)流れに適用できるか? 直接はできない 音波に伴う厳しい CFL 条件のために非常に小さな時間刻みを取らなけれ ばならない('Stiff' な状態). 陰解法を使ったとしても, 収束が遅くまた数値精度に問題が生じる. ただし, 前処理法*を用いることで可能となる. 流速ヤコビアンの固有値を操作することで, 音速の伝播速度を流れの速度 オーダーに置き換えてしまうことで上の問題を解決する. この発表では, Weiss and Smith(1995) の前処理法を紹介する. * Turkel(1987), Choi and Merkle(1993), Weiss and Smith(1995) ほか

8 前処理法による圧縮性(低速)流れ の数値解法 ~ Weiss and Smith(1995)

9 支配方程式 簡単のために２次元 Cartesian 座標系で説明 する. Navier Stokes Eq: （保存系） および 状態方程式
および　状態方程式 保存変数: ｘ 方向の非粘性流束: ｙ 方向の非粘性流束: x 方向の粘性流束: ｙ 方向の粘性流束: (u,v):Cartesian座標における各速度成分, ρ:密度, E:単位質量当たりの全エネルギー P: 圧力, T: 温度, τij: 応力テンソル, κ: 熱伝導係数

10 前処理法：従属変数の変換 Weiss and Smith(1995) では圧力-温度の組み合わせを用 いる.
まず, 依存変数を保存変数からプリミティブ変数に変換す る プリミティブ変数の選択には任意性がある. 圧力-温度, 圧力-エントロピー, 圧力-エンタルピー(Turkel, 1999) Weiss and Smith(1995) では圧力-温度の組み合わせを用 いる.

11 前処理法：従属変数の変換 プリミティブ変数 の選択について 非圧縮流れの数値解法では自然な選択である
プリミティブ変数 の選択について 非圧縮流れの数値解法では自然な選択である 粘性流束内の速度勾配/温度勾配, 非粘性流束内の圧 力勾配を求める際に高次の再構築を用いるが, それを 高精度で行うには, p-T のプリミティブ変数の選択が 有利である. 依存変数として圧力を選ぶと, 系の音波伝播の表現を 表面化することができる.

12 前処理法：従属変数の変換 プリミティブ変数 の選択について 依存変数として圧力を選ぶと, 系の音波伝播の表現を 表面化することができる.
プリミティブ変数 の選択について 依存変数として圧力を選ぶと, 系の音波伝播の表現を 表面化することができる. : 系を保存形式から非保存形式に変換するための変換行列 理想気体の場合, c: 音速

13 前処理法：前処理行列 低マッハ数流れにおいて, 系の各固有値(eg 一次元では u, u±c )のオーダは同じオーダとならず, Stiff な状態を引 き起こす. 前処理法では, 音波伝播に関する項を局所的な速度に置 き換えてしまう. 系の固有値を全て同じオーダする 前処理行列： ここで, Cp:定圧比熱. Ur は参照速度で以下のように定義する. (ε: 十分小さな値) 粘性流の場合は, Ur が拡散速度より小さくならないように制限する.

14 前処理法：前処理行列 保存形式における前処理された系 定常状態において, 方程式は保存特性を満たす.
前処理により方程式系の時間発展は破壊されるので, 非定常計 算には二重時間ステップを用いる. 保存形式における前処理行列: ここで, 全エンタルピー:

15 前処理法：置換された系の固有値について 前処理された系の固有値 Ｘ方向の流速ヤコビアンの固有値
Ur=c のとき α=0となり, 系の固有値は伝統的な形式 u±ｃ をとる. 一方, Ur がゼロに近づくと, α=1/2 となり, 系の全ての固有値は u と同じ オーダになる. => 音速と流速の大きさが離れることによって生じる”Stiff”な問題を解決す ることができる!!

16 前処理法：非定常問題への適用 前処理法を施した方程式を使って, 時間精度の良い解 を得るには, 二重時間ステップを用いる.
擬時間無限大で, もとの方程式に帰着する. ｔ：物理時間, τ: 擬時間(時間発展で用いられる)

17 前処理法：非定常問題への適用 物理時間微分を２次精度後退差分で離散化し, 擬時間微分 を多段解法を使ってゼロとなるまで, 時間発展させる場合
i:多段解法のカウンタ(i=1,..,m), n: 物理時間レベル αi: 多段解法の係数 * ΔQ が十分小さくなったとき, その を持ってして, 次の物理時間レベルの解 を更新する.

18 前処理法：非定常問題への適用 現実的な 3 次元問題に対しては, 擬時間における時間発 展の収束加速が必要である.
陰解法の使用(Δτを大きく取る) 近似因数分解法 LU-SGS 法(Yoon and Jameson, 1988) 陽解法における収束加速 多重格子法, (陽的)残差平均法 局所時間ステップ

19 前処理法による圧縮性流れ数値解法を用い た ２次元低速流れの数値計算

20 数値計算手法 支配方程式系 ２次元圧縮性 Navier-Stokes 方程式 (完全圧縮系) 空間離散化 再構築
(セル中心法による)有限体積法 非粘性流束 Roe の近似リーマン解法 (MUSCL 法を用いて高精度化, 前処理法の適用) 粘性流速 中心差分的評価 時間積分 二重時間ステップ * 物理時間：２次精度後退差分 * 擬時間 * 陽解法: 多段 Runge=Kutta 法 * 陰解法: LU-SGS 法 * 収束加速は未実装(多重格子法を実装予定) 乱流モデル なし 並列化 OpenMP によるスレッド並列

21 再構築 セル境界の各側の物理量を, 近傍の格子点上のセル平 均量データを使って補間する.
再構築によってえられる境界上の依存変数の値を使って, 非粘性数値 流束を評価する. 左右のセル平均値をそのまま用いて, 数値流束を風上化すると一次精 度になる. 高精度化には, MUSCL 法などが用いられる. 解の単調性を維持するために, 制限関数が利用される. 今回は, 再構築に MUSCL 内挿(風上側にバイアスした３次精度)を使 用し, 制限関数は Van Albada の制限関数を用いた. L R : 計算格子点 : コントロールボリューム

22 数値流束の評価 離散化された非粘性流束の評価には, 流束差分離法 (Roe, 1986)を用いる. 前処理法の適用
FR, FL: セル境界の各側の再構築された解 WR, WL を用いて計算された数値流束. また, A は流速ヤコビアン　　　であり, ここで, M は行列 A を対角化する固有行列, L R ここで, MΓ は行列 AΓ を対角化する固有行列,

23 圧縮性流れの数値解法を用いた 数値実験 前処理なしの圧縮性流れ数値解法の実装が, 正しく行われ ているかをテスト. ２次元爆破テスト
計算領域: [-1,-1]×[1,1] 計算格子数: 100×100 初期条件: 境界条件 固体壁：スリップ条件, 断熱壁 ただし,

24 <2次元爆破テストの数値計算の結果>
コンター: 密度 ベクトル: 速度 * なお, 計算不安定を防ぐために, 初期場データの不連続境界に対してスムージングを施している.

25 計算結果の妥当性チェック 密度 圧力 速度 本研究 Xiao et al(2006) Fig 3：
VISAM3 法(CIP 法)による２次元爆破シミュレーションの結果. 実線は, 高精度一次元軸対称モデルによる結果(解析解として利用) 本研究

26 前処理法付き圧縮性流れの数値解法 を用いた数値実験
２次元キャビティ流れ 計算領域: [0,0]×[1,1] 計算格子数: 80×80 初期条件 静止流体. 密度 1, 圧力は 4, 40, 400 と各数値実験で変化させ る. 境界条件 固体壁：滑りなし条件(y=1においてu=1, それ以外の壁速度はゼ ロ), 断熱壁 無次元パラメタ マッハ数: 0.42, 0.13, 0.042 レイノルズ数: 1000

27 各マッハ数における２次元キャビティ流れの 数値計算結果
初期圧力:4 Pa マッハ数: 0.42 初期圧力:40 Pa マッハ数: 0.13 初期圧力:400 Pa マッハ数: 0.042 各数値実験ともレイノルズ数は, 1000 に固定. 上の図は, t=25(準定常状態) における圧力・速度分布.

28 各マッハ数における２次元キャビティ流れの 数値計算結果
X=0.5 における速度 u Y=0.5 における速度 v * 初期圧力: 4*10^4 Pa, マッハ数: * 格子数 80*80 ( 圧縮性流れの Cavity 問題でも, 超低マッハ数にすると, 非圧縮 Cavity 流れの結果と漸近してくることが知られている) + 印は, Ghia(1982) の非圧縮 Cavity 流れの数値計算結果(格子数 128*128).

29 GFD シミュレーションへの適用に向けた 課題

30 前処理法を施した密度に基づく圧縮性流の数値解法による GFD シミュレーションへの適用に向けた課題 1
重力項, コリオリ項を付加する. 今回紹介した数値手法を適用した数値モデルは, 各大気波動(内部重力波, ロスビー波 etc)の特性を 正しく表現できるのか? 計算領域を ３ 次元へ拡張する. 今回の数値計算手法においては, 2 次元モデルに 単に 1 次元付け加えることは(実装上)容易であ る. 3 次元球殻計算領域に対しては, 多面体格子を利 用すればよいだろう. 立方格子(Sadourny, 1972) この場合, 一般曲線座標系に対応する必要があ る. （バネ力学を用いて修正した）正二十面体格子 (Tomita et al, 2001)

31 前処理法を施した密度に基づく圧縮性流の数値解法による GFD シミュレーションへの適用に向けた課題 2
木星型惑星大気の数値シミュレーションを行う ための, 理想気体近似の除去 今回紹介した前処理法は理想気体を前提したも のでは無い. Weiss and Smith 型の前処理法を実在気体に 拡張している例は存在する(例えば, Zong and Yang, 2007).

32 まとめ

33 まとめ 前処理法を適用することで, 大気現象のような 低マッハ数流れの数値シミュレーションに圧縮 性流れの数値解法を用いることができる.
完全圧縮系非静力学モデルへの応用を念頭におきなが ら, 密度に基づく圧縮性流の数値解法を紹介した. 前処理法を適用することで, 大気現象のような 低マッハ数流れの数値シミュレーションに圧縮 性流れの数値解法を用いることができる. 前処理法を施した圧縮性流の数値解法による２次元流 体シミュレーションコードを作成し, それによる数値 計算結果を例に示した. 低マッハ数流れの数値シミュレーションにおい ても, 定量的に正しい数値解が得られた. 前処理法を施した圧縮性流の数値解法の GFD 的数値 シミュレーションへの適用に向けた課題を整理した.

34 参考文献・参考資料 Weiss, M, J, Smith, A, W, 1995: Preconditioning Applied to Variable and Constant Density Flows. AIAA JOURNAL, 33, No 11. Yamamoto, S, 2005: Preconditioning method for condensate fluid and solid coupling problems in general curvilinear coordinates. J. Comput. Phys., 207, 240—250. Xiao, F, et al, 2006: Unified formulation for compressible and imcompressible flows by using multi-integrated moments Ⅱ: Multi-dimensional version for compressible and imcompressible flows. J. Comput. Phys., 213, 31—56. Zong, N and Yang, V, 2007: An efficient preconditioning scheme for real-fluid mixtures using primitive pressure-temperature variables: I. J. Comput. CFD., 21, Li, X and Gu, C:On the Mechanism of Roe-type Schemes for All-Speed Flows 数値流体力学編集委員, 1995:圧縮性流体解析, 東京大学出版会 藤井考藏, 1994:流体力学の数値計算法, 東京大学出版会

35

36