Thermal Dunking Problem Modeling

• Thermal Dunking is the study of transient heat transfer when a solid is suddenly immersed in a media at a different temperature, revealing conduction, boiling, and kinetic phenomena.
• The analysis employs lumped-capacitance models, eigenvalue corrections, and error quantification to predict temperature evolution under low Biot-number conditions.
• Practical implications include safe quenching, spacecraft thermal design, and accurate calibration of conduction and phase-change models across diverse environments.

The thermal dunking problem concerns the transient heat transfer and fluid/phase-change phenomena arising when a solid at one temperature is suddenly immersed (“dunked”) in an environment—air, liquid, or rarefied gas—at another temperature. Dunking protocols appear in thermal management, safety-critical quenching (e.g., cryostats, heat exchangers), spacecraft thermal design, and calibration of conduction/boiling/creep models. Technical analyses focus on the conduction-dominated regime (small Biot number), the boiling transition in liquid environments (critical heat flux and film boiling), and kinetic (thermal-slip) effects in dilute gases.

1. Mathematical Formulation of the Dunking Problem

The core mathematical abstraction is a solid body $\Omega \subset \mathbb{R}^d$ with possibly heterogeneous material properties—spatially varying density $\rho$, specific heat $c$, and thermal conductivity $k$—initially at $T_\text{i}$. For $t > 0$, the body is exposed to an environment at farfield temperature $T_\infty$ with uniform heat transfer coefficient $h$ (encompassing fluid convection or external coupling). The initial-boundary-value problem is governed by the heat equation with Robin boundary condition: $$\rho c \, \partial_t T(x,t) = \nabla \cdot (k(x) \nabla T(x,t)), \quad x \in \Omega, \ t>0$$

$$k(x)\, \partial_n T(x,t) + h(T(x,t)-T_{\infty}) = 0, \quad x \in \partial\Omega$$

where $\partial_n$ denotes the outward normal derivative.

Nondimensionalization introduces the extrinsic scale $\ell$ and the Biot number $B = h \ell/k_{\mathrm{inf}}$ with $k_{\mathrm{inf}} = \inf_{x} k(x)$. The intrinsic length scale $L = |\Omega|/|\partial\Omega|$ yields the intrinsic Biot number $\mathrm{Bi} = h L/ k_{\mathrm{inf}} = B/\gamma$ where $\gamma = |\partial\Omega|/|\Omega|$. The classical regime of interest is small Biot number ($\mathrm{Bi} \ll 1$), in which internal temperature gradients are negligible and lumped models are accurate.

2. Lumped-Capacitance and Higher-Order Models: Error Analysis

The lumped-capacitance (LCM) approximation assumes a spatially uniform temperature in the solid. Integrating the PDE system, the domain-average temperature $u_{\mathrm{avg}}(t)$ evolves as: $$\dfrac{d u_{\mathrm{avg}}}{d t} = -B \gamma u_{\mathrm{avg}}, \quad u_{\mathrm{avg}}(0) = 1$$ with solution

$u_{\mathrm{avg}}(t) = \exp(-B \gamma t)$

or $u_{\mathrm{avg}}(t) = \exp(-\mathrm{Bi} \, t)$.

For increased accuracy, a second-order Padé-type correction leverages Biot-sensitivity analysis. The leading eigenvalue of the associated Robin problem expands as

$\lambda_1(B) = B \gamma - \varphi B^2 + O(B^3)$

where $\varphi$ is a quadratic functional, computed as $\varphi = \int_\Omega \kappa |\nabla \xi|^2\,dx$ for the sensitivity field $\xi$ solving: $$-\nabla \cdot (\kappa \nabla \xi) = |\Omega|^{-1/2} \gamma \sigma, \qquad \kappa \partial_n \xi = -|\Omega|^{-1/2}, \quad \int_\Omega \sigma \xi = 0$$ This leads to the corrected evolution: $$u_{\mathrm{avg}}^{(2)}(t) = \exp\left(-\frac{B \gamma t}{1 + \varphi B / \gamma}\right)$$ Error analysis yields the following bounds:

