【（２）SALOMEでヒートシンクの熱流体固体連成メッシュ作成】chtMultiRegionSimpleFoamで計算実行

こんにちは（@t_kun_kamakiri）

OpenFOAMの熱流体固体連成によるヒートシンクの熱流体解析の設定方法を解説します。

前回メッシュ作成まで行いましたので、メッシュ作成編を見たい方はこちらをどうぞ。

【（１）SALOMEでヒートシンクの熱流体固体連成メッシュ作成】chtMultiRegionSimpleFoam用マルチリージョンモデル構築手順 こんにちは（@t_kun_kamakiri） OpenFOAMの熱流体固体連成によるヒートシンクの熱流体解析の設定方法を解説しま...

使用するソルバはchtMultiRegionSimpleFoamで熱流体と固体を連成させた定常解析用のソルバになります。

OpenFOAMv2412のバージョンにおいて現在は、浮力を考慮した熱流体解析ソルバは以下のものがあります。

ソルバ名	定常 / 非定常	主な特徴	主な用途・備考
chtMultiRegionSimpleFoam	定常	多領域の熱流体連成（非圧縮性）	構造物と流体の熱連成、定常問題向け
chtMultiRegionFoam	非定常	多領域・非圧縮性の時間依存連成解析	熱伝導と流体連成の過渡解析
chtMultiRegionTwoPhaseEulerFoam	非定常	多領域 + 二相流（Euler-Euler） + 熱連成	液体/気体の二相流＋固体との熱連成に対応
buoyantSimpleFoam	定常	単領域の自然対流・熱伝導	密度差による自然対流、定常解析向け
buoyantPimpleFoam	非定常	上記の非定常版	過渡的な自然対流解析
buoyantBoussinesqSimpleFoam	定常	Boussinesq近似（小温度差）	温度変化が小さい自然対流
buoyantBoussinesqPimpleFoam	非定常	上記の時間依存版	過渡的自然対流、簡略化モデル

この中のchtMultiRegion系は流体（液体、気体）・固体の複数領域の熱流体解析を行うために設計されています。
ただ、テキストベースで解析設定しようとすると、めちゃくちゃ大変です。
以下にOpenFOAMが用意しているチュートリアルを解説した記事がありますので、一度読んでみると良いでしょう。

以前にsappyHexMeshでヒートシンクの固体・流体連成の解析を紹介しました。

【ヒートシンクの熱流体(1)】OpenFOAMの熱流体固体連成のメッシュ作成こんにちは（@t_kun_kamakiri） OpenFOAMの熱流体固体連成によるヒートシンクの熱流体解析の設定方法を解説します。 ...

【ヒートシンクの熱流体(2)】OpenFOAMの熱流体固体連成の計算実行 こんにちは（@t_kun_kamakiri） OpenFOAMの熱流体固体連成によるヒートシンクの熱流体解析の設定方法を解...

今回は違った方法でモデル作成を行いヒートシンクの解析まで行います。
領域分割やメッシュ作成はSALOMEを使用します。

OpenFOAM v2412(WSL Ubuntu 22.04)
SALOME 9.15.0

OpenFOAM用の設定ファイルの準備

SALOMEで問題なくメッシュ作成と境界の設定が終わったらOpenFOAMの形式に変換します。

以下のようにunvファイルとして出力します。

出力したunvファイルをOpenFOAM形式に変換するにはダミーのケースファイルが必要です。適当なチュートリアルからコピーしてきます。

今回は固体・熱流体の定常解析で計算するのでchtMultiRegionSimpleFoamからコピーしてきます。

cp -r /home/kamakiri/repo/openfoam/tutorials/heatTransfer/chtMultiRegionSimpleFoam/cpuCabinet/

$FOAM_TUTORIALS = /home/kamakiri/repo/openfoam/tutorials

salomeで作成したunvファイルをケースファイルに保存します。

cp -r model/Mesh_1.unv .

今フォルダ構成を確認してくと以下のようになっています。

