T_distrib.py 技術ドキュメント

プログラムの動作

T_distrib.py は、積層型熱電 (TE) スタックの1次元定常状態温度プロファイルを計算するPythonプログラムです。このプログラムは、スタックを構成する複数の層（電極、TE半導体、界面）を考慮し、外部環境との熱交換モデルに基づいて温度分布を予測します。

主な機能は以下の2つのモードで提供されます。

  1. delta モード: このモードでは、スタックの上面と下面における表面熱抵抗を、明示的な空気境界層（厚さ delta_top / delta_bottom、熱伝導率 k_air）としてモデル化します。これは、単純な熱抵抗ネットワークとして扱われ、総熱抵抗とそこを流れる熱流束を計算します。

  2. conv_rad モード: このモードでは、表面熱抵抗を自然対流と熱放射の物理モデルに基づいて予測します。

    • 自然対流: 水平平板のヌセルト数相関式を用いて、対流熱伝達係数を計算します。

    • 熱放射: ステファン・ボルツマンの法則に基づいて放射熱流束を計算します。 このモードでは、内部スタックと外部環境（空気、放射環境）との間で、熱流束が平衡するような表面温度を数値的に求めます。

両モード共通で、内部スタックは以下の層で構成されます。 上部電極 / 界面1/2 / TE半導体 / 界面2/3 / 下部電極

プログラムは、以下の情報を標準出力に出力します。

  • 各層および界面の熱抵抗

  • 総熱抵抗（または内部スタック熱抵抗および表面実効熱抵抗）

  • 各界面の温度

  • TE半導体層をまたぐ温度降下

また、温度プロファイルを示す図をPNG画像ファイルとして保存します。この図は、物理スケールと模式（変形）スケールの2つのグラフが縦に並べられた形式で表示されます。

現在のプログラムのバグについて: main 関数内の delta モードの処理が2度記述されており、2度目の呼び出しにおいて compute_piecewise_profile 関数の戻り値の数が誤って指定されています。これにより、delta モードでプログラムを実行すると ValueError: too many values to unpack エラーが発生し、正常に動作しません。この技術ドキュメントでは、バグが修正された場合の delta モードの期待される動作について記述します。

原理

1次元定常熱伝導

プログラムは、1次元定常状態の熱伝導を仮定しています。この条件下では、熱流束 \(q\) はスタック全体で一定であり、フーリエの法則によって各層の温度勾配と熱伝導率に関係付けられます。

\[q = -k \frac{dT}{dz}\]

ここで、\(k\) は熱伝導率、\(dT/dz\) は温度勾配です。各層の熱抵抗 \(R\) は、厚さ \(\Delta z\) と熱伝導率 \(k\) を用いて次のように定義されます。

\[R = \frac{\Delta z}{k}\]

層を通過する熱流束 \(q\) と温度差 \(\Delta T\) の関係は、オームの法則のアナロジーとして次のように表されます。

\[q = \frac{\Delta T}{R}\]

delta モードの原理

このモードでは、スタックの各層および上面・下面の空気境界層を独立した熱抵抗として扱います。全体の熱抵抗 \(R_{total}\) は、これらすべての層の熱抵抗の合計として計算されます。

\[R_{total} = R_{top\_air} + R_{top\_electrode} + R_{interface12} + R_{TE} + R_{interface23} + R_{bottom\_electrode} + R_{bottom\_air}\]

熱流束 \(q\) は、上面空気温度 \(T_h\) と下面空気温度 \(T_c\) の差を総熱抵抗で割ることで直接計算されます。

\[q = \frac{T_h - T_c}{R_{total}}\]

この熱流束を用いて、各層を通過する温度降下を順次計算し、スタック全体の温度プロファイルを決定します。

conv_rad モードの原理

このモードでは、スタック内部の熱伝導に加え、上面と下面の熱交換を自然対流と熱放射によってモデル化します。スタックの表面温度 \(T_{s,top}\) と \(T_{s,bottom}\) は未知数となり、スタックを貫流する熱流束と表面からの熱交換が平衡するという条件を数値的に解きます。

自然対流: 水平平板における自然対流の熱伝達係数 \(h_{nat}\) は、無次元数であるレイリー数 \(Ra\) とヌセルト数 \(Nu\) を用いた経験的な相関式から導出されます。

レイリー数 \(Ra\) は次式で計算されます。

\[Ra = \frac{g \beta |\Delta T| L_{char}^3}{\nu \alpha}\]

