2. Heat transfer fundamentals for coupled analysis

Thermal–structural coupling begins with an accurate thermal solution. Before any displacement, strain, or stress can be computed, the temperature field must be established by solving the energy balance across the domain. In finite element analysis (FEA), this means assembling and solving the heat equation with appropriate material properties, source terms, and boundary conditions. This chapter covers the three fundamental heat transfer mechanisms — conduction, convection, and radiation — and the thermal boundary conditions that close the governing equations. The emphasis is on how each mechanism is represented in a coupled FEA workflow and what inputs you must supply.

2.1 The Governing Energy Equation

For a stationary solid with no phase change, the steady-state heat conduction equation is:

∇ · (k ∇T) + Q = 0

where k is the thermal conductivity tensor [W/(m·K)], T is temperature [K or °C], and Q is the volumetric heat generation rate [W/m³]. For transient analysis, the equation becomes:

ρ c_p ∂T/∂t = ∇ · (k ∇T) + Q

with density ρ [kg/m³] and specific heat capacity c_p [J/(kg·K)]. In a coupled analysis, the temperature field from this equation becomes the thermal load for the structural solver. The coupling can be one-way (thermal results drive structural deformation) or two-way (deformation feeds back into the thermal solution through contact conductance or geometry change).

2.2 Conduction

Conduction is the dominant heat transfer mode within solid bodies. Fourier's law defines the local heat flux:

q = −k ∇T

The negative sign indicates that heat flows from high to low temperature — this sign convention must be consistent throughout your model. In FEA, conduction is captured by the element conductivity matrix, which depends on the mesh and the material's k value. For isotropic materials, k is a scalar; for orthotropic materials (e.g., laminated composites, wood), k becomes a 3×3 tensor with directional values k_x, k_y, k_z.

Temperature-dependent conductivity is common in engineering materials. For example, carbon steel at 20 °C has k ≈ 52 W/(m·K), dropping to ≈ 30 W/(m·K) at 600 °C. If your thermal gradient spans several hundred degrees, a constant-k assumption may introduce significant error. Most FEA solvers allow tabular k(T) input, which triggers a nonlinear iteration.

2.3 Convection

Convection transfers heat between a solid surface and an adjacent fluid. In FEA, it is applied as a boundary condition rather than solved within the solid domain (unless you are running a conjugate heat transfer simulation). Newton's law of cooling gives the convective flux:

q_conv = h (T_s − T_∞)

where h is the convective heat transfer coefficient [W/(m²·K)], T_s is the surface temperature, and T_∞ is the reference fluid temperature. The sign convention here is that positive q_conv means heat leaves the surface when T_s > T_∞.

Selecting h is often the most uncertain step. Typical ranges:

These are order-of-magnitude guides. For critical applications, h should be derived from empirical correlations (e.g., Dittus–Boelter for internal pipe flow, Churchill–Chu for external flow over cylinders) or calibrated against measured temperature data.

2.4 Radiation

Thermal radiation becomes significant when temperatures exceed roughly 300 °C, or when surfaces exchange heat across a vacuum or transparent medium. The net radiative flux between a small gray surface and its surroundings is:

q_rad = ε σ (T_s⁴ − T_sur⁴)

where ε is the surface emissivity (0 ≤ ε ≤ 1), σ = 5.67×10⁻⁸ W/(m²·K⁴) is the Stefan–Boltzmann constant, and T_sur is the surrounding temperature. Because radiation scales with the fourth power of absolute temperature, it is highly nonlinear. In FEA, radiation boundary conditions require an iterative solution even for steady-state problems.

For enclosure problems — where multiple surfaces radiate to each other — view factors F_ij must be computed. The net heat transfer to surface i is:

q_i = Σ_j ε_i σ F_ij (T_i⁴ − T_j⁴)

View factor calculation is geometry-dependent and can be expensive for complex 3D models. Many solvers offer a radiosity solver or a Monte Carlo ray-tracing option. For simple configurations (parallel plates, concentric cylinders), analytical view factors are available and should be used to validate the solver's output.

2.5 Thermal Boundary Conditions

A well-posed thermal problem requires at least one boundary condition on every boundary. The four common types are:

• Prescribed temperature (Dirichlet): T = T₀ on the boundary. Used when a surface is clamped to a known temperature (e.g., a heat sink base).
• Prescribed heat flux (Neumann): −k ∂T/∂n = q₀. Used for applied heating (e.g., laser irradiation, resistive heating).
• Convection: −k ∂T/∂n = h(T_s − T_∞). The most common boundary condition in thermal FEA.
• Radiation: −k ∂T/∂n = εσ(T_s⁴ − T_sur⁴). Often combined with convection.

Adiabatic (insulated) boundaries are the natural boundary condition — if no condition is specified, the solver assumes zero flux. This is rarely physical and should be applied intentionally only when justified.

2.6 Coupling to Structural Analysis

Once the temperature field is solved, it is passed to the structural solver as a thermal load. The thermal strain is:

ε_th = α (T − T_ref)

where α is the coefficient of thermal expansion (CTE) and T_ref is the stress-free reference temperature. The choice of T_ref is critical: it should correspond to the temperature at which the part is stress-free, often the manufacturing or assembly temperature. If T_ref is set incorrectly, the entire structure will show a uniform expansion or contraction that does not reflect reality.

For one-way coupling, the thermal solution is computed first, then the structural solution uses the temperature field as input. For two-way coupling, the structural deformation updates the thermal boundary conditions (e.g., contact conductance changes with pressure), and the cycle repeats until convergence.

2.7 Practical Workflow and Verification

A reproducible thermal FEA workflow for coupled analysis follows these steps:

1. Define the thermal problem: Identify the domain, material properties (k, ρ, c_p), and heat sources. Determine whether the analysis is steady-state or transient.
2. Assign boundary conditions: For each external surface, specify the appropriate condition type. Document the source of each h, ε, or q value.
3. Mesh and solve: Use a mesh fine enough to resolve thermal gradients. For radiation-dominated problems, ensure the mesh captures view factor geometry.
4. Verify the thermal solution: Check the energy balance — total heat in should equal total heat out at steady state. Inspect temperature contours for physical plausibility.
5. Pass to structural solver: Map the temperature field to the structural mesh (if different). Set T_ref correctly.
6. Validate: Compare predicted temperatures against thermocouple data or analytical solutions for simplified cases.

A common verification check is the Biot number, Bi = hL/k, which indicates whether internal conduction or surface convection dominates. For Bi < 0.1, the body can be treated as isothermal (lumped capacitance). For Bi > 0.1, spatial temperature gradients are significant and must be resolved by the mesh.