OpenFOAMで異方性熱伝導ソルバー anisotropicLaplacianFoam を作る
はじめに
OpenFOAMで固体内の熱伝導を解くとき、熱伝導率を1つのスカラー値として扱うだけなら比較的シンプルです。
しかし、複合材料、積層材、基板、繊維強化材のように「方向によって熱の伝わりやすさが違う」材料を扱いたい場合があります。このような材料では、熱伝導率をスカラーではなくテンソルとして扱う必要があります。
この記事では、内部向け資料である異方性熱伝導を解くための簡単なOpenFOAMソルバー anisotropicLaplacianFoam を紹介・説明します（あくまで身内向けですので、説明が不足しているところはあると思います）。
このソルバーは、静止した固体中の非定常熱伝導だけを扱います。流速、圧力、移流、放射、流体との連成熱伝達は扱いません。その分、異方性熱伝導の実装を追いやすい小さなソルバーになっています。
対象
この記事は、次のような方を想定しています。
- OpenFOAMで独自ソルバーを作ってみたい
-
laplacianFoamに近い構造のソルバーを読んでみたい - 異方性熱伝導率を
symmTensorとして扱いたい - 固体中の発熱と熱流束をOpenFOAMで確認したい
ソルバーが解く式
anisotropicLaplacianFoam が解く式は次の非定常熱伝導方程式です。
rhoCp * dT/dt = div(kappa & grad(T)) + Q
ここで、各量は次の意味です。
| 記号 | OpenFOAM field | 意味 |
|---|---|---|
T |
T |
温度 |
rhoCp |
rhoCp |
体積熱容量 |
kappa |
kappa |
異方性熱伝導率テンソル |
Q |
fvOptions |
任意の体積発熱 |
熱流束は次の符号規約で定義しています。
q = -kappa & grad(T)
grad(T) は温度が高くなる方向を向きます。一方、熱は高温側から低温側へ流れるため、熱流束 q にはマイナス符号が付きます。
異方性熱伝導率とは
等方性材料では、熱伝導率は1つの値で表せます。
k = 1.0
一方、異方性材料では、方向ごとに熱の伝わりやすさが異なります。そのため、熱伝導率をテンソルとして表します。
[ kxx kxy kxz ]
[ kxy kyy kyz ]
[ kxz kyz kzz ]
OpenFOAMの symmTensor は対称テンソルなので、入力ファイルでは6成分で指定します。
(xx xy xz yy yz zz)
例えば、次の指定は
internalField uniform (1.0 0 0 0.3 0 0.3);
次のテンソルを意味します。
[ 1.0 0 0 ]
[ 0 0.3 0 ]
[ 0 0 0.3 ]
この材料は、x 方向には熱が伝わりやすく、y 方向と z 方向には相対的に伝わりにくい材料です。
ソルバーの構成
リポジトリの主なファイル構成は次の通りです。
anisotropicLaplacianFoam.C
createFields.H
write.H
Make/files
Make/options
tutorials/compositePlate2D
それぞれの役割は次の通りです。
| ファイル | 役割 |
|---|---|
anisotropicLaplacianFoam.C |
メインループと行列方程式の組み立て |
createFields.H |
T, kappa, rhoCp, fvOptions の読み込み |
write.H |
gradT と q の出力 |
Make/files |
ビルド対象の指定 |
Make/options |
include path とリンクライブラリの指定 |
tutorials/compositePlate2D |
動作確認用の2次元複合板ケース |
フィールドの読み込み
createFields.H では、温度 T、熱伝導率テンソル kappa、体積熱容量 rhoCp を読み込みます。
volScalarField T
(
IOobject
(
"T",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
kappa は volSymmTensorField として読み込みます。
volSymmTensorField kappa
(
IOobject
(
"kappa",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
rhoCp は volScalarField です。
volScalarField rhoCp
(
IOobject
(
"rhoCp",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
内部発熱は標準の fvOptions で扱います。
#include "createFvOptions.H"
行列方程式
ソルバーの中心は、anisotropicLaplacianFoam.C の次の部分です。
fvScalarMatrix TEqn
(
fvm::ddt(rhoCp, T) - fvm::laplacian(kappa, T)
==
fvOptions(rhoCp, T)
);
左辺の
fvm::ddt(rhoCp, T)
は非定常項です。
rhoCp * dT/dt
に対応します。
次の項が熱伝導項です。
fvm::laplacian(kappa, T)
kappa が volSymmTensorField なので、異方性熱伝導率テンソルを使ったラプラシアンになります。
コード上では
fvm::ddt(rhoCp, T) - fvm::laplacian(kappa, T) == fvOptions(rhoCp, T)
という形で組まれています。これは式を整理すると
rhoCp * dT/dt = div(kappa & grad(T)) + Q
に対応します。
非直交補正
時間ループの内側では、simpleControl による非直交補正ループを使っています。
while (simple.correctNonOrthogonal())
{
fvScalarMatrix TEqn
(
fvm::ddt(rhoCp, T) - fvm::laplacian(kappa, T)
==
fvOptions(rhoCp, T)
);
fvOptions.constrain(TEqn);
TEqn.solve();
fvOptions.correct(T);
}
この構造は、OpenFOAMの標準的な拡散方程式ソルバーに近い形です。メッシュが直交でない場合にも、fvSolution 側の設定に応じて補正を行えます。
出力される派生フィールド
write.H では、出力時刻ごとに温度勾配 gradT と熱流束 q を計算して書き出します。
volVectorField gradT
(
IOobject
(
"gradT",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
fvc::grad(T)
);
熱流束は次のように計算しています。
volVectorField q
(
IOobject
(
"q",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
(-kappa & gradT)
);
異方性材料では、熱流束ベクトルが必ずしも温度等値線に対して単純な向きになるとは限りません。q をベクトル表示すると、材料ごとの熱の流れ方の違いを確認しやすくなります。
必要な入力フィールド
ケースの初期時刻ディレクトリ、通常は 0/ に、次の3つのフィールドが必要です。
| Field | Type | Dimensions | 意味 |
|---|---|---|---|
T |
volScalarField |
[0 0 0 1 0 0 0] |
温度 |
kappa |
volSymmTensorField |
[1 1 -3 -1 0 0 0] |
熱伝導率テンソル |
rhoCp |
volScalarField |
[1 -1 -2 -1 0 0 0] |
体積熱容量 |
2次元ケースでは、frontAndBack などの奥行き方向パッチに empty を指定します。
frontAndBack
{
type empty;
}
ビルド方法
OpenFOAM v2512 の環境を読み込んだ状態で、リポジトリ直下から wmake を実行します。
wmake
実行ファイルは次の場所に作られます。
$(FOAM_USER_APPBIN)/anisotropicLaplacianFoam
ビルド後、次のコマンドでパスが通っているか確認できます。
which anisotropicLaplacianFoam
チュートリアルケース
付属のチュートリアルは tutorials/compositePlate2D にあります。
これは薄い2次元複合板のケースです。
- サイズ:
0.10 m x 0.05 m x 0.001 m - メッシュ:
100 x 50 x 1 - 左境界:
300 K - 右境界:
330 K - 上下境界: 断熱、つまり
zeroGradient - 前後境界: 2次元計算用の
empty
領域は topoSet で作成します。
| Region | 役割 |
|---|---|
materialCore |
中央の材料領域 |
materialInsert |
中央内部の異方性インサート |
heatSource |
内部発熱を与えるセルゾーン |
材料物性は setFields で割り当てます。
| Region |
kappa components (xx xy xz yy yz zz)
|
rhoCp |
|---|---|---|
| Background | (1.0 0 0 0.3 0 0.3) |
2.0e6 |
materialCore |
(8.0 0 0 0.8 0 0.8) |
2.4e6 |
materialInsert |
(0.4 0 0 3.0 0 0.4) |
1.6e6 |
materialCore は x 方向に熱が伝わりやすい材料です。一方、materialInsert は y 方向に熱が伝わりやすい材料として設定されています。
チュートリアルの実行
まずソルバーをビルドします。
wmake
次にチュートリアルケースへ移動します。
cd tutorials/compositePlate2D
Allrun を実行すると、初期ディレクトリの復元、メッシュ生成、領域作成、物性値設定、ソルバー実行までまとめて行います。
./Allrun
手動で実行する場合は、次の流れです。
rm -rf 0
cp -r 0.orig 0
blockMesh
topoSet
setFields
anisotropicLaplacianFoam
計算後、各時刻ディレクトリに次のフィールドが出力されます。
TgradTq
ParaViewなどで T を表示すると温度分布を確認できます。q をベクトル表示すると、材料ごとに熱流束の向きや大きさが変わる様子を確認できます。
内部発熱の設定
チュートリアルでは、fvOptions の scalarSemiImplicitSource を使って内部発熱を与えています。
heatGeneration
{
type scalarSemiImplicitSource;
active yes;
selectionMode cellZone;
cellZone heatSource;
volumeMode specific;
injectionRateSuSp
{
T (1.0e6 0);
}
}
selectionMode cellZone により、heatSource というセルゾーンだけに発熱を与えています。
T (1.0e6 0) は、温度方程式に対して陽的なソース項を与える指定です。このケースでは、局所的な発熱によって温度分布が変わることを確認できます。
このソルバーでできること
anisotropicLaplacianFoam は小さなソルバーですが、次のような検討に使えます。
- 異方性材料中の温度分布の確認
- 複合材料中の熱流束方向の確認
- 局所発熱がある固体部品の簡易解析
- 独自熱伝導ソルバー作成の学習
-
volSymmTensorFieldを使った拡散項の実装例
特に、kappa をスカラーではなく volSymmTensorField として読み込み、
fvm::laplacian(kappa, T)
に渡している点が重要です。
注意点
このソルバーは意図的にシンプルにしてあります。次の物理モデルは含んでいません。
- 流体の移流
- 圧力方程式
- 速度場
- 放射
- 相変化
- 複数領域の共役熱伝達
- 温度依存物性
そのため、流体と固体の連成や、温度によって物性が変化する問題を扱うには拡張が必要です。
また、異方性熱伝導率テンソルは物理的に妥当な値にする必要があります。少なくとも、熱伝導率テンソルが非物理的な向きや符号になっていないか確認してください。
よくあるエラー
anisotropicLaplacianFoam: command not found
ソルバーがビルドされていないか、OpenFOAM環境が読み込まれていない可能性があります。
wmake
which anisotropicLaplacianFoam
を確認します。
kappa または rhoCp が見つからない
0/ ディレクトリに必要なフィールドがない可能性があります。
チュートリアルでは、次のように 0.orig から復元します。
rm -rf 0
cp -r 0.orig 0
2次元ケースで patch type のエラーが出る
T, kappa, rhoCp のすべてで、前後方向パッチが empty になっているか確認します。
frontAndBack
{
type empty;
}
発熱の効果が見えない
topoSet を実行して heatSource の cellZone が作られているか確認します。
topoSet
また、constant/polyMesh/cellZones に heatSource が存在するか確認します。
まとめ
ポイントは次の3つです。
- 熱伝導率
kappaをvolSymmTensorFieldとして読み込む -
fvm::laplacian(kappa, T)でテンソル熱伝導率を使った拡散項を組む -
q = -kappa & grad(T)を出力し、異方性材料中の熱流束を確認する
最小限の構成なので、独自ソルバーの学習用としても扱いやすいと思います。まずは付属の tutorials/compositePlate2D を動かし、T と q を可視化して、材料ごとの熱の流れ方の違いを確認すると理解しやすいです。