多次元非定常熱伝導

Good question. Actually, if the boundary conditions in each direction are independent in an orthogonal coordinate system, the "product solution" can be used. For example, the temperature of a short rectangular bar (like a brick shape) is

It can be found simply by multiplying the 1D solution for each direction. You can obtain a 3D unsteady temperature field by hand calculation without FEM.

The key is the separation of variables method. The multidimensional heat conduction equation can be decomposed into 1D problems for each direction if the boundary conditions are separable in each axial direction. Mathematically, the partial differential equation reduces to a product of ordinary differential equations for each variable. For example, consider the quenching cooling of an automotive forging (short cylindrical shape). Treat the side surface as a cylinder and the top/bottom surfaces as plates, and multiply their respective 1D solutions.

Exactly right. The product solution is extremely useful as a benchmark (verification solution) for FEM. You can quantitatively evaluate the validity of the mesh and time step by comparing it with an exact analytical solution. Whether you know this in practice makes a huge difference in the confidence level of your analysis results.

The starting point is the multidimensional version of Fourier's heat conduction equation. For an isotropic, homogeneous medium with no internal heat generation:

Here $\alpha = k/(\rho c_p)$ is the thermal diffusivity [m²/s]. This equation is a parabolic partial differential equation; a unique solution is determined by giving the initial condition $T(\mathbf{x}, 0) = T_i$ and boundary conditions.

In cylindrical coordinates $(r, \phi, z)$:

For axisymmetric cases (uniform in the $\phi$ direction) with finite length, it becomes a 2D problem in the $r$ and $z$ directions. This is exactly where the product solution comes into play.

The core lies in the conditions for separation of variables to hold. Assuming a solution of the form $\theta^*(\mathbf{x},t) = X(x,t) \cdot Y(y,t) \cdot Z(z,t)$ and substituting it into the governing equation yields:

If each term depends only on its respective variable, it separates into 1D problems for each direction. The resulting $X$, $Y$, $Z$ become the solutions to the 1D unsteady heat conduction (Heisler solution) for each direction.

Let me list three typical cases:

1. Short Rectangular Bar — Product of 3 directional plate solutions

2. Short Cylinder — Cylinder solution × Plate solution

3. Semi-infinite Solid Corner — Product of 3 directional semi-infinite solid solutions

For example, when considering reflow heating of an electronic component's resin package (rectangular parallelepiped), case 1 can be used directly.

That's an important point. The product solution is not applicable in the following cases:

• When there is internal heat generation ($\dot{q} \neq 0$) — The equation becomes non-homogeneous and variables cannot be separated.
• When material properties are temperature-dependent — The equation becomes nonlinear.
• When boundary conditions are coupled between directions (oblique boundaries, spatial distribution of contact resistance, etc.)
• For shapes with non-orthogonal coordinate systems (L-shape, T-shape, etc.)

For these cases, numerical methods (FEM, FDM, FVM) are necessary. However, many practical verification problems are set up for shapes that can be handled by the product solution, so knowing this method is very valuable.

Yes, the superposition principle. For linear partial differential equations, you can add individual solutions to create a new solution. This is a powerful technique separate from the product solution.

For example, even with asymmetric boundary conditions like one surface being rapidly heated while others are insulated, you can decompose it into symmetric and antisymmetric problems and add them. Single-sided quenching of steel is exactly this pattern.

To summarize:

Remember: product solution is for "decomposition of spatial directions," superposition is for "decomposition of boundary conditions." By combining both, you can handle fairly complex problems analytically.

The shape factor method is primarily used for steady-state multidimensional heat conduction. Defining the conduction shape factor $S$:

Here $Q$ is the heat flow rate per unit time [W]. For example, heat dissipation from a pipe buried underground, or the heat path from a semiconductor chip to a substrate—you can calculate the heat flow simply by looking up the value of $S$ for standard shapes for which analytical solutions exist from a table.

Here are some representative ones:

Incropera's heat transfer textbook lists dozens of shape factor patterns. In practice, it's common to use this for rough estimates in the initial stage and then move to FEM for detailed design.

Exactly. The multidimensional conductive thermal resistance can be expressed as $R_{\text{cond}} = 1/(kS)$. By combining this in series/parallel with convective resistance $R_{\text{conv}} = 1/(hA)$, you can evaluate the total thermal resistance of a complex system by hand calculation using circuit analogy. This is a technique used daily in the thermal design of electronic devices.

Let me explain the procedure using the example of a short cylinder (radius $r_0$, half-length $L$). For the case where all surfaces are cooled by convection ($h$, $T_\infty$):

1. Step 1: For the $r$ direction, calculate Bi$_{\text{cyl}} = hr_0/k$, Fo$_{\text{cyl}} = \alpha t/r_0^2$, and read the center temperature $\theta^*_{\text{cyl}}(0,t)$ from the infinite cylinder Heisler chart.
2. Step 2: For the $z$ direction, calculate Bi$_{\text{wall}} = hL/k$, Fo$_{\text{wall}} = \alpha t/L^2$, and read the midplane temperature $\theta^*_{\text{wall}}(0,t)$ from the infinite plate Heisler chart.
3. Step 3: The center temperature of the short cylinder is found by the product:

If surface temperature or temperature at arbitrary locations is needed, you also use Heisler's second chart (position correction chart) to apply location-specific corrections.

Let's do a specific example. SUS304 short cylinder ($r_0 = 25$ mm, $L = 50$ mm = half-length):

• $k = 14.9$ W/(m-K), $\alpha = 3.95 \times 10^{-6}$ m²/s
• $h = 200$ W/(m²-K), $T_i = 500$°C, $T_\infty = 25$°C
• Desired time: $t = 300$ s

$r$ direction: Bi$_{\text{cyl}} = 200 \times 0.025 / 14.9 = 0.336$, Fo$_{\text{cyl}} = 3.95 \times 10^{-6} \times 300 / 0.025^2 = 1.896$
From the chart: $\theta^*_{\text{cyl}}(0,t) \approx 0.23$

$z$ direction: Bi$_{\text{wall}} = 200 \times 0.05 / 14.9 = 0.671$, Fo$_{\text{wall}} = 3.95 \times 10^{-6} \times 300 / 0.05^2 = 0.474$
From the chart: $\theta^*_{\text{wall}}(0,t) \approx 0.61$

Center temperature: $\theta^* = 0.23 \times 0.61 = 0.140$
$T_{\text{center}} = T_\infty + \theta^*(T_i - T_\infty) = 25 + 0.140 \times 475 = 91.6$°C

The center temperature after 300 seconds is about 92°C. By solving the same conditions with FEM and comparing with this value, you can verify the validity of the mesh and time step.

A key point in FEM formulation is the separation of spatial discretization and temporal discretization. Applying the weak form (Galerkin method) in space yields a semi-discretized system of ordinary differential equations:

Here $[\mathbf{C}]$ is the capacitance matrix (thermal version of mass matrix), $[\mathbf{K}]$ is the conduction matrix (thermal version of stiffness matrix), and $\{F\}$ is the external heat load vector.

For the time direction, time integration schemes like Euler methods (forward/backward) or the Crank-Nicolson method are applied. For 2D, triangular or quadrilateral elements are used; for 3D, tetrahedron or hexahedron elements are used.

It can be written uniformly using the generalized $\theta$ method (here $\theta$ is a parameter, distinct from dimensionless temperature):

In practice, the Galerkin method ($\theta=2/3$) is often used as it balances stability and accuracy well.