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

Slides:



Advertisements
Similar presentations
Division of Process Control & Process Systems Engineering Department of Chemical Engineering, Kyoto University
Advertisements

相対論的場の理論における 散逸モードの微視的同定 斎藤陽平( KEK ) 共同研究者:藤井宏次、板倉数記、森松治.
宇宙ジェット形成シミュレー ションの 可視化 宇宙物理学研究室 木村佳史 03S2015Z. 発表の流れ 1. 本研究の概要・目的・動機 2. モデルの仮定・設定と基礎方程式 3. シンクロトロン放射 1. 放射係数 2. 吸収係数 4. 輻射輸送方程式 5. 結果 6. まとめと今後の発展.
Determining Optical Flow. はじめに オプティカルフローとは画像内の明る さのパターンの動きの見かけの速さの 分布 オプティカルフローは物体の動きの よって変化するため、オプティカルフ ローより速度に関する情報を得ること ができる.
Computational Fluid Dynamics(CFD) 岡永 博夫
有限差分法による 時間発展問題の解法の基礎
Fill-in LevelつきIC分解による 前処理について
概要 基礎理論 1.応力とひずみおよび平衡方程式 2.降伏条件式 3.構成式(応力-ひずみ関係式)
自己重力多体系の 1次元シミュレーション 物理学科4年 宇宙物理学研究室  丸山典宏.
医薬品素材学 I 3 熱力学 3-1 エネルギー 3-2 熱化学 3-3 エントロピー 3-4 ギブズエネルギー 平成28年5月13日.
〜 「バネ力学を用いた正二十面体測地線格 子の改良(Tomita et al, 2001)」
ウェーブレットによる 信号処理と画像処理 宮崎大輔 2004年11月24日(水) PBVセミナー.
(Fri) Astrophysics Laboratory MATSUO Kei
Finger patternのブロック化による 陰的wavelet近似逆行列前処理の 高速化
スペクトル法による数値計算の原理 -一次元線形・非線形移流問題の場合-
高精度有限体積法による 非静力学惑星大気循環モデルの開発 神戸大学 地球および惑星大気科学研究室 河合 佑太
大阪工業大学 情報科学部 情報システム学科 宇宙物理研究室 B 木村悠哉
相対論的輻射流体力学における 速度依存変動エディントン因子 Velocity-Dependent Eddington Factor in Relativistic Photohydrodynamics 福江 純@大阪教育大学.
北海道大学大学院理学研究科地球惑星科学専攻 地球流体力学研究室 M1 山田 由貴子
反応性流体力学特論  -燃焼流れの力学- 燃焼の流体力学 4/22,13 燃焼の熱力学 5/13.
流体のラグランジアンカオスとカオス混合 1.ラグランジアンカオス 定常流や時間周期流のような層流の下での流体の微小部分のカオス的運動
計算アルゴリズム 計算理工学専攻 張研究室 山本有作.
文献名 “Performance Tuning of a CFD Code on the Earth Simulator”
サポートベクターマシン によるパターン認識
応用数学 計算理工学専攻 杉原研究室 山本有作.
速度勾配依存 変動エディントン因子 Velocity-Gradient-Dependent Relativistic Variable Eddington Factor Plane-Parallel Case 福江 純@大阪教育大学.
圧力発展格子ボルツマン法による大規模気液二相流GPUコードの開発 ならびに多孔体浸潤液滴シミュレーション
北大MMCセミナー 第38回 Date: 2015年2月13日(金)16:30~18:00 Speaker: 宮路 智行(明治大学)
ガンマ線バーストジェット内部における輻射輸送計算
非エルミート 量子力学と局在現象 羽田野 直道 D.R. Nelson (Harvard)
ひび割れ面の摩擦接触を考慮した損傷モデル
正規分布における ベーテ近似の解析解と数値解 東京工業大学総合理工学研究科 知能システム科学専攻 渡辺研究室    西山 悠, 渡辺澄夫.
スペクトル法の一部の基礎の初歩への はじめの一歩
電磁流体力学乱流の高精度・高並列LESシミュレーションコード開発研究
半無限領域のスペクトル法による竜巻を模した渦の数値実験に向けた研究開発
「移流輸送の 高次精度かつ高解像度スキームの開発」
数値計算ゼミの進行に 関する提案 ~実りある春にするために           新しい形のゼミを求めて~ 清水慎吾、野村光春.
Jh NAHI 横田 理央 (東京工業大学) Hierarchical low-rank approximation methods on distributed memory and GPUs 背景  H行列、H2行列、HSS行列などの階層的低ランク近似法はO(N2)の要素を持つ密行列をO(N)の要素を持つ行列に圧縮することができる。圧縮された行列を用いることで、行列積、LU分解、固有値計算をO(NlogN)で行うことができるため、従来密行列の解法が用いられてきた分野では階層的低ランク近似法
バルク法について ~deepconv を用いて地球の積雲対流の数値計算をするにあたって~
燃焼の流体力学 4/22 燃焼の熱力学 5/13 燃焼流れの数値解析 5/22
独立成分分析 (ICA:Independent Component Analysis )
システム制御基礎論 システム工学科2年後期.
化学工学基礎 −後半の後半− 第1回目講義 (2009年7月10日) 1 担当 二又裕之 物質工学1号館別館253ー3号室
開放端磁場における低温プラズマジェットに関する研究
連続体とは 連続体(continuum) 密度*が連続関数として定義できる場合
FEM勉強会 (第3回).
2009年7月9日 熱流体力学 第13回 担当教員: 北川輝彦.
課題 1 P. 188.
Chapter 26 Steady-State Molecular Diffusion
円柱座標系の基底関数系を用いたSCF法による 円盤銀河のシミュレーション
第1回、平成22年6月30日 ー FEM解析のための連続体力学入門 - 応力とひずみ 解説者:園田 恵一郎.
4. システムの安定性.
「ICAによる顔画像特徴量抽出とSVMを用いた表情認識」
竜巻状渦を伴う準定常的なスーパーセルの再現に成功
これらの原稿は、原子物理学の講義を受講している
潮流によって形成される海底境界層の不安定とその混合効果
対象:せん断補強筋があるRCはり(約75万要素)
ポッツスピン型隠れ変数による画像領域分割
電気回路学Ⅱ 通信工学コース 5セメ 山田 博仁.
・Bernoulli(ベルヌーイ)の定理
卒論中間発表 2001/12/21 赤道の波動力学の基礎 北海道大学理学部 地球科学科 4年 山田 由貴子.
力覚インタラクションのための 物理ベースモデリング
北大MMCセミナー 第65回 附属社会創造数学センター主催 Date: 2017年4月20日(木) 16:30~18:00
How shall we do “Numerical Simulation”?
確率的フィルタリングを用いた アンサンブル学習の統計力学 三好 誠司 岡田 真人 神 戸 高 専 東 大, 理 研
背景 粒子法(SPH・MPSなど)は大規模流体シミュレーションなどで幅広く利用.一方で,手法の数学的正当化(数値解析)が不十分
? リー・ヤンの零点 これまでの格子QCD計算の結果 今年度の計画 リー・ヤンの零点分布から探る有限密度QCDにおける相構造の研究
固体→液体 液体→固体 ヒント P131  クラペイロンの式 左辺の微分式を有限値で近似すると?
Presentation transcript:

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

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

イントロダクション

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

まとめ

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

参考文献・参考資料 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, 217--230. Li, X and Gu, C:On the Mechanism of Roe-type Schemes for All-Speed Flows 数値流体力学編集委員, 1995:圧縮性流体解析, 東京大学出版会 藤井考藏, 1994:流体力学の数値計算法, 東京大学出版会