ここで、\(g\) は重力加速度、\(|\Delta T|\) は表面と空気の温度差、\(L_{char}\) は代表長さ、\(\nu\) は動粘度、\(\alpha\) は熱拡散率、\(\beta\) は体積膨張係数 (\(1/T_{film,K}\)、ここで \(T_{film,K}\) は膜温度の絶対温度) です。

ヌセルト数 \(Nu\) は、安定性（上向き加熱、下向き冷却などで不安定）に応じて以下の相関式で計算されます。

  • 不安定な場合:

    • \(Ra < 10^7\): $\(Nu = 0.54 Ra^{0.25}\)$

    • \(Ra \ge 10^7\): $\(Nu = 0.15 Ra^{1/3}\)$

  • 安定な場合: $\(Nu = 0.27 Ra^{0.25}\)$

熱伝達係数 \(h_{nat}\) は \(Nu = h_{nat} L_{char} / k_{air}\) の関係から得られます。

熱放射: 表面からの放射熱流束 \(q_{rad}\) は、ステファン・ボルツマンの法則に基づいて計算されます。

\[q_{rad} = \epsilon \sigma (T_{surface,K}^4 - T_{ambient,K}^4)\]

ここで、\(\epsilon\) は表面の放射率、\(\sigma\) はステファン・ボルツマン定数 (\(5.670374419 \times 10^{-8} \text{ W/(m}^2 \text{K}^4)\))、\(T_{surface,K}\) と \(T_{ambient,K}\) はそれぞれ表面と周囲の絶対温度です。

熱流束の平衡: conv_rad モードでは、以下の3つの熱流束が平衡するようなスタック表面温度 \(T_{s,top}\) および \(T_{s,bottom}\) を求めます。

  1. 上部空気環境から上部スタック表面への全熱流束 \(q_{top}\) （対流 + 放射）

  2. 下部スタック表面から下部空気環境への全熱流束 \(q_{bottom}\) （対流 + 放射）

  3. 上部スタック表面から下部スタック表面を介して内部スタックを流れる熱流束 \(q_{internal}\)

これらは全て等しくなるはずです。すなわち、\(q_{top} = q_{internal} = q_{bottom}\)。 \(q_{internal} = (T_{s,top} - T_{s,bottom}) / R_{internal}\) であり、\(R_{internal}\) は内部スタックの総熱抵抗です。 このシステムは、上部表面温度 \(T_{s,top}\) を変数とする残差関数 residual を定義し、二分法を用いて \(q_{top} - q_{bottom} = 0\) を満たす \(T_{s,top}\) を探索することで解かれます。

空気物性値

空気の熱伝導率 \(k_{air}\)、動粘度 \(\nu\)、熱拡散率 \(\alpha\)、およびプラントル数 \(Pr\) は、air_properties_simple 関数によって簡単な温度依存モデルで計算されます。 例: \(k_{air} = k_{air,ref} (T_{film,K} / 300.0)^{0.10}\)

必要な非標準ライブラリとインストール方法

このプログラムの実行には、以下の非標準Pythonライブラリが必要です。

  • numpy: 数値計算（特に配列操作）に使用されます。

  • matplotlib: グラフ描画に使用されます。

これらのライブラリは、Pythonのパッケージインストーラ pip を使ってインストールできます。コマンドプロンプトやターミナルで以下のコマンドを実行してください。

pip install numpy matplotlib

必要な入力ファイル

T_distrib.py は、実行に必要な入力ファイルを必要としません。 すべてのパラメータは、コマンドライン引数として直接指定されます。

生成される出力ファイル

プログラムは、スタック内の温度プロファイルを示す画像を生成し、以下の形式で保存します。

  • ファイル名:

    • デフォルト: layered_temperature_profile.png

    • --save 引数で任意のファイル名を指定できます。

  • ファイル形式: PNG画像ファイル

  • 内容: 生成される画像には、2つの垂直方向に積み重ねられたプロットが含まれます。

    1. 物理スケールの温度プロファイル: スタックの物理的な厚さに対する温度分布が示されます。各層の境界は破線で示され、層の名称がプロット内に表示されます。

    2. 模式（変形）スケールの温度プロファイル: 各層の熱抵抗や表面抵抗を強調するために、層の厚さを任意に調整した模式的なスケールでの温度分布が示されます。これにより、抵抗の大きな層での温度降下が視覚的に分かりやすくなります。各領域の境界は破線で示され、領域の名称がプロット内に表示されます。

コマンドラインでの使用例 (Usage)

T_distrib.py の基本的なコマンドラインでの使用法は以下の通りです。 --help オプションを使用して、利用可能なすべての引数とその説明を確認できます。

