球殻の定常熱伝導

カテゴリ: 熱解析 > 定常熱伝導 | 更新 2026-04-12
Steady-state temperature distribution in a spherical shell showing 1/r profile from inner to outer radius
球殻の定常熱伝導 ― 内面から外面への温度分布は1/rに従う

理論と物理

球殻と円筒の違い

🧑‍🎓

先生、球殻の定常熱伝導って円筒座標の場合と何が違うんですか? どちらも「丸い形」なので似た話になりそうですけど。

🎓

一番大きな違いは伝熱面積の増え方だ。円筒では面積が半径に比例($A = 2\pi r L$)するけど、球では面積が半径の2乗に比例($A = 4\pi r^2$)する。

🧑‍🎓

面積の増え方が違うと、温度分布も変わるんですか?

🎓

その通り。円筒では温度が$\ln r$(対数関数)で変わるのに対して、球殻では$1/r$で変わる。直感的に言うと、球殻の方が外側へ行くほど面積が急激に増えるから、温度降下が急速に減衰するんだ。

形状面積 $A(r)$温度分布熱抵抗の形
平板一定$T \propto r$(線形)$L/(kA)$
円筒$\propto r$$T \propto \ln r$$\ln(r_2/r_1)/(2\pi k L)$
球殻$\propto r^2$$T \propto 1/r$$(1/r_1 - 1/r_2)/(4\pi k)$
🧑‍🎓

なるほど、面積の「べき数」が上がるほど温度降下が緩やかになるんですね。球殻って実務ではどんな場面で使うんですか?

🎓

代表的なのは3つだ。

  • LNG球形タンク(直径40m級)の断熱設計 ― ボイルオフガス量の予測
  • 原子炉圧力容器の肉厚方向温度分布 ― 熱応力評価の入力
  • 球形燃料ペレット(HTGR: 高温ガス炉のTRISO燃料)の中心温度推定

いずれも球対称を仮定できるため、1次元の解析解で設計の第一近似が出せるのが強みだよ。

支配方程式と解析解

🧑‍🎓

じゃあ、球殻の支配方程式を教えてください。円筒のときは$1/r$の項が付きましたよね。

🎓

球座標で定常・径方向のみの1次元熱伝導方程式はこうなる。

$$ \frac{1}{r^2}\frac{d}{dr}\left(k r^2 \frac{dT}{dr}\right) = 0 $$

円筒の$\frac{1}{r}\frac{d}{dr}\left(kr\frac{dT}{dr}\right)=0$と比べると、$r$の指数が1から2に変わっている。これが面積が$r^2$に比例することの反映だ。

🧑‍🎓

$k$が一定なら、どうやって積分するんですか?

🎓

$k$一定の場合、$\frac{d}{dr}\left(r^2 \frac{dT}{dr}\right) = 0$ を2回積分すると

$$ T(r) = -\frac{C_1}{r} + C_2 $$

境界条件 $T(r_1) = T_1$, $T(r_2) = T_2$ を代入すると、温度分布の解析解が得られる。

$$ \boxed{T(r) = T_1 + (T_2 - T_1) \frac{\dfrac{1}{r} - \dfrac{1}{r_1}}{\dfrac{1}{r_2} - \dfrac{1}{r_1}}} $$

円筒の$\ln r$分布とのアナロジーで、$\ln r \to 1/r$に置き換わった形だ。

🧑‍🎓

$1/r$分布って、具体的にはどんな形のグラフになるんですか?

🎓

例えば内径100mm・外径200mmの鋼球殻で $T_1 = 500$°C、$T_2 = 100$°C の場合、

  • $r = 100$ mm: $T = 500$°C(内面)
  • $r = 133$ mm(肉厚1/3地点): $T = 300$°C
  • $r = 200$ mm: $T = 100$°C(外面)

平板なら肉厚1/3地点は$(500-100)\times(1-1/3)+100 = 367$°Cだから、球殻の方が内面側の温度降下が急だ。内面付近の熱流束が高いということでもある。

球殻の熱抵抗

🧑‍🎓

