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

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

前準備(方針の確認)

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

$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の3ソースファイル中で、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 による開発パターン

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

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

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);
    }  

OpenFOAM【chtMultiRegionFoamで輻射:基礎編】

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 - 非圧縮性流体ソルバーにおける温度計算

OpenFOAM【pythonでblockMesh】

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

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

WSL2のインストール

今更ですが。基本的には下サイト通りにやれば、エラーも出ずにうまくいくはずです(私はまずはちゃんとサイトの通りやらなかったので、以下でエラーが出ています。)

docs.microsoft.com

管理者権限でpowershellを開いて実行します。

PS C:\WINDOWS\system32> wsl --set-default-version 2

以下エラーが出ました。

Windows の仮想マシン プラットフォーム機能を有効にして、BIOS で仮想化が有効になっていることを確認してください。
詳細については、https://aka.ms/wsl2-install を参照してください

仮想マシン関係の機能を有効化します。コマンド打ったあとは、一旦再起動。

PS C:\WINDOWS\system32> dism.exe /online /enable-feature /featurename:VirtualMachinePlatform /all /norestart

もう一度バージョンを上げるコマンドを入れると、今度は以下のエラーが出ます。

WSL 2 を実行するには、カーネル コンポーネントの更新が必要です。

素直に、x64 マシン用 WSL2 Linux カーネル更新プログラム パッケージをインストールします。今度こそ。

PS C:\WINDOWS\system32> wsl --set-default-version 2

以下、コマンドが出ます。

WSL 2 との主な違いについては、https://aka.ms/wsl2 を参照してください

Powershellでバージョンを確認します。

PS C:\Users\> wsl -l -v
  NAME      STATE           VERSION
* Ubuntu    Stopped         1

以下コマンドでバージョン変更のようです。

wsl --set-version <NAME> <VERSION>

自分はubuntuをv2に変更

PS C:\Users\> wsl --set-version Ubuntu 2
PS C:\Users\> wsl -l -v

バージョンアップ中(変換)には時間がかかります。寝て待ったらインストールされてた。

参考にしました。

WSL1からWSL2への移行 - Qiita

WSL2の環境構築手順 - Qiita

WSL バージョンを変更する方法 | SEECK.JP サポート

WSL2でXサーバ

前回。

matsuo-san.hatenablog.com

WSL2になったので、ソフトもxming→VcXsrvに変更しました。

kb.seeck.jp

VcXsrvのインストール

以下からソフトウェアをダウンロードします。

sourceforge.net

ダウンロードしたら、ポチポチと設定してインストール(あんまり気にするところはない。)

サーバ接続の設定

ここを参考に設定していきます。

odaryo.hatenablog.com

まずは、x11関係をインストール

$ sudo apt install x11-apps

ディスプレイ関係を、.bashrcに追加しておきます。

export DISPLAY=$(cat /etc/resolv.conf | grep nameserver | awk '{print $2}'):0.0

早速チェックしますが、エラーが出ます。

$ xeyes
Authorization required, but no authorization protocol specified
Error: Can't open display: xxx.xx..0.1:0.0

起動時に「Disable access control」にチェックをいれればokでした。

f:id:matsuo_san:20210124113312p:plain

こちら参照にしました。

rin-ka.net