Deciphering Temperature Fluctuations and Cooling Conditions in Heat Treatment Furnaces Using Partial Differential Equations: 10 Numerical Calculations in Manufacturing Industry (No.081–No.090)

This article uses Heating and cooling of steel plates and flow inside the furnace in a fictional metal parts factory as the subject and connects partial differential equations (PDEs) to determining manufacturing conditions. Check finite difference methods, finite element methods, finite volume methods, thermal conduction, wave analysis, Poissson-Navier–Stokes equations, boundary conditions, and CFL conditions using formulas, Python, tables, and graphs.

The aim is not to replace high-precision CAE, but to understand the model’s assumptions and judgment materials regarding temperature variability, retention time, sensor placement, and computational stability.

[!NOTE] This material is a notebook previously used by Surikoubo (or personally by the representative, Kazuyama), and has been reconstructed, edited, and published with the company’s permission. All data listed is fictional and has no relation whatsoever to real companies, factories, or figures.

Introduction: Practical Challenges in Manufacturing Covered in This Article

In heat treatment, even if the representative furnace temperature reaches the target, temperature differences may remain at the center and edges of the steel plate, as well as between the surface and interior. Temperature variation affects hardness, residual stress, and dimensional accuracy. On the other hand, extending the heating time toward safety increases the burden on electricity, processing capacity, and oxidation.

This article focuses on “which position enters the standard temperature and when,” “how to set the temperature at the boundary,” and “whether the time interval is disrupting the calculation results.”

Common situations on site

• The in-furnace thermocouple is the target temperature, but the core temperature of the product cannot be measured directly.
• Retention time is adjusted based on experience whenever material, board thickness, or loading amount changes
• Avoid overheating at the ends and insufficient heating at the center at the same time
• Although CAE results exist, it is difficult to explain differences caused by boundary conditions or meshes.

Why is this issue so difficult to judge?

Temperature and flow rate change not only over time but also in space. Measurement point values alone cannot guarantee the intervals between points; results also depend on physical properties such as thermal conductivity, convective heat transfer, and boundary conditions like fixed temperature and heat flux. Furthermore, discretized calculations have cutoff errors and stability conditions. Therefore, it is necessary to cross-check measurements, physical models, and numerical calculations.

Overview of Exercise covered this time

No.	Theme	Examples of judgments in manufacturing
081	What is a partial differential equation?	Definition of Quality Fields Including Time and Location
082	finite difference method	Calculation of temperature gradients and curvature on the lattice
083	finite element method	Incorporating shapes and material properties into matrices
084	finite volume method	Maintaining heat balance at the cellular level
085	Heat conduction equation	Evaluation of Holding Time and Temperature Variation
086	wave equation	Grasping the propagation time of shocks and vibrations
087	Poisson equation	Steady-state temperature distribution with internal heating
088	Navier–Stokes Equation	Understanding Internal Reactor Circulation and Retention Zones
089	boundary condition	Hypothetical comparison of furnace walls, convection, and insulation
090	CFL Conditions	Stable Time-Tracking Design

Preparing the Python environment

NumPy performs array calculation, pandas sets the table of decision indicators, SciPy sparse matrices, and matplotlib visualizes. No external data is used; random number seeds are fixed. Graph labels are displayed in English to avoid garbled text caused by environmental differences.

%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'}

Creation of Fictional Data

A 1.2-meter-long steel plate is one-dimensionalized and 21 points of temperature are considered. The reference physical properties are thermal diffusivity lpha=1.2 imes10^{-5}\ \mathrm{m^2/s}, initial temperature 25°C, and furnace temperature 850°C. In practice, temperature-dependent properties, plate thickness direction, jig contact, and radiation are added, but here we use constants to clarify the comparison of methods.

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: What is a Partial Differential Equation?

Meaning in Practice

PDEs describe “fields” dependent on multiple independent variables, such as temperature $T(x,t)$. It is a common language for considering the temperature history of the entire product from single-point furnace temperature monitoring.

