2021-06-06

OpenFOAM【chtMultiRegionFoamで輻射：viewFactorモデルについて】

OpenFOAM FVM CAE

OpenFOAMでchtMultiRegionFoam+viewFactorの導入について検討したのでそのまとめです。

viewFactorモデルの考え方

形態係数に基づいて輻射を計算します。OpenFOAMではfluid側で新しく輻射の熱流束に関する変数qrを追加して計算します。qrの境界条件は、fluid側の境界条件および固体流体の界面であるfluid_to_solidで設定されます。qrはfluid_to_solid、solid_to_fluidの両設定により固体側へ伝えることができます。

viewFactorの実際の導入についての方針

tutorialにはchtMultiRegionFoam(chtMultiRegionSimpleFoamではない)+viewFactorが存在していないので、viewFactorモデルを実際に導入するには

chtMultiRegionFoamで正常に計算できるモデル・条件を設定
viewFactorを追加

するのが正着と思われます。

viewFactorを追加するときに設定する部分

/0.org にqrを追加する
/constant中のファイルは
- /solidにradiationPropertiesを追加する。モデルはopaqueSolid。
- /fluidにradiationPropertiesを追加する。モデルはviewFactor。
- /fluidにboundaryRadiationProperties、viewFactorsDictを追加する。
/system中のファイルは
- changeDictionaryDictを修正。
- /solidのT/solid_to_fluidのtypeをturbulentTemperatureRadCoupledMixedに修正。
- /solidのT/solid_to_fluidにqr, qrNbrを追加。
- /fluidのT/fluid_to_solidのtypeをturbulentTemperatureRadCoupledMixedに修正。
- /fluidのT/fluid_to_solidにqr, qrNbrを追加。
- /fluidにqrを追加
- /fluidの boundaryにinGroupsでwall, viewFactorWallを追加

なお、前処理（メッシュ作成）で

$ faceAgglomerate
$ viewFactorsGen

を実行する必要があります。

デバッグの覚書（トラブルシューティング）

輻射を受ける対象物（固体）の温度が上がり続ける（下がり続ける）

平衡状態に達するハズなので、温度が上がり続けることはないです。直接的には、

壁面の輻射の境界条件が初期条件から更新されていない。

ことが考えられます。paraviewなどで、qrの時間経過を確認しましょう。時間の経過とともに平衡へと向かう熱流束に収束していくはずです。qrはfluid側の境界面およびfluid_to_solidで可視化できます。

見直すべき部分としては

境界条件、物性の設定
fluid_to_solidの伝熱物性の見直し

qrの計算はfluid側で実施されているので、solid側の温度情報がfluid側に反映されていないとqrの更新はされないと推察されます。

温度が振動して、最後は異常終了で計算停止

標準出力で出てくるMin/max Tあたりが大きな－の値として出力されたりします。

qrの値、輻射境界面のfluid温度が振動

したりしてます（paraviewとかで確認）。

輻射での熱流束は平衡に達するよう進行しますが、各タイムステップの1回目においては平衡前の温度を参考に熱流束が計算され、過剰に評価してしまいそのまま温度が暴走することがあります。デバッグとしては、まず輻射のモデルをoffにした条件で適切に計算が実施されることを確認したうえで

時間刻みを短くする
qr, hを緩和する
輻射面を細かくする

が対処ほうで考えられます。緩和はhよりもqrの方が有効な気がしています（体感）。qrの緩和はfluidのfvSolution中でrelaxationFactors -> fieldsで設定できるようです。

relaxationFactors
{
    fields
    {
        rho             1.0;
        p_rgh           0.7;
        qr              0.3;
    }
    equations
    {
        U               0.3;
        h               0.5;
        "(k|epsilon|omega)" 0.7;
        ".*Final"       1;
    }
}

輻射面を細かくするにはviewFactorsDictのnFacesInCoarsestLevel（大きくする）、featureAngle（大きくする）で設定します。

参考にしました

OpenFOAM【chtMultiRegionFoam のpostProcessでsingleGraph】

