Lock in $30 Savings on PRO—Offer Ends Soon! ⏳

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

matsubaraDaisuke
October 29, 2016

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

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

matsubaraDaisuke

October 29, 2016
Tweet

More Decks by matsubaraDaisuke

Other Decks in Science

Transcript

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

    定常/層流・乱流(RANS)/SIMPLE法 • rhoSimplecFOAM定常/層流・乱流(RANS)/SIMPLEC法 • rhoPimpleFOAM 非定常/層流・乱流/HVAC用 出典 http://dot-prototype.appspot.com/OpenFOAM.html
  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. 密度ベースソルバーとは? 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)等を求める
  5. 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などの未知変数は一意に求まるため、残差はゼロと出 力される(バグではない)。
  6. 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スキームが実装されている。
  7. OPENFOAMに実装されている基本的な 圧縮性流体ソルバーの一覧 • rhoCentralFOAM KURGANOVとTADMORによる中心-風上スキームに基づく密度ベースソルバー • sonicFOAM 非定常/遷音速・超音速/層流・乱流 • rhoSimpleFOAM

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

    fvm::ddt(rho) + fvc::div(phi) ) Phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg; = . (. . ) + / (/ / ) 係数 風上成分と風下成分を係数で結合した形となる。 今回は理論的背景は考えない。 例えばある時、. , / が支配的ならば風上的、等分で あれば中心差分となる。運動量、エネルギーの式も同じ 考え方。 運動量方程式とエネルギー方程式については割愛
  9. 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”
  10. 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; } 補間する場の変数と方向を引数にする
  11. 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
  12. 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
  13. 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が一般的
  14. rhoCentralFoamを使う際に用いられるパラメータ • constant/thermophysicalProperties中 Mixture { Specie { nMoles 1; molWeight

    11640.3; } thermodynamics { Cp 2.5; Hf 0; } transport { mu 0; Pr 1; } } 定圧比熱 エンタルピー 0ならば非粘性(invisid) プラントル数 モル数 モル質量