2018」用記事です.はじめに自分は流体系の人間なので、今回も何か流体、特に数値流体の話がしたいなぁと思って考えてました。なったので、今回は空気力学第四のパクリで圧縮性の話をしたいと思います。圧縮性流体を考えるCFDの世界ではどのように衝撃波を捉えるかが重要になります。CFD界隈では銀本と呼ばれる本(流体力学の数値計算法)を見ると、衝撃波捕獲法として重要な技法として二つあげています。流体を少しでも齧ったことのある方ならきっとオイラー方程式くらいは知ってるはず。今回の問題設定で考えるのは、そのオイラー方程式が使えるような状況、かつ一次元の衝撃波菅問題に関して考えたいと思います。当たり前といえば当たり前ですが、実はこの三つを完全に満たすようにCFDを作るのは基本的にはできないはずで、今回は触れないですが、精度とこれらの保存則は割とトレードオフの関係にあります。ここで今回の問題で以上の保存則を表す方程式は以下の三つになります。流体の熱エネルギーであることに注意しましょう。この式、普段古典的な流体力学の授業で紹介されるような式の形をしていません。全エネルギーを括ったからとかいうのはありますが、何より重要なのは、全て「時間微分項」+「空間微分項」の形で表記されている点です。ゆえに、上の式は以下のように書き換えることが可能です。これが「保存形」と呼ばれる流体のシステム方程式の形になります。確かにこれ、保存則の形そのまんまですしね。「流束」と呼ばれます。ここで、式(\ref{qt})の形のままでは何もできません。この式を局所的に線形化します(係数マトリックスを凍結するといいます)。固有値分解すると次の三つの固有値を得ることができます。同じ形であることがわかります。輸送速度に対応していると予測が立ちます。非対角成分を持ちますが、対角化することで独立な三つの波動方程式に分解でき、その三つの波の特性速度がこの行列の固有値として現れているのだと考えることができます。実際、対角化は以下のように行われます。基本的に「超音速」を超えた場合に風上からしか来ないからです。上流に伝播しなくなるからです。このような背景から今回は風上差分法を使います。ここまでをまとめるとこんな感じです。情報は風上側から流れてくるので、風上差分法を用いるのが適切さて、今までの話から、システム方程式の解法では三つの波が独立に存在し、それぞれに対して高解像度差分法を適用することが可能になります。CFDの内部では次のグラフのような形でデータが存在しているとイメージしてください。CFDのイメージこの次の時間ステップを求めるためには各セル境界面での流束を評価する必要があります。このように、二つの領域に渡って不連続な初期値をもつ双曲型の連立方程式を「リーマン問題」と呼びます。リーマン問題このリーマン問題の解から、セル境界面における通過する流束を求めることができる。この式から、各セルにおける物理量分布を求めることが可能になります。計算コストの高いリーマン問題を厳密に解く必要はありません。その近似手法の中でもよく知られた(古典的な)ものは以下の二つになります。流束差分法(FDS)流束分離法(FVS)以下ではそれぞれの手法の説明をした後に、結果を比較したいと思います。Splitting)いわゆる「Roeスキーム」と呼ばれる手法になります。数値流束を評価することで風上差分を行う、ということです。具体的に見てもらった方がわかりやすいかもしれないので、早速説明していきます。Roeは非線形スカラー方程式に対する1次の風上差分を非線形のシステム方程式に拡張しました。このスキームでは1次精度の数値流束は以下の式で評価されます。表される量は以下のRoeの平均を用いて評価されます。意味は全ての行列要素で絶対値をとる、という意味です。非線形関数として定義されています。求めた後に式(\ref{int})に基づいて時間積分を行えば良いです。数値粘性を付加しています。逆にいうと、「跳び」の少ない衝撃波から離れた領域に対しては数値粘性を付加しないので、正確になります。この点がFDSの優れている点です。いう意味でなぜか近似リーマン解法に分類されています。この手法の特徴としては、「物理的な流束を各特性波の方向に応じて差分方向を切り替える」という点です。考え方はとても簡単です。右からくる負の固有値に対応する流束です。正の固有値を持つということは波(情報)が左からくるということなので後退差分を利用します。逆に負の固有値を持つならば波(情報)が右からくるということなので前進差分を利用します。それでは具体的に手法を見ましょう。流速の分解は実際には以下の式で行います。負の固有値を持つもののみからなる。今回は一次精度の流束分離を行い、その時には次のように置くことができる。それぞれ左からくる波、右からくる波に対応しています。それぞれ使用する必要がある。次のように評価されます。求めた後に式(\ref{int})に基づいて時間積分を行うだけです。ここでFDSとFVSの違いに関して軽く述べると、式(\ref{ffds})と(\ref{flux_fvs})を見比べると、最後の項だけ違うことがわかります。二次精度中心差分を形成することがわかります。残った最後の項が拡散項として働き、このスキームに風上差分のような特徴を加えていると考えることができます。安定だが鈍るのが特徴で、その両者の組み合わせをするのが差分スキームの基本です)この拡散項がFDSでは必ず正になり、拡散を引き起こしますが、一方FVSではここが必ずしも正にはならないので、そこが本質的に異なっていると考えることができます。実装ここに示すと醜くなるので、Githubに置いておきました。見てみたい方はみてみるといいかもしれません。適当にパパッと作ったので、クソコードなのは許してください。あと、Eigen使ってるので、そこだけ注意お願いします。github.com結果さて、それではFDSとFVSが実際の問題に適用した時にどのような挙動になるのか見てみましょう。今回は問題として、一次元の衝撃波管問題を考えました。して、頻繁にスキームのテスト問題として用いられます。圧力グラフさて、FDSはやはり特徴として衝撃波付近で鈍っていますね。これはそもそもスキームの段階で数値粘性がそのように働くよう作ったので狙い通りです。うーん...FDSはまだそんなものかと思えますが、FVSの結果が異常に悪いですね...もしかしたら実装を間違えたかもしれないので、そこだけ申し訳ありません。一応、FVSは論文に衝撃波管問題への適用結果があったので、載せておきます。今回私たちが学んだのは本当に最初のFVSです。このあと、FVSというジャンルの中には様々に新しいスキームが誕生しました。最初のFVSは使い物にならなかったからで、もしかしたら、自分の実装は間違ってなかったかもしれません(笑)終わりにどうですか?意外と圧縮性のCFDってそんなに難しくないんですよね。考え方はとても簡単なので皆さんに理解してもらえたなら幸いです。...答えはFVSなんです!なんででしょうね?今圧縮性で頻繁に用いられているのはFVSから派生したAUSM系スキームと呼ばれるものです。最初はあんなに結果が悪くても、のちに誰かが発展させて広まっていくのは素晴らしいなぁと思いました。(2)はてなブログをはじめよう!kane-pleaseさんは、はてなブログを使っています。ブログを報告するスターをつけました引用をストックしましたストック一覧を見る閉じる引用するにはまずログインしてくださいログイン閉じる引用をストックできませんでした。再度お試しください閉じる限定公開記事のため引用できません。