OpenFOAM

OpenFOAMの結果の可視化・グラフ化についてす。

chtMultiRegionFoamで実施したときにうまくいかなかったので。

$ postProcess -func singleGraph

通常コマンドで実行すると、こんな感じのエラーが出まくります。

Time = 15000
--> FOAM Warning : 
    From function virtual bool Foam::functionObjects::fieldSelection::checkSelection()
    in file functionObjects/fieldSelections/fieldSelection/fieldSelection.C at line 166
    Field T not found

Reading fields:

Executing functionObjects
--> FOAM Warning : 
    From function Foam::label Foam::sampledSets::classifyFields()
    in file sampledSet/sampledSets/sampledSetsGrouping.C at line 83
    
Cannot find registered field matching 1(T)
No fields to sample

複数領域/マルチリージョンでのポスト処理をするときには、以下のように領域指定が必要なようです。

$ postProcess -func singleGraph -region solid

解析結果のフォルダ分けが領域ごとになされているのが理由と推察されます。

参考

matsuo-san.hatenablog.com

2021-04-10

OpenFOAM【物性の温度依存性の設定２】

OpenFOAM

自作の物性の追加します。

前準備（方針の確認）

まずはリンクの構成を理解するところから。

$FOAM_SRC/thermophysicalModels/solidSpecie/transport

ここに物性。一例として、exponentialSolidTransport.Hのリンク先を確認。

OpenFOAM: API Guide: src/thermophysicalModels/solidSpecie/transport/exponential/exponentialSolidTransport.H File Reference

$FOAM_SRC/thermophysicalModels/solidThermo/solidThermo/solidThermos.C

ここが物性を呼んでいることを確認します。solidThermos.Cを読むと

makeSolidThermo
(
    solidThermo,
    heSolidThermo,
    pureMixture,
    exponentialSolidTransport,
    sensibleEnthalpy,
    hPowerThermo,
    rhoConst,
    specie
);

と書いてあるので、ここに自作の物性を追加すればよいと理解します。すなわち

$FOAM_SRC/thermophysicalModels/solidSpecie/transport 以下に自製の物性のフォルダ・ソースコードを追加
$FOAM_SRC/thermophysicalModels/solidThermo/solidThermo/solidThermos.C を修正

同一環境をローカルに降ろして、ソースを修正できる環境を作る

元ソースコードをローカルにコピーします。一応、開発用フォルダ(dev-src)を作成。

$ mkdir dev-src
$ cd dev-src
$ cp -r $FOAM_SRC/thermophysicalModels/solidSpecie .
$ cp -r $FOAM_SRC/thermophysicalModels/solidThermo .

git管理を開始しておくと吉です。

solidSpecie/Make/filesの中を修正します。これを。

LIB = $(FOAM_LIBBIN)/libsolidSpecie

こう

LIB = $(FOAM_USER_LIBBIN)/libsolidSpecie

/solidSpecieのフォルダ内（Makeフォルダが配置されている階層）でコンパイルします。

$ wmake libso

念のため、コンパイルされているかを確認

$ ls $FOAM_USER_LIBBIN

続いて、solidThermo内。ローカルのsolidSpecieを参照できるようにします。これを

EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/solidSpecie/lnInclude \

こう

EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I../solidSpecie/lnInclude \

コンパイルには同じようにwmake libsoを使います。

プログラムの修正（準備）

$ cd solidSpecie/transport/
$ cp -r exponential/ usr

exponentialSolidTransportをもとに自製の物性プログラムを作成します。

$ sed -e "s/exponentialSolidTransport/usrSolidTransport/g" exponentialSolidTransport.C > usrSolidTransport.C
$ sed -e "s/exponentialSolidTransport/usrSolidTransport/g" exponentialSolidTransport.H > usrSolidTransport.H
$ sed -e "s/exponentialSolidTransport/usrSolidTransport/g" exponentialSolidTransportI.H > usrSolidTransportI.H

