2018 年 82 巻 6 号 p. 183-187
Cooling of aluminum plate, cylinder and sphere induced by heat transfer was simulated by use of direct solving method of differential equation. Representative sizes of the plate, cylinder and sphere as well as heat transfer coefficient were varied in the simulation. The results of simulation showed that the cooling rate of surface in Newton’s raw was the ratio of surface area times transfer-induced heat flux to volume times specific heat. Difference in temperature between center and surface was found to be a half Biot number multiplied by the surface temperature in common among plate, cylinder and sphere. Temperature distribution between center and surface was found to be always parabola, the top of which lie at center, causing smooth heat conduction within the volume. Difference between the present result of simulation and Heisler Chart was also made clear.
高温の物体の冷却はその表面から外部への熱伝達によって生じ，ニュートンの冷却法則が当てはまる．鋼の熱処理では冷却速度の重要性から，古くから種々研究されてきた1）．Russell2）は種々の平板厚さの1/2と熱伝達率との積について，中心の初期温度に対する相対温度と無次元熱拡散時間および中心から表面までの相対距離との関係を数表で示した．またHeisler3）は板厚中心の相対温度の対数と無次元熱拡散時間との関係をビオー数の逆数を変数とする直線群で示し，平板中の温度分布をビオー数の逆数と正弦曲線の大きさとの関係を図示した4）．これらの両者とも中心の冷却速度と板厚方向または丸棒の半径方向における温度分布とが変数分離法で計算されていて，その妥当性を検討する必要があると考えられる．そこでこの研究では固体の熱伝達冷却の一例として高温のアルミニウム平板，丸棒および球の熱伝達冷却を研究の対象とし，前進差分法5）を使って冷却のシミュレーションを行った．前進差分法では移動境界条件の熱伝導計算が可能であり，表面の熱伝達冷却との同時計算ができるから，前進差分法は変数分離法より優れた計算法であると考えた．シミュレーションでは平板の厚さ，丸棒および球の半径と熱伝達率とを変数とした．シミュレーション結果を解析し，ニュートンの冷却法則の具体的な機構を明らかにしようと試みた．この報告ではシミュレーションの方法と結果の解析について報告する．
アルミニウム平板，丸棒および球の熱伝達冷却のシミュレーションは市販のパソコンで前進差分法を用いて行った．Table 1 にシミュレーションで使った熱的定数の記号と数値6）を示す．
Thermal properties used.
前進差分法を使用するため平板厚さの1/2のx0を10等分し，両端を含んで中心から順に分割点に節点番号0～10を付け中心差分とする．熱拡散率をR＝λ/cρで表し，拡散微小時間をdtとする．Δx＝x0/10と表せば番号iの節点の温度変化dUiは一般に次式で表される5）．
| （1） |
| （2） |
| （3） |
丸棒の半径をr0とし，Δr＝r0/10とする．差分式は次式で表わされる5）．
| （4） |
| （5） |
| （6） |
球の半径をr0，Δr＝r0/10とする．球拡散の差分式は次式である5）．
| （7） |
| （8） |
| （9） |
Fig. 1 は初期温度U*＝600°C，外部温度U0＝0°C，x0＝2 cm，h＝0.1 W/cm2 °C（a）および 0.2 W/cm2 °C（b）の条件でシミュレーションを行い，冷却開始から2秒間隔でアナログ表示させた温度分布図である．この図から，時間の経過とともに温度分布の間隔が少しずつ減少し，同じ温度で比較すれば冷却速度および表面と中心との温度差はhにほぼ比例することがわかる．
Temperature distribution curves during cooling of aluminum plates with different heat transfer coefficient.
Fig. 2 は温度分布曲線に及ぼす板厚の半分x0の影響を示す．熱伝達率hが同じ 0.2 W/cm2°Cで，x0が2倍に増加すれば温度分布曲線群はhが半減した場合と同様であることを示している．
Temperature distribution curves during cooling of aluminum plates with different size.
Fig. 3 はFig. 1（a）および（b）の冷却開始から10秒後の（U0-Ui）1/2 vs. iの関係を示す．両者はいずれも直線関係になり，温度分布は板厚中心を頂点とする曲率負の放物線であることがわかる．放物線の1次微分が定数だから2次微分が0となる．したがってFickの法則によって冷却速度は中心からの距離に関係なく皆同じであり，Fig. 1 の温度分布曲線がほぼ平行移動していることと符合している．
Plots of square root of temperature difference between center and relative distance vs. relative distance from center.
いまFig. 1（a）の描画と同時に記録したデータから表面の分割要素に流入する熱流束と，そこから流出する熱流束とを計算してみる．たとえば冷却開始から10秒後の流入熱流束は 47.636 W/cm2，流出熱流束は 50.071 W/cm2で両者の比は0.9514となった．またこの比は最初の2秒間を除き冷却中ほとんど変化しなかった．しかしhおよびx0が増加すればこの比は少しずつ増加した．
いまこの比を0.95と仮定すれば，それは節点番号0～9の分割要素の合計体積率と同じであり，表面要素の体積率0.05と合わせて全体が同じ速度で冷却することとよく符合する．またこの流入熱流束と流出熱流束との差を表面要素の熱容量cρΔx/2で除した冷却速度は，表面要素の流出熱流束を全熱容量cρx0で除した冷却速度と等しいと考えられる．このような冷却では固体中の局部要素や表面要素に熱が蓄積または欠乏することがなく「定常冷却」であると考えられる．すなわち定常冷却では全体の冷却速度vは次式で表される．
| （10） |
| （11） |
Fig. 4 はFig. 1 と同じ条件のシミュレーションで記憶した節点番号0, 5, 10の要素の温度を時間順に並べた冷却曲線である．同図の（a）と（b）とを比較すれば熱伝達率の増加で冷却速度が増加することおよび節点番号5の要素の温度が中心の温度に接近していることが示されている．
Cooling curves of center, middle and surface of aluminum plates with different heat transfer coefficient.
Fig. 5 ではFig. 1（a），（b）と同じ条件でシミュレーション冷却したときの表面要素の4秒ごとの温度（○印）と式（11）で計算した冷却曲線とを比較した．冷却開始直後の放物線温度分布ができるまではシミュレーション温度のほうが低いが，両者の差はしだいになくなることがわかる．
Comparison of cooling curves of surface division by simulation with that by eq. （11） marked by circles.
定常冷却における平板中心の温度は式（3）を逆算して求められる．すなわち
| （12） |
| （13） |
Fig. 6（a），（b）は熱伝達率hが同じでx0が変化した場合を例として，式（13）で計算した温度分布（○印）とシミュレーションで得た温度分布（曲線）とを比較した．シミュレーションではx0が増加すれば定常冷却から遠ざかるが，それにも拘わらず，両者の温度分布はよく合致している．
Comparison of simulated curves of temperature distribution with that calculated by eq. （13） in plates with different size. h＝0.2 W/cm2°C
Fig. 7（a），（b）は半径r0＝2 cm，U*＝600°C，U0＝0°Cでhを変化させたときの2秒間隔で示した温度分布である．Fig. 1 の同じサイズ，同じh，同じ温度での温度分布と比較すれば，放物線の形（曲率）は同じであるが，放物線同士の間隔が（a），（b）ともに約2倍になっている．
Temperature distribution curves during cooling of aluminum cylinders with different heat transfer coefficient.
またFig. 7（a）の冷却条件で冷却開始後10秒における表面要素に流入，流出する熱流束比を計算してみると0.907となった．この熱流束比の値は節点番号0～9の体積率0.952＝0.9025に近似していて，シミュレーションが定常冷却に近似していることがわかる．
丸棒では体積/表面積の比が0.5r0であるから全体の冷却速度vは次式で表される．
| （14） |
またFig. 7 の温度分布曲線の間隔がFig. 1 の場合の約2倍になったことが式（14）から理解できる．
式（14）からニュートンの冷却法則は次式で表されることがわかる．
| （15） |
定常冷却時の温度分布は，表面要素に流入，流出する熱流束比が（内部/表面）の体積比に等しいことが平板と同様であるから次式で表される．
| （16） |
Fig. 1 とFig. 7 で示した同じ冷却条件での温度分布の形（曲率）が同じであることは，式（13）と式（16）との類似によっていて，冷却速度には無関係なことが理解できる．
3.3 球の計算結果と考察Fig. 8（a），（b）は球の半径r0＝2 cm，U*＝600°C，U0＝0°Cでhを変化させたときの2秒間隔で示した温度分布である．平板の温度分布を示すFig. 1 と比較すると放物線の曲率は同じであるが，間隔が約3倍になっている．またFig. 8（a）の冷却時間5秒における表面要素に流入，流出する熱流束比は0.842となった．この比の数値は0.953＝0.8574に相当すると考えると定常冷却の条件が満足される．
Temperature distribution curves during cooling of aluminum spheres with different heat transfer coefficient.
球では体積/表面積の比がr0/3であるから，定常冷却速度vは次式で表される．
| （17） |
このアルミニウム平板，丸棒および球の冷却シミュレーションで得た結果は他の物質の物体にも汎用できると考えられる．平板，丸棒および球の熱伝達による冷却速度vを表す式（10），（14）および（17）は，それぞれのサイズをn分割し，表面積をS，体積をVで表せば次式に統一して表される．
| （18） |
| （19） |
平板，丸棒および球のいずれでも定常冷却における固体中の温度分布は中心を頂点とする曲率負の放物線であり，中心からi番目の要素の温度は次式で表される．
| （20） |
つぎにこのシミュレーション結果とHeisler線図3）とを比較検討する．Heisler線図は中心の相対温度の対数を縦軸，無次元熱拡散時間のフーリエ数を横軸とし，ビオー数の逆数をパラメータとする直線群で示されている．したがって平板の場合の直線の傾きは次式で表される．
この式の末尾の分数cρx0/hの逆数は式（10），（14）および式（17）では表面要素温度の冷却速度指数であり，中心要素温度の冷却速度指数ではない．しかしHeisler線図の直線の傾きがビオー数の逆数に逆比例して示されていることはよく理解できる．またHeisler線図から平板，丸棒および球の冷却を表す式（19）は次式でも表されることがわかった．
| （21） |
またHeisler線図およびRussellの数表では固体中の温度分布が正弦曲線で表されているが，正弦曲線では放物線のような定常冷却は不可能だと考えられる．またシミュレーションの結果から固体中の温度分布は表面温度とBi/2との積から容易に計算できることがわかった．
前進差分法を用いてアルミニウム平板，丸棒および球の熱伝達冷却のシミュレーションを行った結果，ニュートンの冷却法則の詳細を明らかにできた．すなわち
1）ニュートンの冷却を表す指数関数の指数は平板ではビオー数とフーリエ数との積で表される．丸棒ではその指数は2倍，球では3倍になる．
2）平板，丸棒および球の熱伝達冷却中における中心と表面との温度差は表面温度とビオー数の1/2との積であり，温度分布は中心を頂点とする曲率負の放物線である．
またHeisler線図とこのシミュレーション結果との比較検討を行い，両者の相違を明らかにした．
この研究のシミュレーションには教育用ソフト「十進ベーシック」を使用した．制作者白井一雄氏に深く謝意を表したい．