.
├── 0.orig
│   ├── domain0
│   ├── v_CPU
│   └── v_fins
├── Allclean
├── Allrun
├── Allrun.pre
├── Mesh_1.unv
├── constant
│   ├── domain0
│   ├── g
│   ├── regionProperties
│   ├── v_CPU
│   └── v_fins
├── model
│   ├── Mesh_1.unv
│   ├── Study1.hdf
│   ├── heatsinkl.step
│   └── model.FCStd
└── system
    ├── blockMeshDict
    ├── controlDict
    ├── createBafflesDict
    ├── decomposeParDict
    ├── domain0
    ├── fvSchemes
    ├── fvSolution
    ├── meshQualityDict
    ├── probes
    ├── snappyHexMeshDict
    ├── surfaceFeatureExtractDict
    ├── topoSetDict.f1
    ├── v_CPU
    └── v_fins

0.origをコピーします。

cp -r 0.orig 0

すでに0フォルダにあるフォルダをひとまずoldフォルダに避難させます。

mkdir -p 0/old && mv 0/{domain0,v_CPU,v_fins} 0/old
mkdir -p constant/old && mv constant/{domain0,v_CPU,v_fins} constant/old 
mkdir -p system/old && mv system/{domain0,v_CPU,v_fins} system/old 

ここからがマルチリージョンソルバの手順のコツですが、以下のように0フォルダに計算する物理量のファイルを用意します。

今回はk-ωSSTの乱流モデルで計算したいので必要なのは、U,T,p_rgh,p,k,omega,nut,alphatになります。

0
├── T
├── U
├── alphat
├── cellToRegion
├── epsilon
├── k
├── nut
├── old
├── omega
├── p
├── p_rgh
└── rho

中身はすべて以下としておきます。
例としてUの中身を示しています。物理量によってdimensionsが異なるので注意してください。適当なチュートリアルから各物理量をコピーしてきて、以下のように編集しておくとよいでしょう。

FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0.01 0 0);

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    ".*"
    {
        type    calculated;
        value   $internalField;
    }
}

こうしておくと非常に便利です。

このあと0, constant, systemフォルダに各領域のフォルダを作っていく必要があるのですが、これを手動で行うとめちゃくちゃ面倒になります。

ぜひ、以下の方法を覚えておいてください。

OpenFOAMの形式に変換

ここからSALOMEで作成したメッシュ情報をOpenFOAMに変換します。

以下のコマンドでOpenFOAMの形式に変換します。

ideasUnvToFoam Mesh_1.unv

エラーがなければconstant/polyMeshフォルダが作成されOpenFOAMのメッシュ情報が作られています。

※今回は必要ありませんが、salomeではmm単位のつもりで作成しましたが、実は数値しか意味がなく、OpenFOAMでの計算がm単位なので、必要に応じて以下のコマンドでスケールを1/1000倍にします。

transformPoints -scale '(0.001 0.001 0.001)'

再度、モデルの大きさを確認します。

まだメッシュを変換しただけで領域ごとにメッシュが分かれているわけではありません。以下のコマンドで領域を分割します。

splitMeshRegions -cellZones -overwrite