solidThermos.Cの内部を修正して、自製のusrSolidTransportを呼べるようにします。修正箇所は

#include "usrSolidTransport.H"

ヘッダファイルの追加

makeSolidThermo
(
    solidThermo,
    heSolidThermo,
    pureMixture,
    usrSolidTransport,
    sensibleEnthalpy,
    hPowerThermo,
    rhoConst,
    specie
);

usrSolidTransportを呼べるように。

修正したらコンパイルしなおします。（念のため）wcleanからのコンパイル。

$ wclean
$ wmake libso

プログラムの修正（中身）

usrSolidTransport.H

        static word typeName()
        {
            return "usr<" + Thermo::typeName() + '>';
        }

thermophysicalPropertiesを読み込むときに、"usr"のキーワードで読み込むことができるようになります。

コンパイルは、solidThermoディレクトリでwmake libsoたたけばできます（wcleanはいらない）。solidSpacie以下のディレクトリの変更も反映できます（なので、solidThermoディレクトリに居座ってコンパイルすればよい）

プログラムの修正（詳細）

変数a0を追加してみます。exponentialSolidTransport中のn0を編集する形で追加していきます。

usrSolidTransport.H, usrSolidTransportI.H, usrSolidTransport.Cの３ソースファイル中で、n0を探す⇒コピー⇒a0に編集を機械的に繰り返していくことで追加できます。

この段階で一度コンパイルしてみます。実行すると、変数a0が足りないと怒られるので、thermophysicalPropertiesを修正しておきます。

実行できることを確認した後は本格的に修正。熱物性計算の本体はusrSolidTransportI.H中の以下なので、この辺を修正していきます。

template<class Thermo>
inline Foam::scalar Foam::usrSolidTransport<Thermo>::kappa
(
    const scalar p, const scalar T
) const
{
    return (kappa0_*pow(T/Tref_, n0_));
}


template<class Thermo>
inline Foam::vector Foam::usrSolidTransport<Thermo>::Kappa
(
    const scalar p, const scalar T
) const
{
    const scalar kappa(kappa0_*pow(T/Tref_, n0_));
    return vector(kappa, kappa, kappa);
}

実行

OpenFOAMの実行プログラムは$FOAM_USER_LIBBIN中のライブラリを優先的に参照するようです。ですので、上記手順でソースコードを修正した場合、すでにコンパイル済みの実行ファイルでも修正部分が反映されます。

参考にしました

PENGUINITIS - 輸送特性の追加

PENGUINITIS - OpenFOAM による開発パターン

2021-03-31

pythonで関数のフィッティング

python

多項式で近似する場合

numpyの機能で実行できます。

import numpy as np
x            = np.arange(0, 1000, 1)
y_target = xxx

ndim = 3

# 
Coeffs=np.polyfit(x,y_target,3)

# 
plt.plot( x, np.poly1d(Coeffs)(x) )

一般の関数で近似する場合

scipyのoptimize, curve_fitを使います。

from scipy.optimize import curve_fit

def y_func(x,a,b):
    y = a*x+b
    return y

coefs, pcov = curve_fit(y_func, x, y_target)

plt.plot( x, y_func(x, coefs ) )

参考にしました

NumPy で回帰分析をする - Qiita

Numpy.polyfit を使ったカーブフィッティング - Qiita

2021-03-31

OpenFOAM【物性の温度依存性の設定１】

OpenFOAM FVM

chtMultiFegionFoam中で固体（solid）の物性の温度依存性を反映させるのが目的。熱伝導率（transport）、比熱（thermo）の温度依存性を考慮できます。thermophysicalPropertiesが設定ファイル。まずは、thermoTypeで関数の指定。

thermoType
{
    type            heSolidThermo;
    mixture         pureMixture;
    transport       polynomial;
    thermo          hPolynomial;
    equationOfState rhoConst;
    specie          specie;
    energy          sensibleEnthalpy;
}

（わざと）合わない設定をすると、以下出力されるので参考に設定します。