Approach to Analysis and Modeling

The one-dimensional heat conduction equation is $\frac{\partial T}{\partial t}=\alpha\frac{\partial^2T}{\partial x^2}$ That’s right. The left side shows the temporal temperature change, and the right side shows the spatial temperature bend. The solution is only determined when the initial and boundary conditions are met. Here, we use an example of an analytical solution $T=T_0+A e^{-\alpha k^2t}\sin(kx)$ to confirm that it depends on both space and time.

Check with 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

Reading the results

Even at the same time, the temperature varies depending on the position, and this difference gradually diminishes over time. It is necessary not only to link the furnace temperature to “what temperature the furnace temperature is,” but also to the quality conditions of “which position of the product and what temperature after what minute.” This analytical solution is intended for simplified boundary confirmation and identifies physical properties and heat entry and exit for real-world reactor prediction.

No.082: Finite Difference Method

Meaning in Practice

The finite difference method (FDM) replaces temperature slopes and curvature on regular grids with differences in nearby points. It is suitable for quickly estimating the thermal history of simple shapes and investigating the effects of sensor spacing and calculation grids.

Approach to Analysis and Modeling

In the central difference $\frac{\partial^2T}{\partial x^2}\bigg|_i\approx\frac{T_{i+1}-2T_i+T_{i-1}}{\Delta x^2}.$ Finer lattices usually reduce discretization errors, but also increase computational complexity and narrow the allowable time intervals for explicit solutions.

Check with 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

Reading the results

If the grid width is halved, the error becomes roughly one-quarter, allowing confirmation of the secondary accuracy of the center difference. In practice, the narrowest grid is not considered the correct answer, but the main KPIs (maximum temperature, center temperature, time reaching standards) are confirmed to remain sufficiently stable due to grid subdivision.

No.083: Finite Element Method

Meaning in Practice

The Finite Element Method (FEM) divides shapes into elements and assembles the contributions of each element into a whole matrix. This is the foundation of CAE, which handles complex product shapes, holes, heterogeneous material joints, and local mesh subdivision.

Approach to Analysis and Modeling

In the linear elements of steady-state one-dimensional thermal conduction $-\frac{d}{dx}(k\frac{dT}{dx})=0$, the element stiffness matrix is $K_e=\frac{k}{h_e}\begin{bmatrix}1&-1\\-1&1\end{bmatrix}.$ By varying the thermal conductivity of each element, the thermal resistance of different materials can be expressed. Fixed temperatures on both sides and compare the temperature distribution between the base material and the low-thermal-conductivity insert.

Check with 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

Reading the results

The low thermal conductivity section causes a steep temperature gradient. Local thermal resistance that cannot be seen by average temperature alone can be checked by corresponding to material arrangements. In the production FEM, we verify element quality, contact thermal resistance, two-dimensional and three-dimensional shapes, and the temperature dependence of physical properties, and report mesh convergence.

No.084: Finite Volume Method

Meaning in Practice

The finite volume method (FVM) calculates the inflow and outflow to each cell as the balance of payments. It is easy to specify the conservation of heat, mass, and momentum, and is widely used for calculating in-furnace flow and cooling paths.

Approach to Analysis and Modeling

The heat balance of cell $i$ is $\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}.$ The heat flux between internal cells is stored because the amount emitted from one side enters the adjacent area, so the total heat is stored. Let’s look at an example where a single high-temperature cell diffuses at the insulated end.

Check with 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

Reading the results

The peak temperature decreases and spreads outward, and since it is an adiabatic boundary, the heat corresponding to the sum of all cell temperatures is stored within the range of numerical error. For CFD and thermal balance, acceptance criteria are based not only on color appearance but also on the balance errors at the entrance, outlet, wall, and accumulation.

No.085: Equation of Heat Conduction

Meaning in Practice

