熱処理炉の温度むらと冷却条件を偏微分方程式で読み解く：製造業の数値計算10本ノック（No.081〜No.090）

本記事では、架空の金属部品工場における鋼板の加熱・冷却と炉内流れを題材に、偏微分方程式（PDE）を製造条件の判断へ結び付けます。有限差分法・有限要素法・有限体積法、熱伝導・波動・ポアソン・ナビエ–ストークス方程式、境界条件、CFL条件を、数式、Python、表、グラフで確認します。

狙いは高精度CAEの代替ではなく、温度むら、保持時間、センサー配置、計算の安定性について、モデルの仮定と判断材料を理解することです。

[!NOTE] 本資料は、数理工房 (もしくは代表である和山個人) が過去に企業研修において使用した notebook を企業様の許可を得て再構成・編集のうえ公開しています。 掲載データはすべて架空のものであり、実在する企業・工場・数値とは一切関係ありません。

はじめに：この記事で扱う製造業の実務課題

熱処理では、炉温の代表値が目標に達していても、鋼板の中心と端部、表面と内部に温度差が残ることがあります。温度むらは硬さ、残留応力、寸法精度に影響します。一方、加熱時間を安全側へ延ばせば、電力、処理能力、酸化の負担が増えます。

本記事では「どの位置がいつ規格温度へ入るか」「境界の熱の入り方をどう置くか」「時間刻みが計算結果を壊していないか」を中心に扱います。

現場でよくある状況

• 炉内熱電対は目標温度だが、製品中心温度を直接測れない
• 材質・板厚・装入量が変わるたびに、保持時間を経験で補正している
• 端部の過熱と中心部の昇温不足を同時に避けたい
• CAE結果はあるが、境界条件やメッシュによる差を説明しにくい

なぜこの問題は判断が難しいのか

温度や流速は時間だけでなく空間でも変わります。測定点の値だけでは点と点の間を保証できず、熱伝導率などの物性、対流熱伝達、固定温度・熱流束といった境界条件にも結果が依存します。さらに、離散化した計算には打切り誤差と安定性条件があります。したがって、実測、物理モデル、数値計算の三者を照合する必要があります。

今回扱うノックの全体像

No.	テーマ	製造業での判断例
081	偏微分方程式とは	時間・位置を含む品質場の定義
082	有限差分法	格子上の温度勾配・曲率の計算
083	有限要素法	形状や材料特性を行列へ組み込む
084	有限体積法	熱量収支をセル単位で守る
085	熱伝導方程式	保持時間と温度むらの評価
086	波動方程式	衝撃・振動の伝播時間の把握
087	ポアソン方程式	内部発熱を伴う定常温度分布
088	ナビエ・ストークス方程式	炉内循環と滞留域の把握
089	境界条件	炉壁・対流・断熱の仮定比較
090	CFL条件	安定な時間刻みの設計

Python 環境の準備

NumPy で配列計算、pandas で判断指標の表、SciPy で疎行列、matplotlib で可視化を行います。外部データは使わず、乱数シードを固定します。グラフのラベルは環境差による文字化けを避けるため英語表記とします。

%matplotlib inline
import platform
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from scipy import sparse
from scipy.sparse.linalg import spsolve
from IPython.display import display

rng = np.random.default_rng(81)
plt.rcParams["figure.figsize"] = (7.2, 4.2)
plt.rcParams["axes.grid"] = True
print({"python": platform.python_version(), "numpy": np.__version__,
       "pandas": pd.__version__, "matplotlib": matplotlib.__version__})

{'python': '3.13.1', 'numpy': '2.5.1', 'pandas': '3.0.3', 'matplotlib': '3.11.0'}

架空データの作成

長さ1.2 mの鋼板を一次元化し、21点の温度を考えます。基準物性は熱拡散率 lpha=1.2 imes10^{-5}\ \mathrm{m^2/s}、初期温度25 ℃、炉内温度850 ℃です。実務では温度依存物性、板厚方向、治具接触、放射を追加しますが、ここでは手法の比較を明確にするため定数とします。