Valid solidThermo types :
6(heSolidThermo<multiComponentMixture<constIso<hConst<rhoConst<specie>>,sensibleEnthalpy>>> heSolidThermo<pureMixture<constAnIso<hConst<rhoConst<specie>>,sensibleEnthalpy>>> heSolidThermo<pureMixture<constIso<hConst<rhoConst<specie>>,sensibleEnthalpy>>> heSolidThermo<pureMixture<exponential<hPower<rhoConst<specie>>,sensibleEnthalpy>>> heSolidThermo<pureMixture<polynomial<hPolynomial<icoPolynomial<specie>>,sensibleEnthalpy>>> heSolidThermo<pureMixture<polynomial<hPolynomial<rhoConst<specie>>,sensibleEnthalpy>>>)
type           mixture                transport    thermo       equationOfState  specie  energy            

heSolidThermo  multiComponentMixture  constIso     hConst       rhoConst         specie  sensibleEnthalpy  
heSolidThermo  pureMixture            constAnIso   hConst       rhoConst         specie  sensibleEnthalpy  
heSolidThermo  pureMixture            constIso     hConst       rhoConst         specie  sensibleEnthalpy  
heSolidThermo  pureMixture            exponential  hPower       rhoConst         specie  sensibleEnthalpy  
heSolidThermo  pureMixture            polynomial   hPolynomial  icoPolynomial    specie  sensibleEnthalpy  
heSolidThermo  pureMixture            polynomial   hPolynomial  rhoConst         specie  sensibleEnthalpy

transportの設定

一定値（constIso）のほかに、指数形（exponential）、多項式（polynomial ）を設定できます。

    transport
    {   
        // exponentialの場合
        kappa0 1;
        n0 1.1;
        Tref 300;
        // polynomial の場合
        kappaCoeffs<8> (1e3 0 0 0 0 0 0 0)
    }

具体的な計算は、以下のソースコードを見ればわかります。

/src/thermophysicalModels/solidSpecie/transport/exponential/exponentialSolidTransportI.H
/src/thermophysicalModels/solidSpecie/transport/polynomial/polynomialSolidTransportI.H

thermoの設定

一定値（hConst）のほかに、指数形（hPower）、多項式（hPolynomial）を設定できます。

    thermodynamics
    {   
        Hf      0;  
        Sf      0;
        //  hPower の場合
        //  Cp(T) = c0*(T/Tref)^n0
        C0      1;
        Tref  300;
        n0    1.1;
        //  hPolynomial の場合
        CpCoeffs<8> (2000 0 0 0 0 0 0 0);
    }

2021-02-06

OpenFOAM【chtMultiRegionFoamで輻射：基礎編】

OpenFOAM FVM

chtMultiRegionFoamで輻射を入れる方法についてまとめ。

基本的な考え方

OpenFOAMで輻射を計算する場合、輻射に関する変数（G, IDefault, qr）が追加され、その変数を合わせて解くことになります。

輻射のモデルとしては、P1, viewFactor, FvDOMなどを選択可。viewFactorは形態係数を解析領域から計算させる必要があります（後述）。FvDOMは直接求解しているようで、計算コストを要します。

G: 入射の強度。前ステップの温度場から、MarshakRadiationの境界条件にもとづいて計算される。Marshak Radiationの境界条件の原著は下の論文（みたい）。 Phys. Rev. 71, 443 (1947) - Note on the Spherical Harmonic Method As Applied to the Milne Problem for a Sphere

IDefault：放射の強度。greyDiffusiveRadiationの境界条件から計算される。

qr：輻射による熱交換量。greyDiffusiveRadiationViewFactorの境界条件から計算される。viewFactorモデルの際に使用。

境界条件の設定

OpenFOAMの一般的な境界条件はこちら。

PENGUINITIS - 境界条件タイプリスト

一方で、constant中に配置するboundaryRadiationPropertiesでも輻射に関しての境界を設定できそうです。下のフォーラムによればどちらかで設定すればよさそう、とのことですが、tutorialのmultiRegionHeaterRadiationでは、どちらでも設定されているようです（例えばsystem/topAir/changeDictionaryDictとconstant/topAir/boundaryRadiationProperties）。これらの設定の反映については要確認だが、とりあえずどっちにも設定しとけばよいかも。