python T_distrib.py --help

上記コマンドの出力例:

usage: T_distrib.py [-h] [--surface_mode {delta,conv_rad}] [--Th TH] [--Tc TC] [--k_air K_AIR] [--delta_top DELTA_TOP] [--delta_bottom DELTA_BOTTOM] [--d1 D1] [--k1 K1] [--delta12 DELTA12] [--k12 K12] [--d2 D2] [--k2 K2] [--k23 K23] [--delta23 DELTA23] [--d3 D3] [--k3 K3] [--L_top L_TOP] [--L_bottom L_BOTTOM] [--eps_top EPS_TOP] [--eps_bottom EPS_BOTTOM] [--T_rad_top T_RAD_TOP] [--T_rad_bottom T_RAD_BOTTOM] [--save SAVE] [--npts NPTS] [--schematic_layer_width SCHEMATIC_LAYER_WIDTH] [--schematic_interface_width SCHEMATIC_INTERFACE_WIDTH] [--schematic_surface_width SCHEMATIC_SURFACE_WIDTH]

Steady-state 1D temperature profile of layered TE stack.

options:
  -h, --help            show this help message and exit
  --surface_mode {delta,conv_rad}
                        Surface model: explicit delta-layers or natural convection + radiation.
  --Th TH               Top air temperature [°C]
  --Tc TC               Bottom air temperature [°C]
  --k_air K_AIR         Reference thermal conductivity of air [W/mK]
  --delta_top DELTA_TOP
                        Top air boundary layer thickness [m] (delta mode)
  --delta_bottom DELTA_BOTTOM
                        Bottom air boundary layer thickness [m] (delta mode)
  --d1 D1               Top electrode thickness [m]
  --k1 K1               Top electrode thermal conductivity [W/mK]
  --delta12 DELTA12     Interface 1/2 equivalent thickness [m]
  --k12 K12             Interface 1/2 effective thermal conductivity [W/mK]
  --d2 D2               TE semiconductor thickness [m]
  --k2 K2               TE semiconductor thermal conductivity [W/mK]
  --k23 K23             Interface 2/3 effective thermal conductivity [W/mK]
  --delta23 DELTA23     Interface 2/3 equivalent thickness [m]
  --d3 D3               Bottom electrode thickness [m]
  --k3 K3               Bottom electrode thermal conductivity [W/mK]
  --L_top L_TOP         Characteristic length for top natural convection [m]
  --L_bottom L_BOTTOM
                        Characteristic length for bottom natural convection [m]
  --eps_top EPS_TOP     Top surface emissivity
  --eps_bottom EPS_BOTTOM
                        Bottom surface emissivity
  --T_rad_top T_RAD_TOP
                        Top radiation temperature [°C]
  --T_rad_bottom T_RAD_BOTTOM
                        Bottom radiation temperature [°C]
  --save SAVE           Output PNG filename
  --npts NPTS           Plot points per explicit region
  --schematic_layer_width SCHEMATIC_LAYER_WIDTH
                        Display width for normal layers
  --schematic_interface_width SCHEMATIC_INTERFACE_WIDTH
                        Display width for interface layers
  --schematic_surface_width SCHEMATIC_SURFACE_WIDTH
                        Display width for surface regions

コマンドラインでの具体的な使用例

例1: conv_rad モードでの熱電スタックの温度プロファイル計算

この例では、conv_rad モードを使用して、自然対流と熱放射を考慮したスタックの温度プロファイルを計算します。上面の空気温度を100°C、下面を20°Cとし、トップ表面の放射環境温度を50°Cに設定して、熱放射の影響を明確にします。

python T_distrib.py --surface_mode conv_rad --Th 100 --Tc 20 --T_rad_top 50 --save conv_rad_example.png

実行結果の説明: このコマンドを実行すると、以下のような情報が標準出力に表示され、conv_rad_example.png という画像ファイルが生成されます。

=== Mode ===
surface_mode = conv_rad

=== Input temperatures ===
Top air temperature    Th = 100 °C
Bottom air temperature Tc = 20 °C
Top radiation temperature    = 50 °C
Bottom radiation temperature = 20 °C

=== Natural-convection / radiation parameters ===
L_top    = 0.01 m
L_bottom = 0.01 m
eps_top    = 0.8
eps_bottom = 0.8