L, nx = 1.2, 21
x = np.linspace(0, L, nx)
dx = x[1] - x[0]
alpha, T0, T_furnace = 1.2e-5, 25.0, 850.0
sensor_df = pd.DataFrame({
    "position_m": x,
    "initial_temp_C": T0 + rng.normal(0, 0.35, nx),
    "zone": np.where((x < 0.18) | (x > 1.02), "edge", "center")
})
display(sensor_df.iloc[[0, 1, 9, 10, 11, 19, 20]].round(2))

	position_m	initial_temp_C	zone
0	0.00	25.20	edge
1	0.06	25.27	edge
9	0.54	25.23	center
10	0.60	25.21	center
11	0.66	25.48	center
19	1.14	25.47	edge
20	1.20	25.07	edge

No.081：偏微分方程式とは

実務での意味

PDEは温度 $T(x,t)$ のように、複数の独立変数に依存する「場」を記述します。炉温の一点監視から、製品全体の温度履歴を考えるための共通言語です。

分析・モデル化の考え方

一次元熱伝導方程式は $\frac{\partial T}{\partial t}=\alpha\frac{\partial^2T}{\partial x^2}$ です。左辺は温度の時間変化、右辺は空間的な温度の曲がりです。初期条件と境界条件がそろって初めて解が定まります。ここでは解析解の一例 $T=T_0+A e^{-\alpha k^2t}\sin(kx)$ を使い、空間と時間の両方に依存することを確認します。

Pythonで確認する