What is boundaryRadiationProperties -- CFD Online Discussion Forums

また、multiRegionHeaterRadiationの設定を確認すると、流体側のG,IDefault, qrには輻射に関する境界条件・界面の条件が付与されている一方で、固体側の境界条件には

    {   
        type            calculated;
        value           uniform 0;
    }

で付与されています。

温度Tの境界条件

compressible::turbulentTemperatureRadCoupledMixed を設定することで、輻射熱を温度の境界に反映させることができるようです。流体側と固体側とでqrNbr,qrをそれぞれ設定する必要があるようです。

流体側

        "topAir_to_.*"
        {
            type            compressible::turbulentTemperatureRadCoupledMixed;
            Tnbr            T;
            kappaMethod     fluidThermo;
            qrNbr           none;
            qr              qr;
            value           uniform 300;
        }

固体側

        heater_to_topAir
        {
            type            compressible::turbulentTemperatureRadCoupledMixed;            Tnbr            T;
            kappaMethod     solidThermo;
            qrNbr           qr; 
            qr              none;
            value           uniform 300;
        }

PENGUINITIS - 流体・固体熱連成解析

輻射の境界条件の設定は、turbulentTemperatureRadCoupledMixedのようですが、もとは以下フォルダに格納です。

src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/turbulentTemperatureRadCoupledMixed

Understanding temperature coupling BCs -- CFD Online Discussion Forums

viewFactorモデルの使用について

viewFactorを計算するには、preprocessで形態係数に関する設定が必要です。なので、constant中にviewFactorsDictを配置して、preprocessが必要です。multiRegionHeaterRadiationを参考にすると、以下のコマンドが必要。

# Agglomerate patch faces
for region in bottomAir topAir
do
    runApplication -s $region \
        faceAgglomerate -region $region -dict constant/viewFactorsDict
done

# Generate view factors
for region in bottomAir topAir
do
    runApplication -s $region \
        viewFactorsGen -region $region
done

あわせて、境界条件の設定も必要なようです。

boundary
{
    ".*"
    {
        inGroups        2(wall viewFactorWall);
    }
}

【参考】multiRegion＋輻射に関するチュートリアル

各モデルと初期条件で設定されている変数のまとめ。参照しているのはOpenFOAM-v1912です。

/heatTransfer/buoyantSimpleFoam/hotRadiationRoom

SIMPLE+P1
G T U alphat epsilon k nut p p_rgh

/heatTransfer/buoyantSimpleFoam/hotRadiationRoomFvDOM

SIMPLE+fvDOM
G IDefault T U alphat epsilon k nut p p_rgh

/heatTransfer/chtMultiRegionFoam/multiRegionHeater

multiregion+PIMPLE
T U epsilon k p p_rgh

/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeaterRadiation

multiregion+SIMPLE+viewFactor
G IDefault T U epsilon k p p_rgh qr

/heatTransfer/buoyantPimpleFoam/thermocoupleTestCase

PIMPLE+fvDOM;
IDefault T U p p_rgh

参考にしました。

OpenFOAM全般

PENGUINITIS - 計算の確認と計算結果の処理

http://www.opencae.or.jp/wp-content/uploads/2015/06/C-_%E7%AC%AC2%E5%9B%9E.pdf

http://eddy.pu-toyama.ac.jp/%E3%82%AA%E3%83%BC%E3%83%97%E3%83%B3CAE%E5%8B%89%E5%BC%B7%E4%BC%9A-%E5%AF%8C%E5%B1%B1/?action=cabinet_action_main_download&block_id=99&room_id=1&cabinet_id=1&file_id=40&upload_id=121

伝熱解析全般