• Asymptotic error for LCM: $$|u(T;B) - u_{\mathrm{avg}}^{(1)}(T)| \leq \varphi B/( \gamma e ) + O(B^2)$$
• Strict bound (all $B$): $$|u(T;B) - u_{\mathrm{avg}}^{(1)}(T)| \leq \frac{1}{2} \sqrt{\varphi B/ \gamma}$$
• Second-order approximation error: $O(B^2)$, with coefficient involving additional functionals $\chi,\, \Upsilon$.

The functional $\varphi$ imparts a domain-dependent correction to the conduction length, with canonical values: slab $1/3$, disk $1/2$, sphere $3/5$, and can be much larger (slender or nonconvex domains), invalidating the standard $\mathrm{Bi} \ll 1$ criterion. Empirical and computational investigations (high-$\varphi$ domains) confirm sharpness of these estimators (Kaneko et al., 17 Jun 2024, Kaneko et al., 20 Dec 2024).

3. Hierarchical Modeling: From Conjugate Heat Transfer to Lumped ODE

Realistic dunking protocols require accounting for conjugate heat transfer (CHT): fluid convection (e.g., Navier–Stokes and energy equations) coupled to solid conduction at the interface. In dimensionless form for a fluid/solid pair with interface $\Gamma$:

• Fluid: $\partial_t f + \ldots = \nabla^2 f$ (advection-diffusion)
• Solid: $$\sigma\partial_t s = \nabla \cdot (\kappa \nabla s)$$
• Interface: temperature and flux continuity.

Reduction to LCM is justified under two conditions: (i) spatial uniformity of solid temperature (small Biot, negligible internal gradients), and (ii) time scale separation, with the fluid equilibrating much faster than the solid. The error decomposition (Guo et al., 10 Nov 2025) is: $$|y(t) - u(t; \tilde{B})| \leq \underbrace{|y(t; \eta(t,x)) - u(t; \bar{\eta}(x))|}_{\text{Time-homogenization error}} + \underbrace{|u(t; \bar{\eta}(x)) - u(t; B)|}_{\text{Lumping error}} + \underbrace{|u(t; B) - u(t; \tilde{B})|}_{\text{Biot estimation error}}$$ Lumping error is bounded by $\varphi/( \gamma e ) B$, consistent with certified estimates (Kaneko et al., 17 Jun 2024, Kaneko et al., 20 Dec 2024). Time homogenization error scales as $O(R^{-1})$, $R = \tau_e / \tau_c$ (solid/fluid time constants), and for practical cases (e.g., Reynolds numbers up to $10^4$), direct simulation validates these bounds. Computable upper bounds for $\varphi$ involve only geometry and measurable variances: $$\varphi(\kappa, \sigma, \eta) \leq [ \sqrt{\varphi(1,1,1)} + \sqrt{ \delta^\sigma } + \sqrt{ \delta^\eta } ]^2$$ with $\delta^\sigma, \delta^\eta$ variances and geometric constants.

4. Dunking in Liquids: Boiling Crisis, Drying Transition, and Critical Heat Flux

When dunking occurs into a boiling liquid, heat fluxes may reach the “boiling crisis” threshold (critical heat flux, CHF) at which a quasi-continuous vapor film forms, suppressing further heat transfer (film boiling). The crisis is governed by the non-equilibrium drying (wetting) transition (Nikolayev et al., 2016): vapor recoil forces near a nucleated bubble contact line impart sufficient momentum to modify the apparent contact angle, $u(0^+)= \theta_{\mathrm{app}}>0$, causing the contact line to depin and the dry patch to grow. The CHF is attained when the nondimensional recoil strength $N_r \simeq 1$, $N_r \propto q_S^2$ (local heat flux). The scaling for $q_{\mathrm{CHF}}$ is: $$q_{\mathrm{CHF}} \approx H \sqrt{ \frac{\rho_V \sigma}{R} }$$ where $H$ is latent heat, $\rho_V$ vapor density, $\sigma$ surface tension, $R$ bubble radius. Empirical and numerical results confirm this recoil-induced criterion and its predictive capability for safe dunking protocols. Thus, maintaining $N_r < 1$ ensures continued nucleate boiling (Nikolayev et al., 2016).