times = np.array([0, 900, 3600, 10800])
k = np.pi / L
field = np.array([T0 + 500*np.exp(-alpha*k**2*t)*np.sin(k*x) for t in times])
fig, ax = plt.subplots()
for row, t in zip(field, times): ax.plot(x, row, marker="o", ms=3, label=f"t={t/60:.0f} min")
ax.set(title="Temperature field changes in space and time", xlabel="Position x [m]", ylabel="Temperature [degC]")
ax.grid(True); ax.legend(); fig.tight_layout(); plt.show()
display(pd.DataFrame({"time_min": times/60, "center_temp_C": field[:, nx//2], "range_C": np.ptp(field, axis=1)}).round(2))

	time_min	center_temp_C	range_C
0	0.0	525.00	500.00
1	15.0	489.33	464.33
2	60.0	396.86	371.86
3	180.0	230.68	205.68

結果の読み取り

同じ時刻でも位置によって温度が異なり、その差は時間とともに減衰します。「炉温が何℃か」だけでなく、「製品のどの位置が何分後に何℃か」を品質条件に結び付ける必要があります。この解析解は境界を簡略化した確認用であり、実炉の予測には物性と熱の出入りを同定します。

No.082：有限差分法

実務での意味

有限差分法（FDM）は、規則格子上で温度の傾きや曲率を近傍点の差に置き換えます。単純形状の熱履歴を素早く試算し、センサー間隔や計算格子の影響を調べるのに向きます。

分析・モデル化の考え方

中心差分では $\frac{\partial^2T}{\partial x^2}\bigg|_i\approx\frac{T_{i+1}-2T_i+T_{i-1}}{\Delta x^2}.$ 格子を細かくすると通常は離散化誤差が減りますが、計算量が増え、陽解法の許容時間刻みも小さくなります。

Pythonで確認する

def second_difference(y, h):
    return (y[2:] - 2*y[1:-1] + y[:-2]) / h**2

rows = []
for n in [11, 21, 41, 81]:
    xg = np.linspace(0, L, n); h = xg[1]-xg[0]
    y = np.sin(np.pi*xg/L)
    exact = -(np.pi/L)**2 * np.sin(np.pi*xg[1:-1]/L)
    approx = second_difference(y, h)
    rows.append({"grid_points": n, "dx_m": h, "max_abs_error": np.max(np.abs(approx-exact))})
fdm_error = pd.DataFrame(rows)
display(fdm_error.round(7))
fig, ax = plt.subplots(); ax.loglog(fdm_error.dx_m, fdm_error.max_abs_error, "o-")
ax.set(title="FDM grid convergence", xlabel="Grid spacing dx [m]", ylabel="Maximum absolute error")
ax.grid(True, which="both"); fig.tight_layout(); plt.show()

	grid_points	dx_m	max_abs_error
0	11	0.120	0.056186
1	21	0.060	0.014081
2	41	0.030	0.003523
3	81	0.015	0.000881

結果の読み取り

格子幅を半分にすると誤差がおおむね4分の1となり、中心差分の二次精度が確認できます。実務では最細格子を正解と決めつけず、主要KPI（最高温度、中心温度、規格到達時刻）が格子細分化で十分変わらなくなることを確認します。

No.083：有限要素法

実務での意味

有限要素法（FEM）は形状を要素へ分割し、各要素の寄与を全体行列へ組み立てます。複雑な製品形状、穴、異材接合、局所メッシュ細分化を扱うCAEの基本です。

分析・モデル化の考え方

定常一次元熱伝導 $-\frac{d}{dx}(k\frac{dT}{dx})=0$ の線形要素では、要素剛性行列は $K_e=\frac{k}{h_e}\begin{bmatrix}1&-1\\-1&1\end{bmatrix}.$ 要素ごとの熱伝導率を変えれば、異材の熱抵抗を表現できます。左右温度を固定し、母材と低熱伝導インサートの温度分布を比較します。

Pythonで確認する

nodes = np.array([0.0, 0.18, 0.42, 0.58, 0.76, 1.0, 1.2])
k_elem = np.array([45, 45, 8, 8, 45, 45], dtype=float)
K = np.zeros((len(nodes), len(nodes)))
for e, ke in enumerate(k_elem):
    h = nodes[e+1]-nodes[e]
    K[e:e+2, e:e+2] += ke/h*np.array([[1, -1], [-1, 1]])
T_fem = np.zeros(len(nodes)); T_fem[[0, -1]] = [850, 200]
free = np.arange(1, len(nodes)-1)
T_fem[free] = np.linalg.solve(K[np.ix_(free, free)], -K[np.ix_(free, [0, len(nodes)-1])] @ T_fem[[0, -1]])
fem_df = pd.DataFrame({"node_m": nodes, "temperature_C": T_fem})
display(fem_df.round(1))
fig, ax = plt.subplots(); ax.plot(nodes, T_fem, "o-"); ax.axvspan(0.42, 0.76, alpha=.15, label="low-k insert")
ax.set(title="FEM temperature across a composite plate", xlabel="Position x [m]", ylabel="Temperature [degC]")
ax.grid(True); ax.legend(); fig.tight_layout(); plt.show()

	node_m	temperature_C
0	0.0	850.0
1	0.2	807.8
2	0.4	751.5
3	0.6	540.5
4	0.8	303.2
5	1.0	246.9
6	1.2	200.0

結果の読み取り

低熱伝導部で温度勾配が急になります。平均温度だけでは見えない局所的な熱抵抗を、材料配置と対応付けて確認できます。本番FEMでは要素品質、接触熱抵抗、二次元・三次元形状、物性の温度依存性を検証し、メッシュ収束を報告します。

No.084：有限体積法

実務での意味

有限体積法（FVM）は各セルへの流入・流出を収支として計算します。熱、質量、運動量の保存を明示しやすく、炉内流れや冷却流路の計算で広く使われます。

分析・モデル化の考え方

セル $i$ の熱量収支は $\rho c_p V_i\frac{T_i^{n+1}-T_i^n}{\Delta t}=q_{i-1/2}-q_{i+1/2},\quad q=-kA\frac{\Delta T}{\Delta x}.$ 内部セル間の熱流束は片方から出た量が隣へ入るため、全熱量が保存されます。断熱端で一つの高温セルが拡散する例を確認します。

Pythonで確認する

ncv, dt, steps = 30, 2.0, 600
dxc = L/ncv; r = alpha*dt/dxc**2
T = np.full(ncv, 25.0); T[ncv//2] = 425.0
energy_history = [T.sum()]
snapshots = {0: T.copy()}
for n in range(1, steps+1):
    face_flux = -alpha*(T[1:]-T[:-1])/dxc
    Tn = T.copy()
    Tn[0] += dt*(-face_flux[0])/dxc
    Tn[-1] += dt*(face_flux[-1])/dxc
    Tn[1:-1] += dt*(face_flux[:-1]-face_flux[1:])/dxc
    T = Tn; energy_history.append(T.sum())
    if n in [50, 200, 600]: snapshots[n] = T.copy()
xc = (np.arange(ncv)+.5)*dxc
fig, ax = plt.subplots()
for n, temp in snapshots.items(): ax.plot(xc, temp, label=f"t={n*dt:.0f}s")
ax.set(title="FVM diffusion with insulated boundaries", xlabel="Cell center x [m]", ylabel="Temperature [degC]")
ax.grid(True); ax.legend(); fig.tight_layout(); plt.show()
print(f"Relative heat-balance error: {(energy_history[-1]/energy_history[0]-1):.3e}")

Relative heat-balance error: -4.441e-16

結果の読み取り

ピーク温度は下がりながら周囲へ広がり、断熱境界なので全セル温度和に対応する熱量は数値誤差の範囲で保存されます。CFDや熱収支では、色の見た目だけでなく入口・出口・壁・蓄積の収支誤差を受入基準にします。

No.085：熱伝導方程式

実務での意味

熱伝導方程式により、板中心が規格下限へ到達する時刻と、その時点の温度むらを予測できます。これは保持時間、ライン速度、エネルギーのトレードオフを定量化する基礎です。

分析・モデル化の考え方

両端を炉温に固定し、陽的差分 $T_i^{n+1}=T_i^n+r(T_{i+1}^n-2T_i^n+T_{i-1}^n),\quad r=\alpha\Delta t/\Delta x^2$ を使います。一次元熱拡散では $r\le 1/2$ が安定性の目安です。中心温度が780 ℃へ達するまで計算します。

Pythonで確認する

dt_heat = 120.0
r_heat = alpha*dt_heat/dx**2
T = np.full(nx, T0); T[[0, -1]] = T_furnace
history = [(0.0, T.copy())]; reach_s = None
for n in range(1, 3001):
    T[1:-1] += r_heat*(T[2:]-2*T[1:-1]+T[:-2])
    if n in [5, 20, 60, 120, 240]: history.append((n*dt_heat, T.copy()))
    if reach_s is None and T[nx//2] >= 780: reach_s = n*dt_heat; history.append((n*dt_heat, T.copy())); break
fig, ax = plt.subplots()
for ts, temp in history: ax.plot(x, temp, label=f"{ts/60:.0f} min")
ax.axhline(780, color="black", ls="--", lw=1, label="lower spec")
ax.set(title="Transient heating of the plate", xlabel="Position x [m]", ylabel="Temperature [degC]")
ax.grid(True); ax.legend(ncol=2); fig.tight_layout(); plt.show()
last = history[-1][1]
display(pd.DataFrame([{"r": r_heat, "center_reach_min": reach_s/60, "min_C": last.min(), "max_C": last.max(), "range_C": np.ptp(last)}]).round(2))

	r	center_reach_min	min_C	max_C	range_C
0	0.4	548.0	780.4	850.0	69.6

結果の読み取り

端部は境界条件により直ちに炉温となる一方、中心の規格到達には長い時間が必要です。中心到達時刻だけで運転条件を決めると端部過熱を見落とすため、最低温度、最高温度、規格内滞在時間を併記します。実機では放射・対流境界と実測中心温度で校正してください。

No.086：波動方程式

実務での意味

波動方程式は、プレス衝撃、棒材の打撃検査、設備振動がどの速度で伝わり反射するかを表します。センサー到達時間から異常位置を推定する際の基礎になります。

分析・モデル化の考え方

一次元波動方程式 $\frac{\partial^2u}{\partial t^2}=c^2\frac{\partial^2u}{\partial x^2}$ を中央差分化します。波速 $c$ と格子条件には $c\Delta t/\Delta x\le1$ が必要です。中央の変位パルスが左右へ伝播する様子を計算します。

Pythonで確認する

nw, c = 121, 5000.0
xw = np.linspace(0, L, nw); dxw = xw[1]-xw[0]; dtw = 0.85*dxw/c
u0 = np.exp(-((xw-0.35)/0.035)**2); u0[[0,-1]] = 0
u_prev = u0.copy(); u = u0.copy(); wave_snaps = {0: u.copy()}
targets = [30, 60, 90, 120]
for n in range(1, 121):
    un = np.zeros_like(u)
    un[1:-1] = 2*u[1:-1]-u_prev[1:-1]+(c*dtw/dxw)**2*(u[2:]-2*u[1:-1]+u[:-2])
    u_prev, u = u, un
    if n in targets: wave_snaps[n] = u.copy()
fig, ax = plt.subplots()
for n, z in wave_snaps.items(): ax.plot(xw, z, label=f"{n*dtw*1e6:.0f} us")
ax.set(title="Impact wave propagation and reflection", xlabel="Position x [m]", ylabel="Normalized displacement")
ax.grid(True); ax.legend(); fig.tight_layout(); plt.show()
print(f"Theoretical travel time over 0.60 m: {0.60/c*1e6:.1f} microseconds; Courant number: {c*dtw/dxw:.2f}")

Theoretical travel time over 0.60 m: 120.0 microseconds; Courant number: 0.85

結果の読み取り

初期パルスは左右へ分かれ、固定端で反射します。到達時間は距離÷波速に対応するため、複数センサーの時刻差から発生位置を絞れます。ただし実材では分散、減衰、断面変化、センサー応答があるため、既知位置の打撃試験で波速と検出閾値を校正します。

No.087：ポアソン方程式

実務での意味

内部発熱を伴う定常温度、静電場、圧力補正などはポアソン方程式へ帰着します。ヒーターや反応熱がある部品で、定常時のホットスポットを見積もれます。

分析・モデル化の考え方

二次元定常熱伝導を $\nabla^2T=-q'''/k$ とし、外周温度を固定します。5点差分で疎な連立一次方程式を作り、内部発熱が一様な場合と右上に偏る場合を比較します。

Pythonで確認する

ny2, nx2 = 25, 35
xx = np.linspace(0, 1.4, nx2); yy = np.linspace(0, 1.0, ny2)
hx, hy = xx[1]-xx[0], yy[1]-yy[0]
X, Y = np.meshgrid(xx, yy)
source = 1.0 + 3.0*np.exp(-((X-1.05)**2+(Y-.72)**2)/.035)
N = (nx2-2)*(ny2-2)
A = sparse.lil_matrix((N, N)); b = np.zeros(N)
def idx(j, i): return (j-1)*(nx2-2)+(i-1)
for j in range(1, ny2-1):
    for i in range(1, nx2-1):
        p=idx(j,i); A[p,p] = -2/hx**2-2/hy**2; b[p] = -source[j,i]
        for jj,ii,w in [(j,i-1,1/hx**2),(j,i+1,1/hx**2),(j-1,i,1/hy**2),(j+1,i,1/hy**2)]:
            if 1 <= ii < nx2-1 and 1 <= jj < ny2-1: A[p,idx(jj,ii)] = w
Tpoi = np.zeros((ny2,nx2)); Tpoi[1:-1,1:-1] = spsolve(A.tocsr(), b).reshape(ny2-2, nx2-2)
hot = np.unravel_index(np.argmax(Tpoi), Tpoi.shape)
fig, ax = plt.subplots(); cs=ax.contourf(X,Y,Tpoi,levels=18,cmap="inferno"); fig.colorbar(cs,ax=ax,label="Temperature rise [a.u.]")
ax.plot(xx[hot[1]],yy[hot[0]],"co",label="hot spot"); ax.set(title="Poisson solution with nonuniform heat source",xlabel="x [m]",ylabel="y [m]")
ax.grid(True); ax.legend(); fig.tight_layout(); plt.show()
print({"hotspot_x_m": round(xx[hot[1]],3), "hotspot_y_m": round(yy[hot[0]],3), "max_rise_au": round(Tpoi[hot],3)})

{'hotspot_x_m': np.float64(0.947), 'hotspot_y_m': np.float64(0.625), 'max_rise_au': np.float64(0.136)}

結果の読み取り

最高温度位置は幾何中心ではなく、発熱が偏った側へ移ります。センサーを中央だけに置くと最大値を見逃す可能性があります。発熱分布と境界温度の不確かさを振り、ホットスポット位置がどの範囲で動くかを確認して配置を決めます。

No.088：ナビエ・ストークス方程式

実務での意味

炉内ファン、冷却ノズル、洗浄槽では、流れの偏りが熱・物質移動のむらにつながります。ナビエ–ストークス方程式は速度と圧力の場を記述し、滞留域や循環を評価する基礎です。

分析・モデル化の考え方

非圧縮性流体では $\frac{\partial\mathbf{u}}{\partial t}+(\mathbf{u}\cdot\nabla)\mathbf{u}=-\frac{1}{\rho}\nabla p+\nu\nabla^2\mathbf{u},\qquad \nabla\cdot\mathbf{u}=0.$ ここでは教育用に、上壁が動く正方形キャビティを圧力ポアソン法で解きます。実務CFDより粗い格子・短い反復ですが、循環と低速域の意味を確認できます。

Pythonで確認する

n=31; nt=300; nit=40; rho=1.; nu=.1; dt=.001
dxn=2/(n-1); dyn=dxn
u=np.zeros((n,n)); v=np.zeros_like(u); p=np.zeros_like(u)
for _ in range(nt):
    un=u.copy(); vn=v.copy()
    b=np.zeros_like(p)
    b[1:-1,1:-1]=rho*(1/dt*((un[1:-1,2:]-un[1:-1,:-2])/(2*dxn)+(vn[2:,1:-1]-vn[:-2,1:-1])/(2*dyn))
      -((un[1:-1,2:]-un[1:-1,:-2])/(2*dxn))**2-2*((un[2:,1:-1]-un[:-2,1:-1])/(2*dyn))*((vn[1:-1,2:]-vn[1:-1,:-2])/(2*dxn))-((vn[2:,1:-1]-vn[:-2,1:-1])/(2*dyn))**2)
    for _ in range(nit):
        pn=p.copy(); p[1:-1,1:-1]=((pn[1:-1,2:]+pn[1:-1,:-2])*dyn**2+(pn[2:,1:-1]+pn[:-2,1:-1])*dxn**2-b[1:-1,1:-1]*dxn**2*dyn**2)/(2*(dxn**2+dyn**2))
        p[:,-1]=p[:,-2]; p[:,0]=p[:,1]; p[0,:]=p[1,:]; p[-1,:]=0
    u[1:-1,1:-1]=(un[1:-1,1:-1]-un[1:-1,1:-1]*dt/dxn*(un[1:-1,1:-1]-un[1:-1,:-2])-vn[1:-1,1:-1]*dt/dyn*(un[1:-1,1:-1]-un[:-2,1:-1])-dt/(2*rho*dxn)*(p[1:-1,2:]-p[1:-1,:-2])+nu*dt*((un[1:-1,2:]-2*un[1:-1,1:-1]+un[1:-1,:-2])/dxn**2+(un[2:,1:-1]-2*un[1:-1,1:-1]+un[:-2,1:-1])/dyn**2))
    v[1:-1,1:-1]=(vn[1:-1,1:-1]-un[1:-1,1:-1]*dt/dxn*(vn[1:-1,1:-1]-vn[1:-1,:-2])-vn[1:-1,1:-1]*dt/dyn*(vn[1:-1,1:-1]-vn[:-2,1:-1])-dt/(2*rho*dyn)*(p[2:,1:-1]-p[:-2,1:-1])+nu*dt*((vn[1:-1,2:]-2*vn[1:-1,1:-1]+vn[1:-1,:-2])/dxn**2+(vn[2:,1:-1]-2*vn[1:-1,1:-1]+vn[:-2,1:-1])/dyn**2))
    u[0,:]=0; u[:,0]=0; u[:,-1]=0; u[-1,:]=1; v[0,:]=0; v[-1,:]=0; v[:,0]=0; v[:,-1]=0
xn=np.linspace(0,2,n); yn=np.linspace(0,2,n); speed=np.sqrt(u*u+v*v)
fig,ax=plt.subplots(); cf=ax.contourf(xn,yn,speed,levels=16,cmap="viridis"); ax.streamplot(xn,yn,u,v,color="white",density=1.1,linewidth=.7)
fig.colorbar(cf,ax=ax,label="Speed [a.u.]"); ax.set(title="Lid-driven cavity: recirculating flow",xlabel="x [m]",ylabel="y [m]")
ax.grid(True); fig.tight_layout(); plt.show()
print(f"Mean speed={speed.mean():.3f}, low-speed interior fraction={(speed[1:-1,1:-1] < 0.03).mean():.1%}")

Mean speed=0.117, low-speed interior fraction=28.3%

結果の読み取り

上壁の駆動で主循環が生じる一方、壁際や隅には低速域が残ります。炉なら低速域が熱伝達不足の候補です。ただし、この結果は無次元の教材モデルであり設備設計値ではありません。実務ではレイノルズ数、乱流モデル、入口条件、温度との連成を合わせ、流速測定や温度分布で妥当性を確認します。

No.089：境界条件

実務での意味

同じ方程式でも、表面を「炉温に固定」「一定熱流束」「対流で加熱」とするかで予測は大きく変わります。境界条件は計算設定ではなく、炉と製品の熱の受け渡しに関する業務仮説です。

分析・モデル化の考え方

代表例は、Dirichlet条件 $T=T_b$、Neumann条件 $-k\partial T/\partial n=q''$、Robin条件 $-k\partial T/\partial n=h(T-T_\infty)$ です。ここでは厚さ20 mmの平板を集中熱容量モデルで近似し、対流係数 $h$ の違いが昇温時間へ与える影響を比較します。

Pythonで確認する

rho_s, cp_s, thickness = 7800., 600., .020
tsec=np.linspace(0,3600,361)
rows=[]; fig,ax=plt.subplots()
for h in [25, 80, 200]:
    tau=rho_s*cp_s*thickness/(2*h)
    temp=T_furnace-(T_furnace-T0)*np.exp(-tsec/tau)
    hit=np.argmax(temp>=780) if np.any(temp>=780) else None
    rows.append({"h_W_m2K":h,"time_constant_min":tau/60,"time_to_780_min":None if hit is None else tsec[hit]/60})
    ax.plot(tsec/60,temp,label=f"h={h} W/m2K")
ax.axhline(780,color="black",ls="--",lw=1); ax.set(title="Heating sensitivity to convective boundary",xlabel="Time [min]",ylabel="Mean plate temperature [degC]")
ax.grid(True); ax.legend(); fig.tight_layout(); plt.show(); display(pd.DataFrame(rows).round(1))

	h_W_m2K	time_constant_min	time_to_780_min
0	25	31.2	NaN
1	80	9.8	24.2
2	200	3.9	9.7

結果の読み取り

対流係数の設定だけで規格到達時刻が大きく変わります。境界条件を推測値のまま精密なメッシュで解いても、精密な予測にはなりません。空運転の炉温だけでなく、代表ワークの温度履歴から $h$ を同定し、装入量やファン条件ごとに有効範囲を管理します。集中熱容量モデルは内部温度差が小さい場合に限る点にも注意します。

No.090：CFL条件

実務での意味

時間刻みが大きすぎると、現実には滑らかな温度や濃度が数値的に振動・発散します。計算が完了したことと、結果が信頼できることは別です。CFL条件は格子間を情報が伝わる速度と時間刻みの整合を確認する基準です。

分析・モデル化の考え方

一次元移流方程式 $\partial C/\partial t+u\partial C/\partial x=0$ の風上差分では、Courant数 $\mathrm{Co}=u\Delta t/\Delta x\le1$ が安定性の目安です。安定な $mathrm{Co}=0.8$ と不安定な $mathrm{Co}=1.2$ を比較します。

Pythonで確認する

na=101; xa=np.linspace(0,1,na); dxa=xa[1]-xa[0]; velocity=0.5
C0=np.exp(-((xa-.2)/.045)**2)
fig,ax=plt.subplots(); cfl_rows=[]
for Co in [.8,1.2]:
    dta=Co*dxa/velocity; C=C0.copy()
    nsteps=int(.9/dta)
    for _ in range(nsteps): C[1:]=C[1:]-Co*(C[1:]-C[:-1]); C[0]=0
    cfl_rows.append({"Courant":Co,"dt_s":dta,"steps":nsteps,"min":C.min(),"max":C.max(),"bounded_0_to_1":bool((C>=-1e-9).all() and (C<=1+1e-9).all())})
    ax.plot(xa,C,label=f"Co={Co}")
ax.set(title="Upwind advection: stable and unstable time steps",xlabel="Position x [m]",ylabel="Concentration [a.u.]")
ax.grid(True); ax.legend(); fig.tight_layout(); plt.show(); display(pd.DataFrame(cfl_rows).round(4))

	Courant	dt_s	steps	min	max	bounded_0_to_1
0	0.8	0.016	56	0.0000	0.7281	True
1	1.2	0.024	37	-0.4579	1.8461	False

結果の読み取り

$\mathrm{Co}=0.8$ ではピークが数値拡散で鈍るものの有界です。$\mathrm{Co}=1.2$ では負値や過大値が生じ、物理的に不合理な結果になります。実務では移流CFLだけでなく、熱拡散、反応、メッシュ最小幅に基づく制約を計算全域で監視し、時間刻みを変えた収束確認を行います。

対象ノックを通して見える実務上の示唆

1. PDEは測れない場所を埋める仮説である：温度場や流れ場を推定できますが、物性・初期条件・境界条件の妥当性が前提です。
2. 手法は形状と保存則で選ぶ：規則格子の試算はFDM、複雑形状はFEM、流体と収支重視はFVMが有力です。
3. 結果は品質KPIへ変換する：温度色図だけでなく、最低温度、最高温度、規格到達時刻、温度むら、低速域率で比較します。
4. 細かい計算ほど正しいとは限らない：格子・時間刻み収束に加え、境界条件の感度と実測との誤差を確認します。
5. 保存・有界性・残差を監視する：熱量収支、負温度や負濃度の発生、連続の式の誤差は、可視化より先に確認する品質指標です。

実務導入する場合に必要なこと

• 目的を「保持時間短縮」「温度むら低減」など検証可能なKPIへ落とす
• 材料物性、装入条件、炉壁・治具・ファン条件の版を管理する
• 熱電対、流速、消費電力など独立した実測データで校正・検証する
• 格子収束、時間刻み収束、収支誤差、パラメータ感度を記録する
• 適用範囲外の材質・形状・運転条件を検知し、再検証する
• 操業、品質、設備、解析担当がモデル変更と承認手順を共有する

まとめ

No.081〜No.090では、PDEの意味から三つの離散化、代表的な物理方程式、境界条件、CFL条件までを一つの熱処理課題で確認しました。数値解析の価値は美しい温度分布図ではなく、どの条件なら品質を満たすか、どの仮定が判断を左右するかを説明できることにあります。小さなモデルで収支と感度を理解し、実測で校正してから詳細CAEや運転最適化へ進むことが堅実です。

法人向けのご相談

数理工房では、製造業向けの数値シミュレーション、データ分析、モデル検証、技術研修、PoC設計をご支援します。現場データと物理モデルをどう結び付けるか、既存CAE結果を意思決定へどう翻訳するかといった段階からご相談いただけます。

📩 お問い合わせ: surikobo.co.jp/contact まずはお気軽にご相談ください。