球殻を通過する熱流量はどう計算するんですか? 円筒と同じように熱抵抗の概念が使えます?

🎓

使える。フーリエの法則 $q = -kA(r)\frac{dT}{dr}$ に$A(r) = 4\pi r^2$を代入して積分すると

$$ Q = \frac{4\pi k (T_1 - T_2)}{\dfrac{1}{r_1} - \dfrac{1}{r_2}} $$

$Q = \Delta T / R$ の形に整理すると、球殻の伝導熱抵抗は

$$ \boxed{R_{\text{sphere}} = \frac{1}{4\pi k}\left(\frac{1}{r_1} - \frac{1}{r_2}\right)} $$

円筒と違って長さ$L$が登場しないのが球殻の特徴だ。球は全方向に閉じた形状だから、管の長さに相当するパラメータが不要になる。

🧑‍🎓

内面・外面に対流がある場合は、対流の熱抵抗を足せばいいんですか?

🎓

そう。球面の対流熱抵抗は $R_{\text{conv}} = 1/(4\pi r^2 h)$ だから、内面の流体温度$T_{\infty,1}$から外面の流体温度$T_{\infty,2}$まで

$$ Q = \frac{T_{\infty,1} - T_{\infty,2}}{\dfrac{1}{4\pi r_1^2 h_1} + \dfrac{1}{4\pi k}\left(\dfrac{1}{r_1} - \dfrac{1}{r_2}\right) + \dfrac{1}{4\pi r_2^2 h_2}} $$

対流+伝導+対流の直列回路だ。電気回路のオームの法則と全く同じ構造だね。

多層球殻モデル

🧑‍🎓

LNGタンクみたいに断熱材が挟まっている場合は、多層で計算するんですよね?

🎓

$n$層の球殻で各層の熱伝導率が$k_i$、半径が$r_0, r_1, \ldots, r_n$のとき、全伝導熱抵抗は

$$ R_{\text{total}} = \sum_{i=1}^{n} \frac{1}{4\pi k_i}\left(\frac{1}{r_{i-1}} - \frac{1}{r_i}\right) $$

対流を含めた総括伝熱係数$U$は、基準面積を内面$A_1 = 4\pi r_0^2$に取ると

$$ \frac{1}{U A_1} = \frac{1}{h_1 A_1} + \sum_{i=1}^{n}\frac{1}{4\pi k_i}\left(\frac{1}{r_{i-1}} - \frac{1}{r_i}\right) + \frac{1}{h_2 A_n} $$
🧑‍🎓

これを使えば、断熱材の厚さを変えたときの熱損失をパラメトリックに計算できるわけですね。

🎓

そうだ。LNG球形タンク(直径40m)を例にとると、鋼殻30mm($k = 50$ W/(m·K))+パーライト断熱材500mm($k = 0.04$ W/(m·K))で、断熱材の熱抵抗が鋼殻の約2,500倍を占める。支配的な層がどこかを定量的に把握できるのが熱抵抗法の強みだ。

内部発熱がある場合

🧑‍🎓

原子炉の燃料球みたいに内部で発熱している場合はどうなりますか?

🎓

一様な体積発熱$\dot{q}$ [W/m³]がある中実球(半径$R$、表面温度$T_s$)の場合、支配方程式は

$$ \frac{1}{r^2}\frac{d}{dr}\left(kr^2\frac{dT}{dr}\right) + \dot{q} = 0 $$

中心の対称条件 $dT/dr|_{r=0} = 0$ と$T(R) = T_s$で解くと

$$ T(r) = T_s + \frac{\dot{q}}{6k}(R^2 - r^2) $$

中心温度は $T_{\max} = T_s + \dot{q} R^2 / (6k)$ だ。円筒の場合は分母が$4k$だったから、球の方が中心温度上昇が小さい。表面積が大きいから放熱しやすいんだ。

🧑‍🎓

高温ガス炉のTRISO燃料粒子(直径約1mm)なら、この式で中心温度を見積もれるんですね。

🎓