領域分割によりフォルダ構成は以下となっています。
splitMeshRegionsコマンドにより、0, constant. system内に領域ごとの設定ファイルが作成されますconstant/air(heatsink）, system/air(heatsink）のフォルダは元々作成していましたが、0/air(heatsink)はフォルダが生成されたと思います。

気体領域と固体領域で共通している面は、heatsink_to_airやair_to_heatsinkの名前が付けられます。

境界条件の設定

splitMeshRegionsコマンドにより境界条件の設定ファイルの中身も変わっています。
例として流速の0フォルダを示します。

0/air/U

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0.01 0 0);

boundaryField
{
    YMin
    {
        type            calculated;
        value           uniform (0.01 0 0);
    }
    ZMax
    {
        type            calculated;
        value           uniform (0.01 0 0);
    }
    YMax
    {
        type            calculated;
        value           uniform (0.01 0 0);
    }
    XMin
    {
        type            calculated;
        value           uniform (0.01 0 0);
    }
    ZMin
    {
        type            calculated;
        value           uniform (0.01 0 0);
    }
    XMax
    {
        type            calculated;
        value           uniform (0.01 0 0);
    }
    air_to_heatsink
    {
        type            calculated;
        value           uniform (0 0 0);
    }
}

このように自動的にair領域の境界条件が作成されています。

ここで注意として、流体領域と固体領域で必要になるファイルが異なります。

• 流体領域：0, T, p_rgh, p
• 固体領域：T, p

k-ωSST乱流モデルを使用するのであれば、流体領域にさらにk, omega, nut, alphatが必要になります。

さてここから境界条件を手動で競ってしていきたいところですが、まだ手動で行うには少々面倒なのでchangeDictionaryDictを使用します。

以下では境界条件の設定と同時にconstant/<region>/polyMesh/boundaryの境界のタイプも変更しています。

system/air/changeDictionaryDict

boundary
{
    ".*"
    {
        type    wall;
    }

    ".*_to_.*"
    {
        type    mappedWall;;
    }
    YMin
    {
        type    patch;
    }

    YMax
    {
        type    patch;
    }
}

k
{
    internalField   uniform 0.00375;

    boundaryField
    {
        ".*"
        {
            type    kqRWallFunction;
            value   $internalField;
        }

        YMin
        {
            type            inletOutlet;
            inletValue      uniform 0.00375;
            value           uniform 0.00375;
        }

        YMax
        {
            type            turbulentIntensityKineticEnergyInlet;
            value           $internalField;
            intensity       0.03;
        }
    }
}

omega
{
    internalField   uniform 5.6;////I=0.05,k=3/2(IU)^2, omega=k^0.5/(Cnu^0.25*L),L=20

    boundaryField
    {
        ".*"
        {
            type    omegaWallFunction;
            value   $internalField;
        }

        YMin
        {
            type    inletOutlet;
            inletValue    uniform 0.67082;
            value    uniform 6.7082e-05;
        }

        YMax
        {
            type    turbulentMixingLengthFrequencyInlet;
            mixingLength    0.01;
            value    uniform 0.67082;
        }
    }
}

nut
{
    internalField   uniform 6000;

    boundaryField
    {
        ".*"
        {
            type    nutkWallFunction;
            value   uniform 0;
        }

        YMin
        {
            type    calculated;
            value   uniform 0;
        }

        YMax
        {
            type    calculated;
            value   uniform 0;
        }
    }
}

U
{
    internalField   uniform ( 0 -0.5 0 );

    boundaryField
    {
        ".*"
        {
            type    fixedValue;
            value   uniform (0 0 0);
        }

        YMin
        {
            type            pressureInletOutletVelocity;
            value           $internalField;
        }

        YMax
        {
            type            fixedValue;
            value           $internalField;
        }
    }
}

T
{
    internalField   uniform 293.15;

    boundaryField
    {
        ".*"
        {
            type    fixedValue;
            value   uniform 293.15;
        }

        YMin
        {
            type    fixedValue;
            value   uniform 293.15;
        }

        YMax
        {
            type    zeroGradient;
        }

        ".*_to_.*"
        {
            type            compressible::turbulentTemperatureRadCoupledMixed;
            value           uniform 293.15;
            Tnbr            T;
            kappaMethod     fluidThermo;
            qrNbr           none;
            qr              none;//qr;
            thermalInertia  true;
            kappa           none;
        }
    }
}

p_rgh
{
    internalField   uniform 100000;

    boundaryField
    {
        ".*"
        {
            type        fixedFluxPressure;
            value       uniform 100000;
        }

        YMin
        {
            type            totalPressure;
            value           uniform 100000;
            p0              uniform 100000;
        }

        YMax
        {
            type            fixedFluxPressure;
            value           uniform 100000;
        }
    }
}

p
{
    internalField   uniform 100000;

    boundaryField
    {
        ".*"
        {
            type            calculated;
            value           uniform 100000;
        }
    }
}

alphat
{
    internalField   uniform 0;

    boundaryField
    {
        ".*"
        {
            type    compressible::alphatWallFunction;
            Prt     0.85;
            value   uniform 0;
        }

        YMin
        {
            type    calculated;
            value   uniform 0;
        }

        YMax
        {
            type    calculated;
            value   uniform 0;
        }
    }
}

G
{
    internalField   uniform 0;

    boundaryField
    {
        ".*"
        {
            type    MarshakRadiation;
            value   uniform 0;
        }
    }
}

qr
{
    internalField   uniform 0;

    boundaryField
    {
        ".*"
        {
            type    calculated;
            value   uniform 0;
        }
    }
}

IDefault
{
    internalField   uniform 0;

    boundaryField
    {
        ".*"
        {
            type    greyDiffusiveRadiation;
            value   uniform 0;
        }
    }
}

G, qr, IDefaultなどはふく射の設定時に使うものなのですが、一応書いておきます（ふく射計算する際には使用する場合もあります）

system/heatsink/changeDictionaryDict

boundary
{
    ".*"
    {
        type    wall;
    }

    ".*_to_.*"
    {
        type    mappedWall;;
    }
    YMin
    {
        type    patch;
    }

    YMax
    {
        type    patch;
    }
}

T
{
    internalField   uniform 293.15;  // 初期固体温度

    boundaryField
    {
        // すべての固体外表面は固定温度（必要に応じて zeroGradient に変更可）
        ".*"
        {
            type    fixedValue;
            value   uniform 293.15;
        }
        heatSource
        {
            type            externalWallHeatFluxTemperature;
            value           $internalField;
            kappaMethod     fluidThermo;
            mode            power;
            Ta              $internalField;
            Q               uniform 2000;
            kappaName       none;
        }
        // 流体-固体連成面
        ".*_to_.*"
        {
            type            compressible::turbulentTemperatureRadCoupledMixed;
            Tnbr            T;
            kappaMethod     solidThermo;
            qrNbr           none;
            qr              none;
            kappa           none;
            thermalInertia  true;
            value           uniform 293.15;
        }
    }
}

設定を終わった後に領域を指定する必要があるので、以下の設定ファイルを変更します。

constant/regionProperties

regions
(
    solid ( heatsink  )
    fluid ( air )
);

ここから以下のコマンドで領域ごとに境界条件を指定します。

for region in $(foamListRegions); do     changeDictionary -region "$region"; done

このコマンドは;で区切ることで次のコマンドを打っていることになるので、スクリプトとして書くならばインデントをきれいにして以下のように書きます。

for region in $(foamListRegions)
do
    changeDictionary -region "$region"
done

これにより協会条件が各領域の各面で設定されたのではないかと思います。

物性値の設定

ここから各領域の物性値を設定します。
物性値はconstantフォルダにあります。

cp -r constant/old/domain0/thermophysicalProperties constant/air/
cp -r constant/old/v_CPU/thermophysicalProperties constant/heatsink/

まずは流体領域の空気の物性値を指定します。

constant/air/thermophysicalProperties

thermoType
{
    type            heRhoThermo;
    mixture         pureMixture;
    transport       const;
    thermo          hConst;
    equationOfState perfectGas;
    specie          specie;
    energy          sensibleEnthalpy;
}

mixture
{
    specie
    {
        nMoles          1;
        molWeight       28.966;
    }

    thermodynamics
    {
        Cp              1006.43;
        Hf              0;
    }

    transport
    {
        mu              1.846e-05;
        Pr              0.706414;
    }
}

固体領域はアルミニウム（Al）の固体物性値とします。

constant/heatsink/thermophysicalProperties

thermoType
{
    type            heSolidThermo;
    mixture         pureMixture;
    transport       constIso;
    thermo          hConst;
    equationOfState rhoConst;
    specie          specie;
    energy          sensibleEnthalpy;
}

mixture
{
    specie
    {
        nMoles          1;
        molWeight       26.98;
    }

    thermodynamics
    {
        Hf              0;
        Sf              0;
        Cp              871;
    }

    transport
    {
        kappa           202.4;
    }

    equationOfState
    {
        rho             2719;
    }
}

また流体と固体では状態方程式を扱うモデルが違うことに注意してください。

• 流体領域：equationOfState perfectGas;理想気体の状態方程式
• 固体領域：equationOfState rhoConst;密度一定

乱流モデルの設定

乱流モデルもチュートリアルの設定からコピーしてきます。

cp -r constant/old/domain0/turbulenceProperties constant/air/

乱流モデルの設定はk-ωSSTにしたいので以下のように変更します。

simulationType  RAS;

RAS
{
    RASModel        kOmegaSST;//realizableKE;
    turbulence      on;
    printCoeffs     on;
}

重力の設定

重力はz方向の負の方向に指定します。
これの向きを間違えると浮力流れの向きがおかしくなります。

constant/g

dimensions      [0 1 -2 0 0 0 0];

value           (0 0 -9.81);

離散化・代数ソルバ

こちらもチュートリアルからコピーしてきます。

cp -r system/old/domain0/{fvSchemes,fvSolution} system/air/
cp -r system/old/v_CPU/{fvSchemes,fvSolution} system/heatsink/

離散化スキームや代数ソルバにomegaを解く設定がないため追加します。

system/air/fvSchemes

ddtSchemes
{
    default         steadyState;
}

divSchemes
{
    div(phi,U)      bounded Gauss upwind;
    div(phi,k)      bounded Gauss linearUpwind grad(k);
    div(phi,epsilon) bounded Gauss linearUpwind grad(epsilon);
    div(phi,omega) bounded Gauss linearUpwind grad(omega); //追加
    div((muEff*dev2(T(grad(U))))) Gauss linear;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
    div(phi,K)      bounded Gauss linearUpwind default;
    div(phi,h)      bounded Gauss upwind;
    div(phid,p)     bounded Gauss upwind;
}

laplacianSchemes
{
    default         Gauss linear limited 0.333;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         limited 0.333;
}

gradSchemes
{
    default         Gauss linear;
}

fluxRequired
{
    default         no;
    pCorr           ;
    p_rgh           ;
}

wallDist
{
    method          meshWave;
    nRequired       false;
}

system/air/fvSolution

solvers
{
    rho
    {
        solver          PCG
        preconditioner  DIC;
        tolerance       1e-7;
        relTol          0;
    }

    p_rgh
    {
        // solver           PBiCGStab;
        // preconditioner   FDIC;
        // tolerance        1e-7;
        // relTol           0.01;
        // smoother         GaussSeidel;
        // cacheAgglomeration true;
        solver           GAMG;
        tolerance        1e-7;
        relTol           0.01;
        smoother         GaussSeidel;
    }


   "(U|k|h|epsilon|omega)"
    {
        solver           PBiCGStab;
        preconditioner   DILU;
        tolerance        1e-7;
        relTol           0.05;
    }
}

SIMPLE
{
    momentumPredictor true;
    nNonOrthogonalCorrectors 0;
    frozenFlow      false;

    residualControl
    {
        default         1e-7;
    }
}

relaxationFactors
{
    fields
    {
        p_rgh           0.7;
        rho             1;
    }

    equations
    {
        U               0.3;
        h               0.9;
        k               0.7;
        epsilon         0.7;
    }
}

発散する場合は緩和係数の数値（relaxationFactors）を少し小さめにするなど調整が必要です。

計算制御の設定

今回は定常解析で2000ステップ（endTimeの値）まで、結果の出力は100ステップごと（writeIntervalの値）にします。

application     chtMultiRegionSimpleFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         2000;

deltaT          1;

writeControl    timeStep;

writeInterval   100;

purgeWrite      0;

writeFormat     binary;

writePrecision  8;

writeCompression off;

timeFormat      general;

timePrecision   8;

runTimeModifiable true;

functions
{
    #include    "probes"
}

probeの設定がされていますが、こちらは今回の仕様に変更する必要があるのですが、何も変えずに行きます（つまり効かない）

計算実行

モデルはそこまで大きくないので並列させずに計算します。

chtMultiRegionSimpleFoam

エラーがなく順調に進めばOKです。

一発で設定が上手くいくことがなくエラーで苦戦しましたが、エラーをなくすために以下を注意して修正しました。

• omegaのinternalFieldが大きい
• 流体・固体境界の温度（uniform 293.15;）が周りと違う
• 行列ソルバのp_rghでPBiCGStabだと発散しGAMG だと安定

結果の確認