円筒座標系の熱伝導
円筒座標系の熱伝導の理論基礎
円筒座標系の熱伝導方程式
先生、配管やシャフトの熱解析はデカルト座標だとうまくいかないんですか？
円筒形状ではデカルト座標だと形状表現に無駄が多く、メッシュが複雑になる。円筒座標 $(r, \phi, z)$ を使えば対称性を活かして次元を落とせる。1次元径方向の定常熱伝導方程式は
デカルト座標のフーリエの法則と比べると $1/r$ の項が増えてますね。
体積要素が $r$ に比例するため、面積の変化を反映した項だ。これが臨界断熱半径の物理的原因でもある。
解析解（内部発熱なし）
$\dot{q}_v = 0$ の場合の解はこうなる。
温度分布は対数関数で、平板の線形分布とは異なる。熱流量は
$\ln$ が出るのは面積が変化するからですね。
そう。断面積 $A = 2\pi r L$ が $r$ に比例するため、一定の熱流量 $q$ に対して熱流束 $q'' = q/A$ は $r$ に反比例する。内面の方が熱流束が大きい。
内部発熱がある場合
一様な内部発熱 $\dot{q}_v$ がある場合の解は
中心 $r=0$ で最高温度 $T_{\max} = T_s + \dot{q}_v R^2/(4k)$ となる。
電線のジュール発熱がこのケースですね。
その通り。AWG18の銅線（直径1.02mm）に10Aを流すと約10 W/m の発熱。$\dot{q}_v \approx 1.2 \times 10^7$ W/m$^3$ で、この式から中心温度上昇を見積もれる。
円筒座標Laplacianの由来
円筒座標での熱伝導方程式は∂²T/∂r² + (1/r)∂T/∂r + (1/r²)∂²T/∂θ² + ∂²T/∂z² + q̇/λ = 0と表される。このLaplacianの(1/r)∂T/∂r項は座標変換のヤコビアンに起源を持つ。Lameが1833年に弾性論で円筒座標Laplacianを整備し、その20年後にFourier熱方程式へ適用された歴史的経緯がある。
円筒座標系の熱伝導の数値計算手法
軸対称FEMモデル
円筒の熱伝導をFEMで解く場合、3Dモデルにする必要がありますか？
周方向に対称なら2D軸対称モデルで十分だ。計算コストが1/100以下になる。
| モデルタイプ | 自由度数 | 計算時間 | 精度 |
|---|---|---|---|
| 3D全体 | 〜100万 | 分オーダー | 基準 |
| 3D 1/4対称 | 〜25万 | 秒オーダー | 同等 |
| 2D軸対称 | 〜1000 | 瞬時 | 同等 |
2D軸対称で十分なのに3Dで回すのは無駄ですね。
対称性の活用は計算効率の基本だ。Ansys MechanicalではPLANE55（KEYOPT(3)=1）、AbaqusではDCAX4（軸対称4節点要素）を使う。
有限差分法での離散化
円筒座標の中心差分は注意が必要だ。
$r = 0$（中心軸）では$1/r$が特異点になるため、ロピタルの定理を適用して
とする。
中心の特異点は別扱いが必要なんですね。
FEMでも同様に中心軸の処理が重要だ。軸対称要素では $r = 0$ の節点で自動的に正しい処理が行われるが、3Dモデルで軸上にウェッジ要素を使う場合は注意が必要だ。
多層円筒の解法
多層円筒（配管+断熱材+被覆など）は各層で解析解を接続する。
層間の接触熱抵抗を含める場合は分母に $1/(2\pi r_i h_{c,i} L)$ を追加する。
電気回路の直列抵抗と同じ足し算ですね。
まさに熱抵抗ネットワークの考え方だ。5層以上の多層でも手計算可能な点が円筒座標の強みだ。
対数温度分布の解法ステップ
無限長中空円筒（内半径ri、外半径ro）の定常温度分布はT(r)=Ti + (Ti−To)·ln(r/ro)/ln(ri/ro)で、対数型になる。この解法は変数分離→ODE積分→2境界条件でCを決定の3ステップで完結する。日本の大学入試問題にも類題が出題されており、国内の熱工学テキストでは必ず章頭例題として扱われる。
円筒座標系の熱伝導の実務適用
配管の放熱計算
実務で円筒熱伝導が使われる典型的な場面を教えてください。
最も多いのはプラント配管の熱損失計算だ。蒸気配管（外径114.3mm、SUS304、断熱材50mm）の例を示す。
| 層 | $r$ [mm] | $k$ [W/(m K)] | 熱抵抗 [m K/W] |
|---|---|---|---|
| 配管壁 | 48.6→57.15 | 16.3 | 0.00160 |
| 断熱材 | 57.15→107.15 | 0.05 | 2.016 |
| 内面対流 | — | — | 0.00033 |
| 外面対流 | — | — | 0.0149 |
全熱抵抗 $R_{\text{total}} = 2.033$ m K/W。断熱材が全体の99%を占める。
配管壁の熱抵抗はほぼ無視できるレベルですね。
そう。金属管の熱抵抗は断熱材の1/1000以下で、実質的に温度降下はない。だから配管材質が銅でもSUSでも放熱量にほとんど影響しない。
ログ平均半径
円筒の単位長さあたり熱抵抗 $R = \ln(r_2/r_1)/(2\pi k)$ は、等価平板の $R = t/(k \cdot A_{lm})$ と表現できる。$A_{lm}$ はログ平均面積で
対数平均を使えば平板の式が流用できるんですね。
$r_2/r_1 < 2$ なら算術平均 $(r_1+r_2)/2$ でも誤差4%以下だ。簡易計算では算術平均で十分なことが多い。
メッシュ作成の注意点
円筒の軸対称メッシュで注意すべき点を挙げる。
- 径方向にバイアスメッシュ（内面側を細かく）を使う
- 薄肉管は径方向に最低3要素
- 多層構造では層間でノードを共有させるか、Tied Contactを使用
バイアスメッシュは内面の方が熱流束が大きいからですね。
その通り。$q'' \propto 1/r$ だから内面側の温度勾配が急で、細かいメッシュが必要になる。
原子炉燃料棒の温度分布
原子炉UO₂燃料棒（直径8〜10 mm、λ≈3 W/m·K）の内部発熱は円筒熱伝導で解析され、中心温度は表面より最大1800℃高くなることがある。1979年のスリーマイル島事故では燃料棒中心が融点（2865℃）を超えた部分があり、円筒熱伝導モデルによる事故後の温度再現解析がNRC報告書（NUREG-0772）にまとめられた。
円筒座標系の熱伝導のソフトウェア比較
商用ツールでの実装
円筒の熱伝導解析は各ツールでどう設定しますか？
軸対称モデルの設定方法を比較する。
| ツール | 要素タイプ | 設定方法 |
|---|---|---|
| Ansys Mechanical | PLANE55 (KEYOPT(3)=1) | 2D軸対称として解析タイプを指定 |
| Abaqus | DCAX4, DCAX8 | Part定義時にAxisymmetricを選択 |
| COMSOL | 2D Axisymmetric | モデルウィザードで座標系を選択 |
| STAR-CCM+ | 2D Axisymmetric mesh | メッシュモデルで軸対称を指定 |
APDL実装例
蒸気配管の放熱量計算例を示す。
```
/PREP7
ET,1,PLANE55,,,1 ! 軸対称
! 配管壁 (SUS304)
MP,KXX,1,16.3
! 断熱材
MP,KXX,2,0.05
CYL4,0,0,48.6,0,57.15,90 ! 配管壁
CYL4,0,0,57.15,0,107.15,90 ! 断熱材
AGLUE,ALL
ESIZE,1
AMESH,ALL
/SOL
SFL,内面LINE,CONV,500,180 ! 蒸気側 h=500, T=180℃
SFL,外面LINE,CONV,10,25 ! 外面 h=10, T=25℃
SOLVE
```
内面の $h$ が500で外面が10。桁が全然違いますね。
蒸気の強制対流は $h = 100$〜$1000$ W/(m$^2$ K)、外面の自然対流は $h = 5$〜$15$ 程度。だから温度降下のほとんどは断熱材で生じる。
Abaqusでの設定
Abaqusでは*FILM条件で対流を定義する。
```
*STEP
*HEAT TRANSFER, STEADY STATE
*FILM
inner_surface, F2, 180., 500.
outer_surface, F2, 25., 10.
*END STEP
```
記法が異なるだけで物理的な設定は同じですね。
その通り。ソルバーの記法は異なるが、フーリエの法則に基づく離散化は共通だ。結果も同一になる。
ABAQUS軸対称要素CAX4の威力
Abaqus（Dassault Systèmes）の軸対称固体要素CAX4は、円筒形状の熱伝導解析を2Dメッシュで高精度に解ける。これを使うと完全3Dモデルと比べ節点数を99%削減できる場合がある。石油化学プラント向け厚肉配管（肉厚100 mm超）の熱応力解析でこのアプローチを採用した事例がAbaqusユーザーコンファレンス2019で発表されている。
円筒座標系の熱伝導の先端研究
異方性円筒の熱伝導
CFRP製の円筒だと方向によって $k$ が違いますよね。
円筒座標での異方性熱伝導方程式は
CFRPの巻き方（ヘリカル、フープ、0度）で $k_r$、$k_\phi$、$k_z$ が大きく異なる。
| 方向 | $k$ [W/(m K)] |
|---|---|
| 繊維方向 | 5〜10 |
| 直交方向 | 0.5〜1 |
| 厚さ方向 | 0.3〜0.7 |
10倍以上の異方性ですね。
CFRP圧力容器（水素タンク、ロケットモーターケース）の熱解析では異方性の考慮が必須だ。Ansysでは材料座標系を円筒座標に設定してKXX/KYY/KZZを個別定義する。
接触問題を含む多層円筒
焼きばめや圧入の円筒では、接触圧力が温度によって変化し、接触熱抵抗も変動する。構造-熱の連成解析が必要だ。
温度が上がると熱膨張で接触圧力が変わり、接触熱抵抗が変わり、温度分布が変わる...と連鎖するんですね。
まさにその通り。Ansys MechanicalやAbaqusではSequentially Coupled Thermal-Stress解析でこの連成を扱える。1回の解析で両方の物理を同時に解くFully Coupled解析も可能だ。
回転体の熱伝導
回転機械（タービンディスク、電動機ロータ）では回転による遠心力で応力が生じ、さらに摩擦熱や風損も発生する。回転座標系での熱伝導と流体の対流をCHT解析で連成させる。
ガスタービンのディスクは典型的な円筒座標の問題ですね。
Ansys CFXやSTAR-CCM+ではRotating Reference Frameを設定してCHT解析が可能だ。ディスクの温度分布がクリープ寿命に直結するため、精度の高い解析が要求される。
回転する円筒での熱伝導と遠心力
高速回転するガスタービンディスクでは、遠心力による弾性歪みと熱膨張が連成するため、温度分布を求める際に熱伝導方程式と弾性方程式を同時に解く必要がある。GE Aviationは1980年代にT700エンジンタービンディスクの軸対称熱-構造連成解析を構築し、この手法が今日のAnsys Coupled Field解析の原型になったとされる。
円筒座標系の熱伝導のトラブル対応
よくあるトラブルと対策
円筒の熱伝導解析で注意すべき点を教えてください。
頻出トラブルを整理しよう。
1. 2Dと3Dで結果が合わない
原因: 2Dモデルで軸対称設定を忘れている。平面応力/平面ひずみのデフォルトでは $r$ の効果が入らない。
対策: Ansys PLANE55ならKEYOPT(3)=1、Abaqusなら軸対称要素を選択。COMSOLならモデル作成時に「2D Axisymmetric」を選ぶ。
2. 多層構造で温度分布が不連続
各層の接続がうまくいかないことがありますか？
原因: 層間のノードが共有されていない、またはTied Contact/Bonded Contactの設定漏れ。
対策: 共通面でノードマージするか、接触条件を適切に設定。意図的に接触熱抵抗を入れる場合はGap Conductanceを設定する。
3. 内面対流係数の見積もりミス
配管内面の $h$ は流速、流体物性、管径で大きく変わる。
| 流体条件 | $h$ [W/(m$^2$ K)] |
|---|---|
| 空気・自然対流 | 5〜25 |
| 空気・強制対流 | 25〜250 |
| 水・強制対流 | 500〜10,000 |
| 蒸気凝縮 | 5,000〜50,000 |
| 沸騰 | 2,500〜75,000 |
水と空気で100倍以上違うんですね。
Dittus-Boelter式やGnielinski式で算出するのが基本だ。$\text{Nu} = 0.023 \text{Re}^{0.8} \text{Pr}^{0.4}$（Dittus-Boelter）。Re数を間違えると $h$ が桁で狂う。
4. 単位系の混乱
SI系（m）とmm系で $k$ の値が変わる。$k = 16.3$ W/(m K) はmm系では $k = 0.0163$ W/(mm K) だ。Ansys Mechanical のWorkbench環境はmm系がデフォルトなので要注意。
単位系は本当に要注意ですね。
理論解で検算する習慣をつければ、単位系のミスはすぐに発見できる。$q = 2\pi k L \Delta T / \ln(r_2/r_1)$ を手計算して結果と照合するだけで良い。
内径・外径の単位ミスが大事故に
配管の円筒熱伝導解析でriとroをmm入力すべきところをmで入力し、壁厚が1000倍に計算されてしまうヒューマンエラーは現場で繰り返される。1999年のNASA Mars Climate Orbiter墜落事故（ポンド力とニュートンの単位混在）を教訓に、CAEソルバーが入力単位を明示的に確認するダイアログを追加したメーカーが複数ある。
関連トピック
なった
詳しく
報告