そう。UO₂カーネル($k \approx 3$ W/(m·K)、$\dot{q} \approx 2 \times 10^8$ W/m³、$R = 0.25$ mm)なら $\Delta T_{\max} = 2\times10^8 \times (0.25\times10^{-3})^2 / (6 \times 3) \approx 0.7$°C。粒子が小さいので温度上昇もわずかだ。ただし実際は多層被覆(SiC層等)の接触熱抵抗を加味する必要がある。

Coffee Break よもやま話

球殻熱伝導と地球科学

球殻の定常熱伝導の理論は19世紀のポワソン(1820年代)やケルビン卿に遡る。ケルビンは地球を冷却する球体とみなし、定常熱伝導の$1/r$分布から地球の年齢を推定しようとした(約1億年と見積もったが、放射性崩壊による内部発熱を知らなかったため大幅に過小評価した)。現代でも地球のマントル対流モデルの初期近似に球殻熱伝導が使われている。

球座標ラプラシアンの導出
  • ラプラシアン$\nabla^2 T$(球座標):直交座標の$\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}+\frac{\partial^2 T}{\partial z^2}$を球座標$(r,\theta,\phi)$に変換すると $$\nabla^2 T = \frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial T}{\partial r}\right) + \frac{1}{r^2\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial T}{\partial\theta}\right) + \frac{1}{r^2\sin^2\theta}\frac{\partial^2 T}{\partial\phi^2}$$ 球対称($T$が$r$のみの関数)なら$\theta$項と$\phi$項が消え、$\frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{dT}{dr}\right) = 0$が残る。
  • $r^2$の物理的意味:半径$r$の位置での球面の面積$4\pi r^2$に由来する。一定の熱流量$Q$が通過する場合、熱流束$q'' = Q/(4\pi r^2)$は$r^2$に反比例して減少する。
仮定条件と適用限界
  • 定常状態:時間変化なし($\partial T/\partial t = 0$)。蓄熱項が無視できる十分な時間が経過した後の温度場。
  • 球対称:温度が$r$のみの関数($\theta$, $\phi$に依存しない)。局所的な熱源や非均一な対流がある場合は適用不可。
  • 等方性・均質:各層内で熱伝導率$k$が一定。温度依存性が大きい場合(鋼の$k$は0〜800°Cで50→25 W/(m·K)に変化)は非線形解析が必要。
  • フーリエの法則:$q = -k \nabla T$。極低温(超流動ヘリウム)やフェムト秒レーザー加熱では非フーリエモデルが必要。
次元解析と主要パラメータ
変数SI単位典型値・注意点
温度 $T$K または °C放射を含む場合は必ず絶対温度[K]を使用
熱伝導率 $k$W/(m·K)鋼≈50、断熱材≈0.03〜0.05、UO₂≈3
内径 $r_1$m球殻は$r_1 > 0$(中実球は別の解)
球殻熱抵抗 $R$K/W$(1/r_1 - 1/r_2)/(4\pi k)$ ― 長さ$L$を含まない
体積発熱 $\dot{q}$W/m³核燃料≈$10^8$、ジュール発熱≈$10^5$〜$10^7$

数値解法と実装

軸対称FEMモデル

🧑‍🎓

球殻の熱伝導をFEMで解くとき、3Dの球体モデルを作らないとダメですか?

🎓

球対称なら2D軸対称モデルで十分だ。断面の半円を描いて、対称軸を中心に回転させたモデルに相当する。計算コストは3Dモデルの1/1000以下になる。

モデルタイプ要素数目安自由度計算時間
3D全体(球全体)50万〜50万〜分オーダー
3D 1/8対称6万〜6万〜秒〜分
2D軸対称200〜200〜瞬時
1D解析解なしなし手計算
🧑‍🎓

軸対称モデルでは何に注意すればいいですか?

🎓

