NEE 研究会第 18 回講演討論会 OpenFOAM への計算機能追加連続的データ同化法 (VCA 法 ) の実装 大阪大学大学院工学研究科博士後期課程松尾智仁 内容 1.OpenFOAM を使う理由 1.1 OpenFOAMの特徴 1.2 OpenFOAMを使うにあたって 2.OpenFOAM

Transcription

1 NEE 研究会第 18 回講演討論会 OpenFOAM への計算機能追加連続的データ同化法 (VCA 法 ) の実装 大阪大学大学院工学研究科博士後期課程松尾智仁 内容 1.1 OpenFOAMの特徴 1.2 OpenFOAMを使うにあたって 2.OpenFOAM への計算機能追加 2.1 計算機能の追加の方法 VCA 法とは 計算例 於大阪大学中之島センター OpenFOAM の特徴 オープンソースである => 内部でどのような計算が行われているかわかる => 式や処理を自由に改変することができる 研究で用いる場合 どのような計算条件 ( 離散化の方法 乱流モデル 差分スキーム ) で計算が行われているのかを明らかにする必要がある 従来とは違う計算手法を試したくなることがある 1.1 OpenFOAM の特徴 機能が多彩である => たいていの計算であれば コードを自作しなくても実行できる => たとえば差分スキームの違いによる結果の差異などが容易に評価できる 研究で用いる場合 適切なスキームは条件によって異なる コーディングは主な目的ではない 3 4