Using the heat conduction equation, it is possible to predict the time when the plate center will reach the lower limit of the standard and the temperature variation at that point. This forms the basis for quantifying the trade-offs of hold time, line speed, and energy.

Approach to Analysis and Modeling

Fix both ends to the furnace temperature, and the positive difference $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$ We use it. For one-dimensional thermal diffusion, $r\le 1/2$ is the benchmark for stability. Calculations are made until the core temperature reaches 780 °C.

Check with 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

Reading the results

The ends immediately reach furnace temperature due to boundary conditions, while it takes a long time to reach the central standard. If operating conditions are determined solely by the time the center reaches the center, overheating at the ends will be overlooked, so the minimum temperature, maximum temperature, and the specified residence time will be listed together. On the actual machine, calibrate using the radiation/convection boundary and the measured center temperature.

No.086: Wave Equation

Meaning in Practice

The wave equation represents the speed at which press shock, bar impact inspection, and equipment vibrations propagate and reflect back. This serves as the foundation for estimating abnormal positions based on sensor arrival time.

Approach to Analysis and Modeling

one-dimensional wave equation $\frac{\partial^2u}{\partial t^2}=c^2\frac{\partial^2u}{\partial x^2}$ is centered differentiated. $c\Delta t/\Delta x\le1$ is required for wave velocity $c$ and lattice conditions. We calculate how the central displacement pulse propagates left and right.

Check with 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

Reading the results

The initial pulse splits left and right and is reflected at the fixed end. The arrival time corresponds to distance ÷ wave speed, allowing you to narrow down the location of the occurrence based on the time differences between multiple sensors. However, since the actual material has dispersion, attenuation, cross-sectional changes, and sensor response, wave velocity and detection thresholds are calibrated during impact tests at known locations.

No.087: Poisson Equation

Meaning in Practice

Steady-state temperature with internal heat, electrostatic fields, and pressure correction all culminate in the Poisson equation. Components with heaters or reaction heat can be estimated for constant hotspot conditions.

Approach to Analysis and Modeling