注意点は3つだ。

  • 要素タイプの選択:Ansys MechanicalならPLANE55(線形)/ PLANE77(2次)のKEYOPT(3)=1(軸対称)。AbaqusならDCAX4 / DCAX8。
  • 対称軸上の境界条件:$r=0$軸上には断熱条件($\partial T/\partial r = 0$)が自動的に適用されるが、ソフトによっては明示的な設定が必要。
  • メッシュのバイアス:内面近傍は温度勾配が急なので細かく、外面側は粗くするのが効率的。内面側要素サイズを外面側の1/3〜1/5にするバイアスメッシュが推奨。

有限差分法での離散化

🧑‍🎓

FEMじゃなくて、有限差分法で自前実装したい場合はどうすればいいですか? Pythonで書いてみたいんですが。

🎓

球座標の定常熱伝導方程式を展開すると

$$ \frac{d^2T}{dr^2} + \frac{2}{r}\frac{dT}{dr} = 0 $$

これを中心差分で離散化すると、節点$i$(位置$r_i$)について

$$ \frac{T_{i+1} - 2T_i + T_{i-1}}{\Delta r^2} + \frac{2}{r_i}\frac{T_{i+1} - T_{i-1}}{2\Delta r} = 0 $$

整理すると三重対角行列の連立方程式になる。

$$ \left(1 - \frac{\Delta r}{r_i}\right)T_{i-1} - 2T_i + \left(1 + \frac{\Delta r}{r_i}\right)T_{i+1} = 0 $$

$r_i$が小さい(内面付近)ほど$\Delta r / r_i$の影響が大きくなるから、等間隔メッシュだと内面近傍の精度が落ちる。不等間隔メッシュか、変数変換$s = 1/r$を使うと精度が改善する。

🧑‍🎓

$s = 1/r$の変換って面白いですね。変換するとどうなるんですか?

🎓

$s = 1/r$と置くと、球殻の方程式が$\frac{d^2T}{ds^2} = 0$になる。つまり$s$空間では温度分布が線形になるんだ。等間隔の$s$グリッドで離散化すれば、中心差分が厳密解を与える。この変換は球殻の熱伝導に特有のテクニックだ。

メッシュ設計の指針

🧑‍🎓

FEMでメッシュを切るとき、何を基準に要素サイズを決めればいいですか?

🎓

球殻の定常熱伝導では以下を指針にするといい。

項目推奨値理由
肉厚方向の分割数10層以上$1/r$分布を正確に再現するため
内面側のバイアス比1:3〜1:5内面の温度勾配が急峻
要素タイプ2次要素推奨$1/r$分布は線形要素では近似が粗い
多層界面節点を一致させる界面で温度勾配が不連続になるため
アスペクト比5以下過度に扁平な要素は精度低下
🧑‍🎓

メッシュ収束性の確認はどうやるんですか?

🎓

3水準以上のメッシュ密度で計算して、注目点の温度が収束しているか確認する。球殻の場合は解析解があるから、それと比較するのが最も確実だ。例えば肉厚方向5, 10, 20分割で内面の熱流束を比較して、変化率が1%以下になるまで細分化するのが一般的なプラクティスだ。

検証用ベンチマーク

🧑‍🎓

自分のFEMモデルが正しいか検証するためのベンチマーク問題はありますか?

🎓

典型的なベンチマーク問題を紹介しよう。

パラメータ
内径 $r_1$0.1 m
外径 $r_2$0.2 m
熱伝導率 $k$50 W/(m·K)
内面温度 $T_1$500°C
外面温度 $T_2$100°C

解析解:$r = 0.15$ m での温度は

$$ T(0.15) = 500 + (100-500)\frac{1/0.15 - 1/0.1}{1/0.2 - 1/0.1} = 500 + (-400)\frac{6.667 - 10}{5 - 10} = 500 - 266.7 = 233.3 \text{°C} $$

熱流量は $Q = 4\pi \times 50 \times 400 / (1/0.1 - 1/0.2) = 4\pi \times 50 \times 400 / 5 = 50{,}265$ W だ。FEMの結果がこれと0.1%以内で一致すれば検証OKだ。

実践ガイド

LNGタンクの断熱設計

🧑‍🎓

LNG球形タンクの断熱設計で、球殻の熱伝導計算は実際にどう使われているんですか?

