Upgrade to Pro — share decks privately, control downloads, hide ads and more …

OpenFOAMの圧縮性流体􏰀ソルバー􏰀解説

 OpenFOAMの圧縮性流体􏰀ソルバー􏰀解説

訂正)
4ページ目に全ての解法が陽解法かの様に説明していますが,OpenFOAMのこのソルバーが陽解法という意味で記載しています.
一般的に密度ベースの解法には陽解法と陰解法の両方を適用することができます.
(2019.6.23)

849be8a73074bccd6f3d5ff3fbbb940f?s=128

matsubaraDaisuke

October 29, 2016
Tweet

More Decks by matsubaraDaisuke

Other Decks in Science

Transcript

  1. 圧縮性流体のソルバーの解説1 松原 大輔 OPENFOAM 2.5 5 7.5 10 12 p

    0.0339 12.7 第52回OpenCAE勉強会@関西
  2. OPENFOAMに実装されている基本的な 圧縮性流体ソルバーの一覧 • rhoCentralFOAM KURGANOVとTADMORによる中心-風上スキームに基づく密度ベースソルバー • sonicFOAM 非定常/遷音速・超音速/層流・乱流 • rhoSimpleFOAM

    定常/層流・乱流(RANS)/SIMPLE法 • rhoSimplecFOAM定常/層流・乱流(RANS)/SIMPLEC法 • rhoPimpleFOAM 非定常/層流・乱流/HVAC用 出典 http://dot-prototype.appspot.com/OpenFOAM.html
  3. OPENFOAMに実装されている基本的な 圧縮性流体ソルバーの一覧 • rhoCentralFOAM KURGANOVとTADMORによる中心-風上スキームに基づく密度ベースソルバー • sonicFOAM 非定常/遷音速・超音速/層流・乱流 • rhoSimpleFOAM

    定常/層流・乱流(RANS)/SIMPLE法 • rhoSimplecFOAM定常/層流・乱流(RANS)/SIMPLEC法 • rhoPimpleFOAM 非定常/層流・乱流/HVAC用 出典 http://dot-prototype.appspot.com/OpenFOAM.html
  4. OPENFOAMに実装されている基本的な 圧縮性流体ソルバーの一覧 • rhoCentralFOAM KURGANOVとTADMORによる中心-風上スキームに基づく密度ベースソルバー • sonicFOAM 非定常/遷音速・超音速/層流・乱流 • rhoSimpleFOAM

    定常/層流・乱流(RANS)/SIMPLE法 • rhoSimplecFOAM定常/層流・乱流(RANS)/SIMPLEC法 • rhoPimpleFOAM 非定常/層流・乱流/HVAC用 出典 http://dot-prototype.appspot.com/OpenFOAM.html
  5. 密度ベースソルバーとは? • 密度ベースソルバーは、物体の周りの遷音速や高音速流れに用いられる方法。 • 圧力ベースソルバーは、非圧縮流体や低マッハ数の流れに対して用いられる方法。 密度ベースソルバーVS圧力ベースソルバー 時間差分の部分で誤差が蓄積されていくので、非定常計算の場合、慎重に スキームを選ぶ必要がある。 質量保存則、運動量保存則、エネルギー保存則から保存量を陽的に求め、 保存量から密度や圧力などのプリミティブな変数を求める。

    陰解法を含めることができるため、密度ベースと比較してΔTを大きく取ることができる。 密度ベースに比べて流体シミュレーションでは主流の方法である。 緩和係数や連立方程式など設定するパラメータが多く、経験も必要
  6. 密度ベースソルバーとは? 2.質量保存則から次の時間ステップの密度を求める 3.運動量保存則から次の時間ステップの運動量を求める 5.エネルギー保存則から次の時間ステップの全エネルギーを求める 4.ここで2で求めた運動量を1の密度で割ることで、次のステップの速度が求まる 6.ここで4で求めたエネルギーを1の密度と3の速度を用いて、次のステップの内 部エネルギーや温度、圧力などが求まる // --- Solve

    density solve(fvm::ddt(rho) + fvc::div(phi)); // --- Solve momentum solve(fvm::ddt(rhoU) + fvc::div(phiUp)); U.ref() = rhoU() /rho(); // --- Solve energy solve(fvm::ddt(rhoE) + fvc::div(phiEp) - fvc::div(sigmaDotU) ); e = rhoE/rho - 0.5*magSqr(U); e.correctBoundaryConditions(); thermo.correct(); p.ref() = rho() /psi(); rhoCentralFoam.C 1.補間を用いてセル界面の物理量(phi,phiUp,phiEp)等を求める
  7. rhoCentralFoamを使う際に用いられるパラメータ • system/fvSolution中 Solvers { ”(rho|rhoU|rhoE)” { solver diagonal; }

    } 右辺(緑側)は陽的に求めるため未知変数に対して緩和法 などを使って連立方程式を解く必要がない。 // --- Solve density solve(fvm::ddt(rho) + fvc::div(phi)); // --- Solve momentum solve(fvm::ddt(rhoU) + fvc::div(phiUp)); // --- Solve energy solve(fvm::ddt(rhoE) + fvc::div(phiEp) - fvc::div(sigmaDotU) ); 既知の変数(密度など)で構成 (fvc::) 未知 (fvm::) 粘性項は今回は割愛 粘性を考慮する場合、ラプラシアン項は隠的に与えるためソルバーを設定する 必要がある。 そしてrhoEなどの未知変数は一意に求まるため、残差はゼロと出 力される(バグではない)。
  8. rhoCentralFoam • KURGANOVとTADMORによる中心-風上スキームに基づく密度ベースソルバー // --- Solve density solve(fvm::ddt(rho) + fvc::div(phi));

    // --- Solve momentum solve(fvm::ddt(rhoU) + fvc::div(phiUp)); // --- Solve energy solve(fvm::ddt(rhoE) + fvc::div(phiEp) - fvc::div(sigmaDotU) ); 既知の変数(密度など)で構成 (fvc::) phi、phiUp、phiEpの決め方ごとに様々なスキームが考案されている。 Roeスキーム、HLLEスキーム、HLLCスキーム..etc rhoCentralFoamは、 KurganovとTadmorスキームが実装されている。
  9. OPENFOAMに実装されている基本的な 圧縮性流体ソルバーの一覧 • rhoCentralFOAM KURGANOVとTADMORによる中心-風上スキームに基づく密度ベースソルバー • sonicFOAM 非定常/遷音速・超音速/層流・乱流 • rhoSimpleFOAM

    定常/層流・乱流(RANS)/SIMPLE法 • rhoSimplecFOAM定常/層流・乱流(RANS)/SIMPLEC法 • rhoPimpleFOAM 非定常/層流・乱流/HVAC用 出典 http://dot-prototype.appspot.com/OpenFOAM.html
  10. • 風上と中心スキームをフラックス上で切り替える rhoCentralFoamのソースコード解読 t + ' () • 支配方程式(例:質量保存則) Solve(

    fvm::ddt(rho) + fvc::div(phi) ) Phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg; = . (. . ) + / (/ / ) 係数 風上成分と風下成分を係数で結合した形となる。 今回は理論的背景は考えない。 例えばある時、. , / が支配的ならば風上的、等分で あれば中心差分となる。運動量、エネルギーの式も同じ 考え方。 運動量方程式とエネルギー方程式については割愛
  11. interpolate(rho,pos) = . (. . ) + / (/ /

    ) この時、± 等はどのように決めるのか? rhoCentralFoam.C // --- Directed interpolation of primitive fields onto faces surfaceScalarField rho_pos(interpolate(rho, pos)); surfaceScalarField rho_neg(interpolate(rho, neg)); createField.H surfaceScalarField pos ( IOobject ( “pos”, runTime.timeName(), mesh ), mesh, dimensionedScalar(“pos”, dimless, 1.0) ); negなら”-1.0”
  12. interpolate(rho,pos) rho_pos(interpolate(rho, pos))を探索 directionInterpolate.H //- Interpolate field vf according to

    direction dir template<class Type> tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> interpolate ( const GeometricField<Type, fvPatchField, volMesh>& vf, const surfaceScalarField& dir, const word& reconFieldName = word::null ) { tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tsf ( fvc::interpolate ( vf, dir, "reconstruct(" + (reconFieldName != word::null ? reconFieldName : vf.name()) + ')' ) ); GeometricField<Type, fvsPatchField, surfaceMesh>& sf = tsf.ref(); sf.rename(vf.name() + '_' + dir.name()); return tsf; } 補間する場の変数と方向を引数にする
  13. interpolate(rho,pos) Interpolate(GF,ssf,word )を探索 interpolate( const GeometricField<Type, fvPatchField volMesh>&, const surfaceScalarField&

    , const word& ) これらが引数のinterpolate( ) を探すと //- Interpolate field onto faces using scheme given by name in fvSchemes template<class Type> static tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > interpolate ( const GeometricField<Type, fvPatchField, volMesh>& tvf, const tmp<surfaceScalarField>& faceFlux, const word& name ); surfaceInterpolate.H
  14. interpolate(rho,pos) Interpolate(GF,ssf,word )を探索 //- Interpolate field onto faces using scheme

    given by name in fvSchemes template<class Type> static tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > interpolate ( const GeometricField<Type, fvPatchField, volMesh>& tvf, const tmp<surfaceScalarField>& faceFlux, const word& name ); surfaceInterpolate.C
  15. interpolate(rho,pos) Interpolate(GF,ssf,word )を探索 セル界面 . / 補間の傾き 傾きの選び方は複数用意されて いる =

    . (. . ) + / (/ / ) 他の物理量やパラメータも同じ 考え方
  16. rhoCentralFoamを使う際に用いられるパラメータ • system/fvSchemes中 fluxScheme Kuraganov; interpolationSchemes { default linear; reconstruct(rho)

    vanLeer; reconstruct(U) vanLeerV; reconstruct(T) vanLeer; } KurganovかTadmorを選択 傾きの選び方を選択 vanLeer superBee vanAlbada等 ベクトルの場合は末尾にVを追記すると成分の中で最も 制限の厳しい値をベクトル全体に採用する。 →単調性が高くなり1次の風上的になる vanLeerが一般的
  17. rhoCentralFoamを使う際に用いられるパラメータ • constant/thermophysicalProperties中 Mixture { Specie { nMoles 1; molWeight

    11640.3; } thermodynamics { Cp 2.5; Hf 0; } transport { mu 0; Pr 1; } } 定圧比熱 エンタルピー 0ならば非粘性(invisid) プラントル数 モル数 モル質量
  18. まとめ • rhoCentralFoam(密度ベースソルバー)の概要を説明した。 • 今後やりたいこと 1.補間法の違いによる影響を調べる 2.sonicFoam等の圧力ベースのソルバーとの比較 3. Kurganov Tadmor以外の密度ベースソルバーを実装、比較していきたい(HLLE等)

    問い合わせ先 matsubara.daisuke.52z@kyoto-u.jp