http://opencae.gifu-nct.ac.jp/pukiwiki/index.php?plugin=attach&pcmd=open&file=buoyantBoussinesqSimleFoam1.pdf&refer=%C2%E8%A3%B1%A3%B8%B2%F3%CA%D9%B6%AF%B2%F1%A1%A7H241103

OpenFOAM 標準ユーティリティー一覧 - XSim

chtMultiRegionFoam情報

PENGUINITIS - 流体・固体熱連成解析

chtMultiRegionSimpleFoam：流体-固体熱連成解析（定常） - Onion Salad - OpenFOAMチュートリアルドキュメント作成プロジェクト

chtMultiRegionSimpleFoam:熱交換器とファン - Onion Salad - OpenFOAMチュートリアルドキュメント作成プロジェクト

OpenFOAMのマルチリージョン（標準ソルバー編その1 ケースディレクトリ） - Qiita

OpenFOAM マルチリージョンの取り扱い（境界条件と並列処理） - Qiita

OpenFOAM-v1712のCHT部分の変更点（sub-looping） - Qiita

OpenFOAMでの輻射解析

https://scielo.conicyt.cl/pdf/ingeniare/v26n4/0718-3305-ingeniare-26-04-00546.pdf

輻射を考慮した直方体領域内の自然対流（fvDOM モデル） - XSim

http://www.tfd.chalmers.se/~hani/kurser/OS_CFD_2009/AlexeyVdovin/Radiation_in_OpenFoam_final.pdf

OpenFOAM: API Guide: Radiation models

How to deal with Conduction, Convection and radiation in OpenFOAM | CFD WITH A MISSION

https://www.foamacademy.com/wp-content/uploads/2018/03/thermo_training_handout_public.pdf

Difference between IDefault and G (radiation) -- CFD Online Discussion Forums

chtMultiRegionFoam: crash with fvDOM and unstructured mesh -- CFD Online Discussion Forums

PENGUINITIS - 非圧縮性流体ソルバーにおける温度計算

2021-01-26

OpenFOAM【pythonでblockMesh】

OpenFOAM python

penguinitis.g1.xrea.com

春日氏作成のpythonプログラムblockMesh.pyを使ってみました。ので、簡単なメモ（OpenFOAMも久しぶりなので）。

$ wget http://penguinitis.g1.xrea.com/study/OpenFOAM/blockMeshPython/blockMeshPython-v1912-20200612.tar.gz 
$ tar -zxvf blockMeshPython-v1912-20200612.tar.gz

でダウンロード＋解凍。チュートリアルプログラムが二つほど入っています。createMesh.pyとcreateMesh2.py。とりあえず、createMesh.pyを動かしてみます。

ディレクトリとして、/constant、/systemを用意します。

$ python createMesh.py

をたたくことで、blockMeshDictを作成します。作成場所はcreateMesh.pyの中に記載。/system中にcontrolDictを配置します。内容は適当でもよい。

$ blockMesh

をたたくと/constantフォルダの中に、polyMeshディレクトリ・meshファイルが作成されています。

$ paraFoam

でメッシュを確認できます。

viewFactorモデルの考え方

viewFactorの実際の導入についての方針

viewFactorを追加するときに設定する部分

デバッグの覚書（トラブルシューティング）

輻射を受ける対象物（固体）の温度が上がり続ける（下がり続ける）

温度が振動して、最後は異常終了で計算停止

参考にしました

OpenFOAM全般

chtMultiRegionFoamに関して

OpenFOAMの輻射モデルについて

フォーラム情報

参考

前準備（方針の確認）

同一環境をローカルに降ろして、ソースを修正できる環境を作る

プログラムの修正（準備）

プログラムの修正（中身）

プログラムの修正（詳細）

実行

参考にしました

多項式で近似する場合

一般の関数で近似する場合

参考にしました

transportの設定

thermoの設定

基本的な考え方

境界条件の設定

温度Tの境界条件

viewFactorモデルの使用について

【参考】multiRegion＋輻射に関するチュートリアル

参考にしました。

OpenFOAM全般

伝熱解析全般

chtMultiRegionFoam情報

OpenFOAMでの輻射解析