🎓

LNG球形タンク(Moss型、直径40m級)の設計では、日量ボイルオフガス(BOG)を貯蔵容量の0.05%以下に抑えることが要求される。これを満たす断熱材厚さを球殻の多層熱伝導計算で決定する。

  • 鋼殻:9%Ni鋼、肉厚30mm、$k = 15$ W/(m·K)
  • 断熱材:パーライト粉末、厚さ500mm、$k = 0.04$ W/(m·K)
  • 内面条件:LNG温度 $-162$°C、$h_1 = 500$ W/(m²·K)(沸騰)
  • 外面条件:大気温度 $+35$°C、$h_2 = 10$ W/(m²·K)(自然対流)

この条件で多層球殻の公式を適用すると、侵入熱量は約80 kWとなり、BOG率0.04%程度に収まる。断熱材を400mmに減らすと100 kWを超えてBOG率が要求値を超過する。

🧑‍🎓

手計算で80kWと出た後、FEMで確認する意味はあるんですか?

🎓

実際のタンクは完全な球対称ではない。支持構造(スカート)からの伝熱パス、溶接部の断熱材の隙間、上下の温度差による自然対流など、1Dモデルでは捉えられない3D効果がある。FEMはこれらの局所的な影響を評価するために使う。手計算は全体像を把握する「設計のたたき台」、FEMは「詳細検討」という役割分担だ。

原子炉圧力容器の温度評価

🧑‍🎓

原子炉の設計では球殻の熱伝導はどう使われますか?

🎓

PWR(加圧水型軽水炉)の圧力容器下鏡は半球殻形状だ。通常運転時の肉厚方向温度分布は球殻の定常解で評価できる。ここから熱応力を求めて、圧力による膜応力と合算して応力強度を評価する。

さらに重要なのは加圧熱衝撃(PTS: Pressurized Thermal Shock)の評価だ。緊急冷却水注入時に圧力容器内面が急冷されると、内面側に引張応力が生じて既存の欠陥が成長するリスクがある。この非定常問題の初期条件として定常温度分布が使われる。

境界条件の設定

🧑‍🎓

球殻の熱伝導解析で、境界条件の設定で注意すべきことは?

🎓
境界条件タイプ数学的表現物理的意味注意点
第1種(Dirichlet)$T(r_i) = T_0$表面温度指定最も単純だが、現実には対流を含むことが多い
第2種(Neumann)$-k\frac{dT}{dr}\bigg|_{r_i} = q_0$熱流束指定断熱条件は$q_0 = 0$の特殊ケース
第3種(Robin)$-k\frac{dT}{dr}\bigg|_{r_i} = h(T - T_\infty)$対流熱伝達$h$の値が結果に大きく影響。文献値の確認が重要
第4種(輻射)$-k\frac{dT}{dr}\bigg|_{r_i} = \epsilon\sigma(T^4 - T_\infty^4)$熱放射非線形。高温(数百°C以上)で支配的になる
🧑‍🎓

$h$(熱伝達係数)の値をどう決めるかが実務上の大きな課題ですよね。

🎓

まさにそこが経験が必要な部分だ。球体外面の自然対流なら Churchill & Chu の相関式($\text{Nu} = 2 + 0.589 \text{Ra}^{1/4}/[1+(0.469/\text{Pr})^{9/16}]^{4/9}$)が使えるが、結果は$h$の仮定に対して線形的に効いてくる。$h$を20%間違えると侵入熱量も20%ずれる。感度分析をして、$h$の不確かさが結果にどう影響するか確認するのが必須だ。

品質保証チェックリスト

🧑‍🎓

球殻の熱伝導解析のレビューで見るべきポイントをまとめてもらえますか?

