流体・固体熱連成解析

2016年12月3日

はじめに

流体・固体熱連成解析の手順について。

使用バージョン

OpenFOAM 4.x (blueCFD)

ファイル

モデル

次のようなモデルを想定する。

2 次元矩形領域を考え、上を固体、下を流体とする。流体は左から 500 K で流入し、右から出て行く。固体の上端は 300 K 固定とする。流体・固体間の熱伝達を考慮する。

ソルバーには chtMultiRegionSimpleFoam を用いる。

設定手順

  1. $FOAM_TUTORIALS/heatTransfer/chtMultiRegionSimpleFoam からケースをコピーする。
  2. 通常通りメッシュを作成し、気体・固体各領域にゾーンを設定する。
    ここでは air と solid とする。
  3. splitMeshRegions により領域を分ける。
    $ splitMeshRegions -cellZones -overwrite
    
    system、constant、0 の中に領域ごとのディレクトリができる。
  4. constant/regionProperties を書き換える。
    FoamFile
    {
        version     2.0;
        format      ascii;
        class       dictionary;
        location    "constant";
        object      regionProperties;
    }
    
    regions
    (
        fluid       (air)
        solid       (solid)
    );
    
  5. 領域ごとに境界条件を設定する。ここでは直接設定するのではなく、changeDictionary を用いる。system の各領域のディレクトリに changeDictionaryDict を用意する。
    air/changeDictionaryDict
    FoamFile
    {
        version     2.0;
        format      ascii;
        class       dictionary;
        object      changeDictionaryDict;
    }
    
    ...
    
    T
    {
        internalField   uniform 300;
        boundaryField
        {
            ".*"
            {
                type            zeroGradient;
            }
            empty
            {
                type            empty;
                ~value;
            }
            
            ...
    
            "air_to_.*"
            {
                type            compressible::turbulentTemperatureCoupledBaffleMixed;
                Tnbr            T;
                kappaMethod     fluidThermo;
                kappa           none;
                value           uniform 300;
            }
        }
    }
    
    ...
    

    solid/changeDictionaryDict
    FoamFile
    {
        version     2.0;
        format      ascii;
        class       dictionary;
        object      changeDictionaryDict;
    }
    
    boundary
    {
        top
        {
            type            patch;
        }
    }
    
    T
    {
        internalField   uniform 300;
        boundaryField
        {
            ".*"
            {
                type            zeroGradient;
            }
            empty
            {
                type            empty;
                ~value;
            }
            
            ...
    
            "solid_to_.*"
            {
                type            compressible::turbulentTemperatureCoupledBaffleMixed;
                Tnbr            T;
                kappaMethod     solidThermo;
                kappa           none;
                value           uniform 300;
            }
        }
    }
    
    領域間の界面では "air_to_solid" や "solid_to_air" といった名前の境界ができる。界面では連成用の境界条件を指定する。
    changeDictionary は領域ごとに次のように実行する。
    $ changeDictionary -region air
    $ changeDictionary -region solid
    
  6. constant の中の設定を領域ごとに行う。
  7. system の中の設定を領域ごとに行う。fvSchemes、fvSolution、decomposeParDict は領域ごとに必要であるが、system 直下にも一応用意しておく必要がある。
  8. chtMultiRegionSimpleFoam を実行。
    $ chtMultiRegionSimpleFoam
    
  9. 計算できたら、ポスト処理用に領域ごとの .OpenFOAM を作成する。
    $ paraFoam -region air -touch
    $ paraFoam -region solid -touch
    
    ここでは "simpleCHT{air}.OpenFOAM"、"simpleCHT{solid}.OpenFOAM" といったファイルができる。
  10. ParaView でポスト処理を行う。paraFoam ではなく paraview を実行し、領域ごとの .OpenFOAM ファイルを読み込む (複数選択で一度に読み込める)。それぞれの領域に対してフィルタを設定するか、全体に設定したければ Group Datasets フィルタで一つにまとめる。

計算結果

熱伝達を評価してみる。固体の温度は、計算ログから 300.764 K とする。

Min/max T:min(T) [0 0 0 1 0 0 0] 300 max(T) [0 0 0 1 0 0 0] 300.764

壁面熱流束は、wallHeatFlux により以下の通りである。

$ wallHeatFlux -region air -latestTime

Wall heat fluxes [W]
bottom
    convective: 0
    radiative:  -0
    total:      0
air_to_solid
    convective: -20.1839
    radiative:  -0
    total:      -20.1839



$ wallHeatFlux -region solid -latestTime

Wall heat fluxes [W]
solid_to_air
    convective: 20.176
    radiative:  -0
    total:      20.176

壁面は幅 1 m だが、奥行を 0.1 m としているので、熱流束は 201.8 W/m2 である。これより、流体から固体への熱伝達率は、流体の温度を 500 K として 201.8/(500 - 300.764) = 1.013 W/m2-K である。

function object で熱拡散率 thermo:alpha (熱伝導率/比熱) を出力している。流体の thermo:alpha の値が 2.57143e-5 で、比熱は 1000 であるので、熱伝導率 kf は 2.57e-2 である。固体の熱伝導率 ks は 80 であるから、壁面熱伝導率 kw は kf と ks の調和平均として以下のようになる。

kw = 2*kf*ks/(kf + ks) = 5.14e-2

界面を挟んだセル中心間距離は 1/20 なので、熱伝達率は 5.14e-2*20 = 1.028 W/m2-K で、上の計算とだいたい合う。

熱伝達境界条件

