流体・固体熱連成解析
2016年12月3日
はじめに
流体・固体熱連成解析の手順について。
使用バージョン
OpenFOAM 4.x (blueCFD)
ファイル
- simpleCHT-chtMultiRegionSimpleFoam-4.x-20161203.tar.gz (chtMultiRegionSimpleFoam 用ケースファイル)
- square3dCHT-chtMultiRegionFoam-4.x-20161203.tar.gz (chtMultiRegionFoam 用ケースファイル)
- simpleCHT-kappa-chtMultiRegionSimpleFoam-4.x-20161203.tar.gz (chtMultiRegionSimpleFoam 用ケースファイル、任意熱伝達率)
モデル
次のようなモデルを想定する。
2 次元矩形領域を考え、上を固体、下を流体とする。流体は左から 500 K で流入し、右から出て行く。固体の上端は 300 K 固定とする。流体・固体間の熱伝達を考慮する。
ソルバーには chtMultiRegionSimpleFoam を用いる。
設定手順
- $FOAM_TUTORIALS/heatTransfer/chtMultiRegionSimpleFoam からケースをコピーする。
- 通常通りメッシュを作成し、気体・固体各領域にゾーンを設定する。
ここでは air と solid とする。 - splitMeshRegions により領域を分ける。
$ splitMeshRegions -cellZones -overwrite
system、constant、0 の中に領域ごとのディレクトリができる。 - constant/regionProperties を書き換える。
FoamFile { version 2.0; format ascii; class dictionary; location "constant"; object regionProperties; } regions ( fluid (air) solid (solid) ); - 領域ごとに境界条件を設定する。ここでは直接設定するのではなく、changeDictionary を用いる。system の各領域のディレクトリに changeDictionaryDict を用意する。
air/changeDictionaryDictFoamFile { 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/changeDictionaryDictFoamFile { 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
- constant の中の設定を領域ごとに行う。
- system の中の設定を領域ごとに行う。fvSchemes、fvSolution、decomposeParDict は領域ごとに必要であるが、system 直下にも一応用意しておく必要がある。
- chtMultiRegionSimpleFoam を実行。
$ chtMultiRegionSimpleFoam
- 計算できたら、ポスト処理用に領域ごとの .OpenFOAM を作成する。
$ paraFoam -region air -touch $ paraFoam -region solid -touch
ここでは "simpleCHT{air}.OpenFOAM"、"simpleCHT{solid}.OpenFOAM" といったファイルができる。 - 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;
}
}