🎓
  • エネルギー保存:内面からの侵入熱量 = 外面からの放出熱量 ± 内部発熱 ― 残差1%以内か確認
  • 解析解との比較:均一物性・球対称のサブケースで解析解と一致するか
  • メッシュ収束性:3水準以上で注目量(温度、熱流束)の変化率が1%以内
  • 物性値の妥当性:温度範囲に対して適切な$k$値を使用しているか(LNG温度域の断熱材$k$はRT値と異なる)
  • 境界条件の整合性:対流係数$h$の算出根拠が明確か。相関式の適用範囲内か
  • 単位系の一貫性:特にmm/m系の混在がないか($k$は通常W/(m·K)だがモデルがmmの場合は注意)
  • Coffee Break よもやま話

    Moss型LNGタンクの断熱設計

    Moss型LNGタンク(ノルウェーKvaerner社が1970年代に開発)は直径40mの球形アルミタンクをパーライト断熱材で覆った構造だ。断熱材厚さ500mmで侵入熱量を約80kWに抑えるが、それでも1日あたり約30トンのLNGが気化する(BOG率0.04%)。このBOGはボイラー燃料として消費されるため無駄にはならないが、BOG率が大きすぎるとタンク内圧が上昇して安全弁が作動する問題が生じる。

    球殻の臨界断熱半径

    円筒に断熱材を巻くと、ある半径までは表面積増大による対流促進効果が断熱効果を上回り、かえって放熱が増える「臨界断熱半径」$r_c = k/h$が存在する。球殻の場合は$r_c = 2k/h$となり、円筒の2倍だ。実際には断熱材の$k$が0.04 W/(m·K)で$h = 10$ W/(m²·K)程度なら$r_c = 8$mmと極めて小さく、工業的に問題になることはほぼない。しかし、細い熱電対の保護管に断熱テープを巻くような場面では注意が必要だ。

    ソフトウェア比較

    商用ツールでの実装

    🧑‍🎓

    球殻の定常熱伝導を解けるCAEソフトにはどんな選択肢がありますか?

    🎓

    主要な商用ソフトと球殻モデリングでの特徴を比較しよう。

    ソフトウェア軸対称要素多層対応放射連成スクリプト
    Ansys MechanicalPLANE55/77材料番号で層分けSURF151/152APDL / PyAnsys
    AbaqusDCAX4/DCAX8Section割り当て*RADIATIONPython scripting
    COMSOL2D軸対称GUIでドメイン分割Surface-to-AmbientMATLAB LiveLink
    Code_AsterAXIS_FOURIER材料割り当てECHANGE_PAROIPython / comm file

    APDL実装例

    🧑‍🎓

    Ansys APDLで球殻の熱伝導を解く最小限のスクリプトを見せてもらえますか?

    🎓

    球殻($r_1 = 0.1$ m、$r_2 = 0.2$ m、$k = 50$ W/(m·K))の定常温度場を求めるAPDLの例だ。

    /PREP7
    ET,1,PLANE55,,,1    ! 軸対称熱伝導要素
    MP,KXX,1,50         ! 熱伝導率 50 W/(m·K)
    
    ! 球殻断面(半円)を作成
    K,1, 0.1, 0         ! 内径点(下端)
    K,2, 0.2, 0         ! 外径点(下端)
    K,3, 0, 0.2         ! 外径点(上端)
    K,4, 0, 0.1         ! 内径点(上端)
    LARC,1,4,0,0.1      ! 内面円弧
    L,4,3               ! 上端直線
    LARC,3,2,0,0.2      ! 外面円弧
    L,2,1               ! 下端直線
    AL,1,2,3,4          ! 面を定義
    
    ESIZE,,20            ! 肉厚方向20分割
    AMESH,ALL
    
    /SOLU
    ANTYPE,STATIC
    ! 境界条件
    NSEL,S,LOC,X,0.1    ! 内面節点
    D,ALL,TEMP,500       ! 内面温度 500°C
    NSEL,S,LOC,X,0.2    ! 外面節点
    D,ALL,TEMP,100       ! 外面温度 100°C
    ALLSEL
    SOLVE

    $r = 0.15$ m の温度が233.3°Cと出れば検証OKだ。

    Abaqusでの設定

    🧑‍🎓

    Abaqusだとどう設定しますか?

    🎓

    AbaqusのINPファイルでの設定例だ。

    *HEADING
    Spherical shell steady-state conduction
    
    *NODE
    ! 軸対称モデルの節点座標 (r, z)
    1, 0.1, 0.0
    ...
    
    *ELEMENT, TYPE=DCAX4, ELSET=SHELL
    ! 軸対称4節点拡散要素
    1, 1, 2, 12, 11
    ...
    
    *MATERIAL, NAME=STEEL
    *CONDUCTIVITY
    50.0
    
    *SOLID SECTION, ELSET=SHELL, MATERIAL=STEEL
    
    *STEP, NAME=STEADY
    *HEAT TRANSFER, STEADY STATE
    
    *BOUNDARY
    INNER_NODES, 11, 11, 500.0
    OUTER_NODES, 11, 11, 100.0
    
    *OUTPUT, FIELD
    *NODE OUTPUT
    NT
    *ELEMENT OUTPUT
    HFL
    
    *END STEP

    ポイントは要素タイプがDCAX4(軸対称4節点拡散要素)で、自由度番号11が温度であること。構造要素のDOF 1〜6とは異なるから注意しよう。

    オープンソースの選択肢

    🧑‍🎓

    商用ソフトのライセンスがない場合、フリーで使えるツールはありますか?

    🎓
    ソフトライセンス球殻対応備考
    Code_AsterGPL軸対称要素ありEDF開発。フランス語ドキュメントが多いが英語版も充実
    FEniCS / FEniCSxLGPLPython API で柔軟に定義可変分形式を直接記述。研究用途に最適
    ElmerGPL2D/3D熱伝導対応CSC(フィンランド)開発。GUIあり
    自前Python実装有限差分/SciPy教育目的なら最もシンプル。解析解で即座に検証可能

    ツール選定の3つの問い

    • 「純粋な球対称か?」:完全な球対称なら解析解で十分。FEMは不要。支持構造や局所的な熱源がある場合に初めてFEMの出番。
    • 「連成解析が必要か?」:熱応力評価が必要なら、熱→構造の片方向連成ができるソルバーを選ぶ。Ansys / Abaqusはこのワークフローが確立されている。
    • 「結果の信頼性をどう担保するか?」:規制対応(原子力、圧力容器)では、V&V(検証と妥当性確認)の文書化が必須。ASME V&V 10やNAFEMS ベンチマーク問題との比較記録を残すこと。

    トラブルシューティング

    温度分布が理論解と合わない

    🧑‍🎓

    先生、球殻のFEMモデルを作ったんですけど、解析解と5%くらいずれるんです。何が悪いんでしょう…

    🎓

    球殻で5%ずれるのは大きい。典型的な原因を順にチェックしよう。

    • 要素タイプの設定ミス:軸対称(KEYOPT(3)=1)を忘れて平面応力/ひずみモデルになっている。これが最も多い原因。$r^2$の面積効果が反映されないので、結果が全く違ってしまう。
    • 対称軸の配置:Ansys APDLでは$r$方向が$X$軸に対応するため、球殻断面をYZ平面に描いてしまうとモデルが回転しない。Abaqusも同様に$r=X$, $z=Y$。
    • 線形要素の精度不足:$1/r$分布を4節点の線形要素で近似すると、特に薄肉球殻($r_2/r_1 > 3$)で誤差が大きくなる。8節点の2次要素に切り替えるか、分割数を増やす。
    🧑‍🎓

    なるほど…軸対称の設定は確認します。他にありがちなミスは?

    🎓
    • $r = 0$付近の節点:軸対称モデルで$r = 0$上に要素があると、ヤコビアンが零になり精度が崩壊する。中実球を解く場合は$r = 0$から微小距離($10^{-6}$ m程度)オフセットする。
    • 対称面の境界条件:半球モデル($z \geq 0$のみ)の場合、$z = 0$面に断熱条件を忘れると結果がおかしくなる。ただし多くのソルバーは「何も指定しない面 = 断熱」なので、明示的に不要な場合もある。
    • 熱流束の不連続

      🧑‍🎓

      球殻の半径方向に沿って熱流束を出力したら、要素境界でガタガタしているんですが…

      🎓

      温度場は連続(C⁰連続)だが、温度勾配(=熱流束)は要素境界で不連続になるのがFEMの特徴だ。対策は3つ。

      • 節点平均:隣接要素の値を平均する。ほとんどのソフトがデフォルトで行う。
      • 2次要素の使用:勾配の近似精度が上がり、不連続幅が減少する。
      • メッシュ細分化:十分細かければ不連続幅は無視できるレベルになる。

      ポスト処理で「Unaveraged」(非平均)の熱流束を出力して、隣接要素間の差が5%以上ある箇所はメッシュ不足の兆候だ。

      多層モデルの接触熱抵抗

      🧑‍🎓

      多層球殻モデルで、層間の接触熱抵抗はどう扱うんですか? 手計算では完全接触を仮定しましたけど、実際は隙間がありますよね。

      🎓

      良い質問だ。層間の接触熱抵抗$R_c$ [m²·K/W]を考慮する方法は3つある。

      • 等価薄層:接触面に$k_c = \delta / R_c$の薄い仮想層を挿入する。$\delta$は層厚(例:0.1mm)。手計算と相性がよい。
      • 接触要素:Ansysの CONTA/TARGE 要素、AbaqusのTIE + GAP CONDUCTANCE で接触熱コンダクタンスを定義する。
      • 接触コンダクタンス:COMSOLでは「薄い層」(Thin Layer)機能で接触熱コンダクタンス$h_c = 1/R_c$ [W/(m²·K)]を直接入力する。

      接触面の仕上げ(研磨、グリス塗布等)により$R_c$は$10^{-5}$〜$10^{-2}$ m²·K/Wの範囲で変動する。面圧が上がるとマイクロ凸部が潰れて$R_c$が減少する。

      よくある単位ミス

      🧑‍🎓

      先生、先輩が「単位を間違えて3日無駄にした」って言ってたんですけど、球殻の計算で特にありがちな単位ミスは?

      🎓
      症状原因対策
      熱流量が1000倍大きい半径をmmで入力したのに$k$がW/(m·K)全てSI(m系)で統一。mm系なら$k$を mW/(mm·K) に換算
      温度分布が平板的(線形)軸対称モデルのKEYOPTを設定し忘れ解析解との比較で即座に気づける
      多層モデルの界面温度が不連続節点が共有されていない(メッシュが分離)Merge nodes / Tie constraint で接続
      外面熱流束がゼロ対流境界条件のFilm Coefficient $h$の単位間違い$h$が W/(m²·K) か確認。mm系なら 10⁻³ を乗じた値になる
      🧑‍🎓

      やっぱり解析解と比較するのが最初のチェックですね。球殻は解析解があるのが助かります。

      🎓

      その通り。解析解のある問題は「自分のモデルが正しいか」を検証できる貴重な機会だ。球殻の$T(r) = 1/r$分布で検証してから、複雑な3Dモデルに進むのがベストプラクティスだよ。

      デバッグの鉄則:「まず最小モデルで検証」

      球殻の熱伝導解析でトラブルが発生したら、まず複雑な3Dモデルではなく、2D軸対称の単層球殻モデル(解析解がある条件)を作って検証する。このモデルで解析解と一致すれば、要素タイプ・境界条件・単位系は正しい。そこから段階的に複雑化(多層化→3D化→連成化)していくことで、問題の原因を絞り込める。いきなり複雑なモデルでデバッグするのは「100変数の連立方程式を暗算で解く」ようなものだ。

      関連シミュレーター

      この分野のインタラクティブシミュレーターで理論を体感しよう

      シミュレーター一覧

      関連する分野

      構造解析流体解析製造プロセス解析
      この記事の評価
      ご回答ありがとうございます!
      参考に
      なった
      もっと
      詳しく
      誤りを
      報告
      参考になった
      0
      もっと詳しく
      0
      誤りを報告
      0
      Written by NovaSolver Contributors
      Anonymous Engineers & AI — サイトマップ
      プロフィールを見る