OpenFOAMで総括熱伝達率（U値）を計算する方法。chtMultiRegionFoam対応
changeDictionaryDictによる境界条件の設定
system/fluid_1/changeDictionaryDict（system/fluid_2/changeDictionaryDictも同様）
k
{
internalField uniform 0.00375;
boundaryField
{
".*"
{
type kqRWallFunction;
value $internalField;
}
defaultFaces
{
type empty;
}
}
}
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;
}
defaultFaces
{
type empty;
}
}
}
nut
{
internalField uniform 6000;
boundaryField
{
".*"
{
type nutkWallFunction;
value uniform 0;
}
defaultFaces
{
type empty;
}
}
}
U
{
internalField uniform ( 0 0 0 );
boundaryField
{
".*"
{
type fixedValue;
value uniform (0 0 0);
}
defaultFaces
{
type empty;
}
}
}
T
{
internalField uniform 293.15;
boundaryField
{
".*"
{
type fixedValue;
value uniform 293.15;
}
defaultFaces
{
type empty;
}
".*_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;
}
defaultFaces
{
type empty;
}
}
}
p
{
internalField uniform 100000;
boundaryField
{
".*"
{
type calculated;
value uniform 100000;
}
defaultFaces
{
type empty;
}
}
}
alphat
{
internalField uniform 0;
boundaryField
{
".*"
{
type compressible::alphatWallFunction;
Prt 0.85;
value uniform 0;
}
defaultFaces
{
type empty;
}
}
}
G
{
internalField uniform 0;
boundaryField
{
".*"
{
type MarshakRadiation;
value uniform 0;
}
defaultFaces
{
type empty;
}
}
}
qr
{
internalField uniform 0;
boundaryField
{
".*"
{
type calculated;
value uniform 0;
}
defaultFaces
{
type empty;
}
}
}
IDefault
{
internalField uniform 0;
boundaryField
{
".*"
{
type greyDiffusiveRadiation;
value uniform 0;
}
defaultFaces
{
type empty;
}
}
}
// ************************************************************************* //
system/solid_1/changeDictionaryDict（system/solid_2/changeDictionaryDictも同様）
T
{
internalField uniform 293.15;
boundaryField
{
".*"
{
type fixedValue;
value uniform 293.15;
}
defaultFaces
{
type empty;
}
".*_to_.*"
{
type compressible::turbulentTemperatureRadCoupledMixed;
value uniform 293.15;
Tnbr T;
kappaMethod fluidThermo;
qrNbr none;
qr none;//qr;
thermalInertia true;
kappa none;
}
}
}
p
{
internalField uniform 100000;
boundaryField
{
".*"
{
type calculated;
value uniform 100000;
}
defaultFaces
{
type empty;
}
}
}
領域の指定と境界条件の設定
次に領域を指定します。
constant/regionProperties
regions
(
fluid (fluid_1 fluid_2)
solid (solid_1 solid_2)
);
領域を正しく設定しておかないと以下のコマンドの$regionでエラーが出ます。
for region in $(foamListRegions); do changeDictionary -region "$region"; done
スクリプトで書くと以下となります。
#!/bin/bash
for region in $(foamListRegions)
do
echo "=== changeDictionary for region: $region ==="
changeDictionary -region "$region"
done
実行すると以下のようになります。
0/fluid_1/U
dimensions [ 0 1 -1 0 0 0 0 ];
internalField uniform ( 0 0 0 );
boundaryField
{
f1_minX
{
type fixedValue;
value uniform ( 0 0 0 );
}
f1_minY
{
type fixedValue;
value uniform ( 0 0 0 );
}
f1_maxY
{
type fixedValue;
value uniform ( 0 0 0 );
}
defaultFaces
{
type empty;
}
fluid_1_to_solid_1
{
type fixedValue;
value uniform ( 0 0 0 );
}
}
0/fluid_1/T
dimensions [ 0 0 0 1 0 0 0 ];
internalField uniform 293.15;
boundaryField
{
f1_minX
{
type fixedValue;
value uniform 333.15;
}
f1_minY
{
type fixedValue;
value uniform 293.15;
}
f1_maxY
{
type fixedValue;
value uniform 293.15;
}
defaultFaces
{
type empty;
}
fluid_1_to_solid_1
{
type compressible::turbulentTemperatureRadCoupledMixed;
value uniform 293.15;
Tnbr T;
kappaMethod fluidThermo;
qrNbr none;
qr none;
thermalInertia true;
kappa none;
}
}
0/fluid_1/p_rgh
dimensions [ 1 -1 -2 0 0 0 0 ];
internalField uniform 100000;
boundaryField
{
f1_minX
{
type fixedFluxPressure;
value uniform 100000;
}
f1_minY
{
type fixedFluxPressure;
value uniform 100000;
}
f1_maxY
{
type fixedFluxPressure;
value uniform 100000;
}
defaultFaces
{
type empty;
}
fluid_1_to_solid_1
{
type fixedFluxPressure;
value uniform 100000;
}
}
物性値
流体領域の物性値は空気にしておきます。
constant/fluid_1/thermophysicalProperties（constant/fluid_1/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;
}
}
固体領域の物性値は以下とします。
今回は、デフォルトで仮想の材料で鉄鋼材っぽい値が設定されているので、そのまま使います。
constant/solid_1/thermophysicalProperties（constant/solid_2/thermophysicalPropertiesも同様）
thermoType
{
type heSolidThermo;
mixture pureMixture;
transport constIso;
thermo hConst;
equationOfState rhoConst;
specie specie;
energy sensibleEnthalpy;
}
mixture
{
specie
{
molWeight 50;
}
transport
{
kappa 80;
}
thermodynamics
{
Hf 0;
Cp 450;
}
equationOfState
{
rho 8000;
}
}
乱流モデルの設定
乱流モデルの設定は流体領域だけに行います。
constant/fluid_1/turbulenceProperties（constant/fluid_2/turbulencePropertiesも同様）
simulationType laminar;
今回は乱流モデルは使用しません。
計算の制御
system/controlDict
application chtMultiRegionFoam;
startFrom latestTime;
startTime 0.001;
stopAt endTime;
endTime 3600;
deltaT 0.001;
writeControl adjustable;
writeInterval 100;
purgeWrite 0;
writeFormat ascii;
writePrecision 8;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable yes;
maxCo 10.0;
// Maximum diffusion number
maxDi 10.0;
adjustTimeStep yes;
functions
{
// #include "vtkWrite"
wallHeatFlux_fluid_1
{
type wallHeatFlux;
libs ("libfieldFunctionObjects.so");
region fluid_1;
patches (fluid_1_to_solid_1);
writeControl writeTime;
}
wallHeatFlux_solid_1
{
type wallHeatFlux;
libs ("libfieldFunctionObjects.so");
region solid_1;
patches (solid_1_to_fluid_1);
writeControl writeTime;
}
qAvg_fluid_1
{
type surfaceFieldValue;
libs ("libfieldFunctionObjects.so");
region fluid_1;
writeControl writeTime;
writeFields false; // 面上の場を書き出さず、集計値だけ出す
regionType patch;
name fluid_1_to_solid_1;
operation areaAverage;
fields (wallHeatFlux);
}
qInt_fluid_1
{
type surfaceFieldValue;
libs ("libfieldFunctionObjects.so");
region fluid_1;
writeControl writeTime;
writeFields false; // 面上の場を書き出さず、集計値だけ出す
regionType patch;
name fluid_1_to_solid_1;
operation areaIntegrate;
fields (wallHeatFlux);
}
qAvg_solid_1
{
type surfaceFieldValue;
libs ("libfieldFunctionObjects.so");
region solid_1;
writeControl writeTime;
writeFields false; // 面上の場を書き出さず、集計値だけ出す
regionType patch;
name solid_1_to_fluid_1;
operation areaAverage;
fields (wallHeatFlux);
}
qInt_solid_1
{
type surfaceFieldValue;
libs ("libfieldFunctionObjects.so");
region solid_1;
writeControl writeTime;
writeFields false; // 面上の場を書き出さず、集計値だけ出す
regionType patch;
name solid_1_to_fluid_1;
operation areaIntegrate;
fields (wallHeatFlux);
}
}
今回最大クーラン数を流体側も固体側も10にしました。
計算実行
2次元計算でメッシュ数も多くないので並列計算させずにいきます。
chtMultiRegionFoam > log.chtMultiRegionFoam &
今回最大クーラン数を流体側も固体側も10にしています。
ログを見ると最大クーラン数は固体側で10になっており、それによって時間刻みが決められているようです。
Region: fluid_1 Courant Number mean: 0.91343604 max: 4.4381642 Region: fluid_2 Courant Number mean: 0.071744007 max: 0.18203249 Region: solid_1 Diffusion Number mean: 1.2466978 max: 9.9983571 Region: solid_2 Diffusion Number mean: 1.2466978 max: 9.9983571 deltaT = 0.11248152 Time = 50.3957
陰解法では最大クーラン数が 1 を超えても計算が安定する場合がありますが、これは過渡解の精度を保証するものではありません。
特に一次精度の時間離散では時間刻みの増大に伴って数値拡散が強まり、温度応答が平滑化される可能性があります。
- 温度の立ち上がりを議論したい
- PID制御や時定数を見たい
- 短時間の熱応答を比較したい
- 流れ場と温度場の瞬時変動が重要
- 実験の「時系列」まで合わせたい
このような場合に、クーラン数が大きすぎると精度的に良くない場合もあります。
そのため、最大クーラン数を大きく設定する場合は、目的量に対する時間刻み依存性を別途確認する必要があります。
setMultiRegionDeltaT.Hに最大クーラン数に対するアルゴリズムの記載があります。
if (adjustTimeStep)
{
if (CoNum == -GREAT)
{
CoNum = SMALL;
}
if (DiNum == -GREAT)
{
DiNum = SMALL;
}
scalar maxDeltaTFluid = maxCo/(CoNum + SMALL);
scalar maxDeltaTSolid = maxDi/(DiNum + SMALL);
scalar deltaTFluid =
min
(
min(maxDeltaTFluid, 1.0 + 0.1*maxDeltaTFluid),
1.2
);
runTime.setDeltaT
(
min
(
min(deltaTFluid, maxDeltaTSolid)*runTime.deltaTValue(),
maxDeltaT
)
);
Info<< "deltaT = " << runTime.deltaTValue() << endl;
}
adjustTimeStep yesのときに動く- 流体側は最大クーラン数
CoNumを見る - 固体側は最大拡散数
DiNumを見る - それぞれから「時間刻みを何倍できるか」を計算する
- 流体と固体のうち厳しい方に合わせて
deltaTを決める - ただし時間刻みは一気に増やしすぎないよう制限する
- 最後に
maxDeltaTを超えないようにする
結果
U値の計算
postProcessing/fluid_1/qAvg_fluid_1/0/surfaceFieldValue.dat
# Region type : patch fluid_1_to_solid_1 # Faces : 60 # Area : 3.00000000e-03 # Scale factor : 1.00000000e+00 # Time areaAverage(wallHeatFlux) 100 -1.36060750e+01 200 -1.35502813e+01 300 -1.35035271e+01 400 -1.34654931e+01 ...（省略）... 3200 -1.33065801e+01 3300 -1.33065652e+01 3400 -1.33065784e+01 3500 -1.33065749e+01 3600 -1.33065780e+01
# Region type : patch solid_1_to_fluid_1 # Faces : 60 # Area : 3.00000000e-03 # Scale factor : 1.00000000e+00 # Time areaAverage(wallHeatFlux) 100 1.28302150e+01 200 1.29376556e+01 300 1.30110614e+01 400 1.30683247e+01 ...（省略）... 3200 1.33065801e+01 3300 1.33065653e+01 3400 1.33065782e+01 3500 1.33065747e+01 3600 1.33065778e+01
ここから定常に落ち着いたであろう最終時間でのU値を求めます。
1 2
Follow @t_kun_kamakiri