2 1.2 OpenFOAM を使うにあたって 導入のハードルが高い 日本語文献がまだ多くはない TutorialやUserGuideがあまり親切ではない Linux 上で動作する GUIが無い 1.2 OpenFOAM を使うにあたって 参考サイトなど OpenCAE 学会 ( OpenFOAM のユーザガイド翻訳など PENGUINITIS( AM/index.html) 計算実行 ソース解読などについての多くの情報 OpenFOAM wiki ( 英語 ) ( OpenFOAM に関する情報多数 ( 非公式 ) 5 Google( 大抵のことは検索すれば見つけられます OpenFOAM を使うにあたって 使用環境 仮想環境 : VMware player (Windows 上で動作 ) 作業環境 : ワークステーション上の Linux(ssh 接続 ) Linux : Ubuntu bit OpenFOAM : OpenFOAM ソースコード操作 :Eclipse ( with Pleiades ) 結果処理 : paraview 2.OpenFOAM への計算機能追加 ソルバの場所 : /opt/openfoam2.3.0/applications/solver 以下 ソースコードは C++ で書かれている いわゆるオブジェクト指向で クラスの継承 オーバーロード ( 関数 / 演算子の多重定義 ) オーバーライド ( 関数の上書き ) といった C++ に特有の機能が使われている 高度にモジュール化されている 機能追加は比較的簡単に行なえる 7 8

3 2.OpenFOAM への計算機能追加 クラスや関数の定義はすでになされているため 主としてそれらを組み合わせることで任意のソルバを作成する たいていは main 関数だけを作り変えれば足りる 例 :simplefoam の main 関数 simplefoam : 非圧縮流体の定常流れを計算するソルバ int main(int argc, char *argv[]) { #include "setrootcase.h" #include "createtime.h" #include "createmesh.h" #include "createfields.h" #include "createfvoptions.h" #include "initcontinuityerrs.h" 計算ケースの読み込み 計算条件 ( 時間 ) の読み込み メッシュの作成 変数 ( 速度 圧力など ) の作成 オプションの設定 初期残差の計算 9 simplecontrol simple(mesh); SIMPLE の繰返しの制御の作成 例 :simplefoam の main 関数 ( つづき ) Info<< " nstarting time loop n" << endl; while (simple.loop()) { Info<< "Time = " << runtime.timename() << nl << endl; 方程式を追加する場合 微分方程式は行列の形で表現される } { } #include "UEqn.H" #include "peqn.h" turbulence->correct(); runtime.write(); 速度の計算圧力の計算 乱流に関する計算結果ファイル書き出し Info<< "ExecutionTime = " << runtime.elapsedcputime() << " s" << " ClockTime = " << runtime.elapsedclocktime() << " s" << nl << endl; } Info<< "End n" << endl; return 0; SIMPLE のループ 例 : 非圧縮定常流れの運動量方程式 ( ν U ) = p ( UU ) + eff fvscalarmatrix UEqn( fvm::div(phi, U) + turbulence->divdevreff( U) == fvoptions(u) ); solve( UEqn() == -fvc::grad(p) ); 対流項 粘性項 オプション ( 何もしない ) 圧力項 12

4 方程式を追加する場合 微分方程式は行列の形で表現される 例 : 非定常の物質の保存式 C + ( ρuc) = ( Γ C) + S t fvscalarmatrix Equation( fvm::ddt(conc) +fvm::div(phi,conc) == fvm::laplacian(turbulence->nueff(), conc) +fvoption(conc) ); 非定常項対流項 拡散項オプション ( ソース項 ) 13 VCA 法とは Variational Continuous Assimilation method データ同化法の一種で 観測データを用いて CFD モデルを修正しながら計算を行う手法 Derber(1989) A variational continuous assimilation technique 数値解析 観測データ データ同化 修正された数値解析 14 VCA 法の概要 - CFD モデルの修正式を定義 n A n n 1 Ψ = Ψ + λ n φ - 評価関数を用いて観測値との誤差を評価 1 P p ~ p T p ~ p I = 2 ( Ψ Ψ ) ( Ψ Ψ ) p= 1 - 誤差を最小化する修正項を求める P p I Ψ = φ p= 1 φ T p ~ p ( Ψ Ψ ) = 0 CFD の支配方程式 修正項 計算値と観測値の差分 I φ 修正項 φ が評価関数 I を最小化するとき I の φ に関する勾配が 0 になる VCA 法の計算フロー 計算開始 CFD 解析 観測値の読み込み CFD と観測値の誤差を評価 計算終了 修正項を更新

5 計算例 : 未知の汚染放出源の放出源位置と濃度の推定 室内に未知の放出源があるケースを想定 吸込み口の濃度を観測データとして与える 流れ場 / 放出時刻は既知 観測データ取得位置における濃度のデータから 放出源の位置と放出強度を推定 観測データ取得位置 U [m/s] CFD 計算条件 計算条件 流れ場 非圧縮等温定常流れ 吸込口 1.0 m/s 吹出し口 自由開放 放出源大きさ 0.01 m2 ( 室内中央 ) 放出強度 g/m2/s 放出時間 1.0 秒間 速度 - 圧力解法 SIMPLE 時間の差分スキーム 1 次精度オイラー法 (Euler) 対流の差分スキーム 2 次精度中心差分 +TVD 制限 (Gauss limitedlinear) 勾配の差分スキーム 線形補間 (Gauss Linear) 拡散の差分スキーム 線形補間 (Gauss linear corrected) 0.0 濃度分布計算結果 VCA 計算条件 計算条件 流れ場 既知 放出源 未知 ( 放出源がない場合の数値計算からスタート ) 観測点数 10 点 観測間隔 0.1 sec 観測点位置 吸込口 U [m/s] 2.0 観測データ取得位置

6 濃度分布 正解値との比較 正解値 推定値 修正項 φ の分布 φ [g/ms/s] VCA 計算条件 濃度 放出源位置ともに 実際よりも拡散している 観測される濃度が拡散していることを反映 制約条件として 放出源強度は負にならない という条件を与えてみる

7 正解値との比較 修正項 φ の分布 φ [g/ms/s] 正解値 推定値 制約条件がない場合に比べ 推定精度が向上? 適切な制約条件を与えることで 推定精度を高めることができる? まとめ 前半では オープンソース CFD ライブラリである OpenFOAM を用いて数値解析をおこなう理由と方法について簡単に触れました 後半では 実際に OpenFOAM に計算機能を追加して行った計算を紹介しました VCA 法を用いたデータ同化計算機能を追加しました 計算例では 未知の物質放出源の位置 / 強度 および濃度分布の推定をおこなった例を紹介しました 28