領域界面 (上の例では air_to_solid や solit_to_air) では熱伝達境界条件を設定できる。熱伝達境界条件には以下のものがある。

  • compressible::turbulentTemperatureCoupledBaffleMixed
  • compressible::turbulentTemperatureRadCoupledMixed

設定は以下のように行う。

            "air_to_.*"
            {
                type            compressible::turbulentTemperatureCoupledBaffleMixed;
                Tnbr            T;
                kappaMethod     fluidThermo;
                kappa           none;
                value           uniform 300;
            }

ここで kappa は熱伝導率を与える方法を指定する。ふつう、流体領域なら "fluidThermo"、固体領域なら "solidThermo" を指定する。任意のフィールドを熱伝導率として指定したければ "lookup" を指定し、kappaName でフィールド名を指定する。

デフォルトでは界面での熱流束は単純に領域の熱伝導率で決まるが、界面に任意の熱伝導率の層を挟むこともできる。

        thicknessLayers (0.1);
        kappaLayers     (0.1);

thicknessLayers で各層の厚みを指定し、kappaLayers で各層の熱伝導率を指定する。それぞれの熱伝導率を層の厚みで割ったものが熱伝達率になる。界面には 2 つの境界があるが、両方に設定すると設定が 2 倍になるので、どちらか一方に設定すればよい。

並列計算

並列計算の場合は、以下のようにする。

$ decomposePar -allRegions
$ mpirun -np 4 chtMultiRegionSimpleFoam -parallel
$ reconstructPar -allRegions

decomposePar、reconstructPar にそれぞれ "-allRegions" オプションをつける。

3次元非定常計算の例

3次元非定常計算の例を下図に示す。非定常計算の場合は chtMultiRegionFoam を用いる。

任意の熱伝達率を与える方法

領域界面で任意の熱伝達を与えるには、境界条件 compressible::turbulentTemperatureCoupledBaffleMixed の kappa を "lookup" を指定し、kappaName に任意のフィールドを指定する。

まず、任意の熱伝導率を指定するためのフィールドを扱えるように、chtMultiRegionSimpleFoam を改造する。

chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H

    PtrList<volScalarField> kappaIntFluid(fluidRegions.size());

    ...

        kappaIntFluid.set
        (
            i,
            new volScalarField
            (
                IOobject
                (
                    "kappaInt",
                    runTime.timeName(),
                    fluidRegions[i],
                    IOobject::READ_IF_PRESENT,
                    IOobject::AUTO_WRITE
                ),
                thermoFluid[i].kappa()
            )
        );

chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H

        ...
        
        kappaIntFluid[i].correctBoundaryConditions();
    }

ケースではたとえば次のように熱伝導率を設定する。熱伝達率を h、熱伝導率を k、距離を d とすると、h = k/d なので、k = h*d である (deltaCoeffs() は d の逆数を与える)。

0/air/kappaInt

    air_to_solid
    {
        type            codedFixedValue;
        value           uniform 0;
        redirectType    air_to_solid_kappa;

        code
        #{
            const scalar h = 10.;

            scalarField f(patch().size(), 0.);
            const scalarField& rdelta = patch().deltaCoeffs();

            forAll(f, i)
            {
                f[i] = h/rdelta[i];
            }

            operator==(f);
        #};
    }

kappaInt を用いた場合の壁面熱流束を計算する function object は、以下のようにする。

functions
{
    wallHeatFlux
    {
        type coded;
        redirectType wallHeatFlux;
        region air;

        codeInclude
        #{
            #include "basicThermo.H"
            #include "turbulentFluidThermoModel.H"
            #include "wallFvPatch.H"
        #};

        codeOptions
        #{
            -I$(LIB_SRC)/transportModels/compressible/lnInclude \
            -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
            -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
            -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude
        #};

        codeLibs
        #{
            -lfluidThermophysicalModels \
            -lcompressibleTurbulenceModels
        #};

        codeExecute
        #{
            const basicThermo& thermo =
                mesh().lookupObject<basicThermo>("thermophysicalProperties");

            const compressible::turbulenceModel& turbulence =
                mesh().lookupObject<compressible::turbulenceModel>
                (
                    "turbulenceProperties"
                );

            const volScalarField& h = thermo.he();

            tmp<volScalarField> kappa(thermo.kappa());
            const volScalarField& kappaInt =
                mesh().lookupObject<volScalarField>("kappaInt");
            volScalarField kappaMod(kappa());

            forAll(kappaMod.boundaryField(), patchi)
            {
                forAll(kappaMod.boundaryField()[patchi], i)
                {
                    if (kappaInt.boundaryField()[patchi][i] > 0.)
                    {
                        kappaMod.boundaryFieldRef()[patchi][i]
                            = kappaInt.boundaryField()[patchi][i];
                    }
                }
            }

            surfaceScalarField heatFlux
            (
                fvc::interpolate
                (
                    turbulence.alphaEff()/kappa*kappaMod
                )*fvc::snGrad(h)
            );

            const surfaceScalarField::Boundary& patchHeatFlux =
                heatFlux.boundaryField();

            Info<< "\nWall heat fluxes [W]" << endl;
            forAll(patchHeatFlux, patchi)
            {   
                if (isA<wallFvPatch>(mesh().boundary()[patchi]))
                {   
                    Info<< mesh().boundary()[patchi].name()
                        << " " 
                        << gSum
                            (   
                                mesh().magSf().boundaryField()[patchi]
                                *patchHeatFlux[patchi]
                            )   
                        << endl;
                }   
            }   
            Info<< endl;
        #};

        codeEnd $codeExecute;

        writeControl timeStep;
        writeInterval 1;
    }
}

古い