5. Rarefied Gas Dunking: Thermal-Slip Phenomena

In rarefied gas environments, thermal sliding or creep occurs along the plate when a tangential temperature gradient is imposed. The BGK kinetic equation with Maxwell mirror-diffuse boundary conditions describes the thermal-slip velocity $U_{\mathrm{sl}}$ (Latyshev et al., 2016): $$\mu \frac{\partial h}{\partial x} + h(x, \mu) = \frac{1}{\sqrt{\pi}} \int e^{-t^2} h(x, t) dt + G_T (\mu^2 - \tfrac{5}{2})$$ where $G_T = d \ln T / dy$ is the temperature gradient, and $h(x,\mu)$ the perturbation of the Maxwellian. The boundary condition reads: $h(0, \mu) = (1 - \alpha) h(0, -\mu), \quad \mu > 0$ with accommodation coefficient $\alpha$, $0 \leq \alpha \leq 1$. The method of a source leads to a Fredholm equation for the spectral density, yielding the closed-form slip velocity expansion: $$U_{\mathrm{sl}}(\alpha) = G_T \left( 0.5 + 0.28566\alpha - 0.021192\alpha^2 \right )$$ In dimensional notation: $$u_{\mathrm{sl}}(\alpha) = K(\alpha) \nu \frac{d\ln T}{dy}, \quad K(\alpha) = 0.75(1+0.5714\alpha - 0.0422\alpha^2)$$ The expansion yields sub-percent accuracy for all $\alpha$, with robust agreement with classical theory.

6. Computational Methods and Practical Guidance

Evaluation of certified error bounds and sensitive functional outputs requires (i) adaptive finite-element solution of the sensitivity PDE for $\xi$, (ii) integration of functionals $\varphi, \chi, \Upsilon$, and (iii) direct numerical evolution of the full heat equation. For realistic domains, one can estimate geometry-dependent corrections using known values for canonical shapes, or via eigenvalue computation for arbitrary $\Omega$. The practical guideline, validated across direct simulations up to $\mathrm{Re} = 10^4$ (Guo et al., 10 Nov 2025), prescribes:

• Estimate Biot number via ISO/correlation or empirical data.
• Compute geometry constants and property variances.
• Bound lumped, time-homogenization, and Biot errors.
• If all contributions are below design tolerances, LCM/ISO models are justified; otherwise, employ higher-fidelity CHT or solve full heat equation.

7. Domain Geometry and Validity of Lumped Models

A key insight from recent analyses is the deficiency of the standard criterion $\mathrm{Bi} \ll 1$ in predicting lumped-model accuracy. The relevant correction is $\varphi \cdot \mathrm{Bi}$; domains with geometric “amplifiers” (e.g., slender fins, gear shapes, sharp angles) yield $\varphi \gg 1$, invalidating the classical guideline. Empirical stability of $\varphi$ under small domain perturbations enables surrogate-based estimation for nearby shapes (Kaneko et al., 17 Jun 2024, Kaneko et al., 20 Dec 2024).

Conclusion

The thermal dunking problem encompasses a wide range of coupled conduction, convection, boiling, and kinetic phenomena. Modern theoretical and computational frameworks provide certified bounds and error predictors for lumped and semi-lumped models, reveal the significance of geometry-corrected conduction lengths, and establish the pivotal role of phase-change thresholds in boiling scenarios. The domain functional $\varphi$ replaces the textbook Biot-number criterion, ensuring that approximation accuracy is a computable, geometry- and property-dependent quantity. These advances enable rigorous design and safety assessment for both classical and extreme-geometry dunking protocols.