Two-dimensional steady-state heat conduction $\nabla^2T=-q'''/k$ and fix the outer circumference temperature. Generate a sparse system of linear equations with a 5-point difference and compare cases where internal heat generation is uniform and skewed to the upper right.

Check with 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)}

Reading the results

The highest temperature position is not at the geometric center, but rather shifts to the side where heat is generated. If you place the sensor only in the center, you might miss the maximum value. Eliminate uncertainty in heat distribution and boundary temperature, and determine the placement by checking the range within which hotspot positions move.

No.088: Navier–Stokes Equation

Meaning in Practice

In furnace fans, cooling nozzles, and cleaning tanks, uneven flow leads to uneven heat and material transfer. The Navier–Stokes equation describes velocity and pressure fields and forms the basis for evaluating retention zones and circulation.

Approach to Analysis and Modeling

In non-compressible fluids, $\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.$ Here, for educational purposes, we will solve the square cavity with a moving upper wall using the pressure Poisson method. Although the grid is coarser and the repetition is shorter than in practical CFD, you can confirm the meaning of circulation and low-speed ranges.

Check with 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%

Reading the results

While the main circulation is generated by driving the upper wall, low-speed ranges remain along the walls and in corners. For furnaces, low-speed ranges are candidates for insufficient heat transfer. However, these results are dimensionless teaching material models and not equipment design values. In practice, validity is confirmed by combining Reynolds number, turbulence model, inlet conditions, and temperature coupling, using flow velocity measurements and temperature distribution.

No.089: Boundary Conditions

Meaning in Practice

Even with the same equation, predictions can vary greatly depending on whether the surface is “fixed at furnace temperature,” “constant heat flux,” or “heated by convection.” Boundary conditions are not calculation settings but operational hypotheses regarding heat transfer between the furnace and the product.

Approach to Analysis and Modeling

Representative examples include Dirichlet $T=T_b$, Neumann $-k\partial T/\partial n=q''$, and Robin $-k\partial T/\partial n=h(T-T_\infty)$. Here, a 20 mm thick plate is approximated using a concentrated heat capacity model, and the effect of the difference in convection coefficient $h$ on heating time is compared.

Check with 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

Reading the results

Simply setting the convection coefficient can greatly change the time it takes to reach the standard. Even if you solve boundary conditions with precise meshes while keeping the estimated values, you won’t achieve precise predictions. Not only the idle furnace temperature but also the temperature history of representative workpieces is used to identify $h$, and the effective range is managed according to loading volume and fan conditions. Note that concentrated heat capacity models are limited to cases where the internal temperature difference is small.

No.090: CFL Conditions

Meaning in Practice

If the time increments are too large, the actual smooth temperature and concentration will vibrate and dissipate numerically. Completing the calculations and the reliability of the results are separate matters. CFL conditions are criteria that confirm the consistency between the speed and timing of information transmission between grids.

Approach to Analysis and Modeling

In the windward difference of the one-dimensional advection equation $\partial C/\partial t+u\partial C/\partial x=0$, the Courant number $\mathrm{Co}=u\Delta t/\Delta x\le1$ This is the guideline for stability. Compare stable $mathrm{Co}=0.8$ with unstable $mathrm{Co}=1.2$.

Check with 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

Reading the results

In $\mathrm{Co}=0.8$, the peak is dulled by numerical diffusion but is bounded. $\mathrm{Co}=1.2$ negative or excessive values occur, resulting in physically irrational outcomes. In practice, not only advection CFL but also thermal diffusion, reaction, and constraints based on mesh minimum width are monitored across the entire calculation range, and convergence confirmation is performed at different time intervals.

Practical Implications Seen Through Target Exercise

1. PDEis a hypothesis that fills in places that cannot be measured.: It is possible to estimate temperature and flow fields, but the validity of physical properties, initial conditions, and boundary conditions is a prerequisite.
2. Choose methods based on shape and preservation rules: FDM is the most effective method for estimating regular grids, FEM for complex shapes, and FVM for fluid and balance considerations.
3. The result is quality.KPIconvert to: Compare not only temperature color diagrams but also minimum and maximum temperatures, time reaching standards, temperature unevenness, and low-speed frequency range.
4. Fine-grained calculations are not always correct: In addition to lattice and time-scale convergence, we check the sensitivity of boundary conditions and the error between the measured measurements.
5. Monitoring conservation, boundality, and residuals: Calorific balance, occurrence of negative temperatures or concentrations, and errors in continuous formulas are quality indicators to be checked before visualization.

What is necessary for practical implementation

• Objectives are set to verifiable KPIs such as “shortening retention time” and “reducing temperature variability.”
• Managing plates for material properties, loading conditions, furnace walls, jigs, and fan conditions
• Calibration and verification using independent measured data such as thermocouples, flow velocity, and power consumption
• Records lattice convergence, time-step convergence, balance error, and parameter sensitivity
• Detects and re-examines materials, shapes, and operating conditions outside the applicable scope
• Operational, quality, equipment, and analysis teams share model changes and approval procedures.

Conclusion

From No.081 to No.090, we examined the meaning of PDE, three discretizations, representative physical equations, boundary conditions, and CFL conditions under a single heat treatment problem. The value of numerical analysis lies not in a beautiful temperature distribution chart, but in explaining which conditions meet quality and which assumptions influence decisions. It is reliable to understand the balance and sensitivity with small models, calibrate them through actual measurements, and then proceed to detailed CAE and operational optimization.

Consultations for Corporations

At Surikoubo, we support numerical simulation, data analysis, model validation, technical training, and PoC design for manufacturing industries. You can consult with us from the stage on how to link on-site data with physical models or how to translate existing CAE results into decision-making.

📩 Contact Us: surikobo.co.jp/contact Please feel free to consult us first.