=== Internal regions and thermal resistances (top -> bottom) ===
Region               type         thickness [m]   k [W/mK]    R=t/k [m^2K/W]
Top electrode        layer          2.000000e-07    100.000000  2.000000e-09
Interface 1/2        interface      5.000000e-08      0.300000  1.666667e-07
TE semiconductor     layer          2.000000e-06      0.500000  4.000000e-06
Interface 2/3        interface      5.000000e-08      0.300000  1.666667e-07
Bottom electrode     layer          2.000000e-07    100.000000  2.000000e-09

Internal thermal resistance R_internal = 4.335333e-06 m^2 K / W
Effective top surface resistance       = 1.156828e-01 m^2 K / W
Effective bottom surface resistance    = 1.258797e-01 m^2 K / W
Total effective resistance             = 2.415617e-01 m^2 K / W
Heat flux q                            = 3.311746e+02 W / m^2

=== Surface heat-flux decomposition ===
Top    convection q_conv_top  = 1.135649e+02 W / m^2
Top    radiation  q_rad_top   = 2.176097e+02 W / m^2
Bottom convection q_conv_bot  = 2.872295e+02 W / m^2
Bottom radiation  q_rad_bot   = 4.394503e+01 W / m^2

=== Surface / interface temperatures ===
Location                   Temperature [°C]
Top air                             100.000000
Top sample surface                   61.737199
After Top electrode                  61.736537
After Interface 1/2                  61.681177
After TE semiconductor               60.370502
After Interface 2/3                  60.315142
Bottom sample surface                60.314480
Bottom air                           20.000000

=== TE layer summary ===
T at TE top    = 61.681177 °C
T at TE bottom = 60.315142 °C
Delta T across TE layer = 1.366035 K

Saved plot: conv_rad_example.png

出力から、上部空気温度 (100°C) よりも上部サンプル表面温度 (約61.7°C) が低いことが分かります。これは、熱放射環境温度を50°Cに設定したことで、上部表面から放射による放熱も起こり、熱流束が減少しているためです。また、各層の抵抗と温度降下、TE半導体層にかかる温度差も計算されています。conv_rad_example.png には、これらの温度プロファイルが可視化されます。

例2: delta モードでの熱電スタックの温度プロファイル計算

注意: 現在の T_distrib.py のコードにはバグが存在するため、以下のコマンドを実行しても ValueError: too many values to unpack エラーが発生し、正常に動作しません。このドキュメントでは、バグが修正された場合の期待される動作を記述します。

この例では、delta モードを使用し、明示的な空気境界層を設定してスタックの温度プロファイルを計算します。上面の空気温度を100°C、下面を20°Cとします。

python T_distrib.py --surface_mode delta --Th 100 --Tc 20 --delta_top 1e-4 --delta_bottom 1e-4 --save delta_example.png

期待される実行結果の説明 (バグ修正後の動作): このコマンドが正常に実行されると、以下のような情報が標準出力に表示され、delta_example.png という画像ファイルが生成されます。

=== Mode ===
surface_mode = delta

=== Input temperatures ===
Top air temperature    Th = 100.0000 °C
Bottom air temperature Tc = 20.0000 °C

=== Regions and thermal resistances (top -> bottom) ===
Region               type         thickness [m]   k [W/mK]    R=t/k [m^2K/W]
Top air boundary     surface        1.000000e-04      0.026000  3.846154e-03
Top electrode        layer          2.000000e-07    100.000000  2.000000e-09
Interface 1/2        interface      5.000000e-08      0.300000  1.666667e-07
TE semiconductor     layer          2.000000e-06      0.500000  4.000000e-06
Interface 2/3        interface      5.000000e-08      0.300000  1.666667e-07
Bottom electrode     layer          2.000000e-07    100.000000  2.000000e-09
Bottom air boundary  surface        1.000000e-04      0.026000  3.846154e-03

Total thermal resistance R_total = 7.863641e-03 m^2 K / W
Heat flux q = (Th - Tc)/R_total = 1.017336e+04 W / m^2

=== Interface temperatures ===
Location                   Temperature [°C]
Top air side                        100.000000
After Top air boundary               60.840733
After Top electrode                  60.840712
After Interface 1/2                  60.756208
After TE semiconductor               52.610515
After Interface 2/3                  52.526011
After Bottom electrode               52.525990
After Bottom air boundary            20.000000

=== TE layer summary ===
T at TE top    = 60.756208 °C
T at TE bottom = 52.526011 °C
Delta T across TE layer = 8.230197 K

Saved plot: delta_example.png

出力から、上下の空気境界層が最も大きな熱抵抗を提供し、その層で最も大きな温度降下が発生していることが分かります。熱流束は単純な総熱抵抗から計算され、各界面温度もそれに従って決定されます。delta_example.png には、これらの温度プロファイルが可視化されます。