Thermal Spallation Simulation

M1 focuses on transient heat transfer from a continuous-wave Gaussian laser beam onto a rock half-space (Bäuerle 2011). M2 focuses on the quasi-static plane-strain thermoelastic response of a uniformly heated rock with an edge crack, as modelled by Kant et al. (2017). M3 combines the two — a localised Gaussian laser drives a transient temperature field that produces a non-uniform hoop stress along an edge crack, and the simulation runs until the stress intensity factor at the crack tip reaches the rock's fracture toughness, marking the onset of thermal spallation.

Setup

mamba env create -f environment.yml
mamba activate thermal_spallation_simulation
jupyter notebook notebooks/m1_prototype.ipynb

M1 — Axisymmetric Laser Heating

Transient heat conduction on a rock half-space in axisymmetric $(r, z)$ coordinates. A Gaussian Neumann flux is applied on the irradiated surface:

$$-k\frac{\partial T}{\partial z}\bigg|_{z=0} = I_0, e^{-2r^2/R^2}, \qquad I_0 = \frac{2P_t}{\pi R^2}$$

The heat equation $\rho C_v \partial T/\partial t = \nabla \cdot (k\nabla T)$ is discretised with backward Euler in time and Lagrange-1 finite elements (r-weighted weak form to account for the cylindrical Jacobian). The FEM solution is validated at the beam centre against the analytical centreline temperature

$$T(0,0,t) = T_i + \frac{I_0 w_0}{k\sqrt{\pi}}, \arctan!\left(\frac{\sqrt{4\alpha t}}{w_0}\right), \qquad w_0 = R/\sqrt{2}$$

(Bäuerle 2011, Eq 7.5.2, surface-absorption limit). Validation is not restricted to the centreline: the temperature field is sampled across the entire domain on a 30 × 30 probe grid and compared pointwise against the full analytical solution (Bäuerle 2011, Eq 6.2.2). At $t = 60$ s the L2 error over all probes is $0.04$ °C and the maximum pointwise error is $0.81$ °C — the FEM solution reproduces the analytical prediction across the whole domain.

Prototype: notebooks/m1_prototype.ipynb.

M1 Validation. Left: centre temperature (analytical vs FEM). Middle: absolute error convergence. Right: L2 and max error across full 30×30 probe grid. Bottom: 2D temperature field at t=60s (FEM, analytical, pointwise difference). Max error 0.8°C over entire domain.

M2 — Plane-Strain Thermoelasticity

Quasi-static plane strain under a uniform temperature jump $\Theta = T - T_i$. The displacement field $u_i(x,z)$ satisfies $\partial_j \sigma_{ij} = 0$ with the Duhamel–Neumann constitutive law

$$\sigma_{ij} = \lambda\varepsilon_{kk}\delta_{ij} + 2\mu\varepsilon_{ij}- \frac{E\alpha}{1-2\nu},\Theta,\delta_{ij} $$

(Landau & Lifshitz, Theory of Elasticity, §6). The computational domain is a Cartesian half-strip with kinematic constraints: $u_x = 0$ on the lateral boundaries, $u_z = 0$ at the base, and a traction-free top surface.

The FEM solution is validated three ways against Kant et al. (2017):

1. $\sigma_{xx}$ is uniform to machine precision under uniform $\Theta$,
2. $\sigma_{xx} = -E\alpha\Theta/(1-\nu)$ (Kant Eq 5) matched to $4.44\times10^{-12}$,
3. the onset $\Theta$ for $K_I \geq K_{IC}$ (Kant Eq 15) matched to $2\times10^{-9}$ °C.

Prototype: notebooks/m2_prototype.ipynb.

M2.5 rebuilds the same physics in cylindrical $(r, z)$ coordinates with the hoop stress $\sigma_{\theta\theta}$ in the role Kant assigns to $\sigma_{xx}$. This is the formulation that generalises to non-uniform heating: Kant's "plane strain" is really a uniform lateral-confinement assumption ($\varepsilon_{xx} = \varepsilon_{yy} = 0$ around the crack), which only survives a Gaussian laser spot under cylindrical symmetry. The three Kant checks pass to $\sim 10^{-10}$ ($\sigma_{\theta\theta}$ vs Eq 5) and $1.5\times 10^{-7}$ °C (onset vs Eq 15). Prototype: notebooks/m2_5_prototype.ipynb.

M3 — Localised Gaussian Heating, Coupled Thermal–Elastic

M3 couples M1's transient axisymmetric thermal solver with M2.5's axisymmetric thermoelastic solver on a single graded mesh: 2 µm cells inside a $1,\text{mm} \times 30,\mu\text{m}$ box covering the edge crack at the origin, smoothly ramping to $\sim 3,\text{mm}$ at the outer boundary of a $150 \times 150$ mm domain. The coupling is one-way and quasi-static: at each timestep $\Delta t = 0.1$ s, the thermal solver advances $T(r, z, t)$, the elastic solver consumes $\Theta(r, z, t) = T - T_i$ as a fem.Function load, and the stress intensity factor

$$K_I(t) = 2\sqrt{\frac{a}{\pi}}\int_0^1 \left[-\sigma_{\theta\theta}(0, a\xi, t)\right] f_K(\xi),d\xi$$

is evaluated by reusing analytical.stress_intensity_factor unchanged from M2 — the integrand $\sigma_{\theta\theta}(0, z, t)$ is sampled from the FEM solution via cell lookup plus point evaluation. The simulation halts at the first $t^$ where $K_I(t^) \geq K_{IC}$.

For the granite baseline (300 W, 25 mm beam radius), crack length $a = 20$ µm, the simulation reports $t^ = 39.3$ s with surface temperature $T(0, 0, t^*) = 827$ °C*.

The result encodes the physics that distinguishes localised heating from Kant's uniform model. The 827 °C surface temperature is far above Kant's uniform-$\Theta$ onset of 446 °C, but $\sigma_{\theta\theta}$ at the crack ($\approx -169$ MPa) precisely matches what Kant Eq 5 predicts for $\Theta = 446$ °C — only $\sim 55%$ of what the formula would give for the local temperature at the crack itself. The cold deep rock relaxes the lateral confinement that Kant's uniform-slab geometry treats as perfect, so a hotter surface is needed to drive the crack open by the same amount. M3 captures this; M2.5 cannot.

The approach was to validate each layer separately before coupling — the axisymmetric formulation was first introduced as M2.5 and matched against the M2 Cartesian result, which gave confidence that M3 only had to add the non-uniform $\Theta$ load on top of an already-validated solver. The stress intensity weight function $f_K$ in analytical.stress_intensity_factor did not need to change; it accepts any depth-dependent stress callable, and the edge-crack form is sufficient for the axisymmetric radial-crack geometry at the order of approximation Kant himself uses.

Prototype: notebooks/m3_prototype.ipynb.

Rock Properties (Granite Baseline)

Property	Value	Source
Density	2250 kg/m³	Han et al. (2017) Table 1
Young's modulus	35 GPa	Han et al. (2017) Table 1
Poisson's ratio	0.26	Han et al. (2017) Table 1
Tensile strength	2.1 MPa	Han et al. (2017) Table 1
Cohesion	26 MPa	Han et al. (2017) Table 1
Friction angle	26.4°	Han et al. (2017) Table 1
Thermal conductivity	2.5 W/(m·K)	Han et al. (2017) Table 1
Specific heat	920 J/(kg·K)	Han et al. (2017) Table 1
Thermal expansion	8.0 × 10⁻⁶ /K	Han et al. (2017) Table 1
Fracture toughness KIC	1.5 MPa·m¹/²	Kant et al. (2017) Table 2
Crack length	20 µm	Kant et al. (2017) Table 2

References

• Bäuerle, D. (2011). Laser Processing and Chemistry, 4th ed. Springer. §§6–7: heat conduction under CW Gaussian illumination; Eqs (6.2.2) and (7.5.2).

• Han, Y., Fang, Y., San-Roman-Alerigi, D. P., & Batarseh, S. I. (2017). Numerical modeling of thermal-mechanical interaction process in high power electromagnetic heating of rocks. SPE-183836-MS. Table 1 (granite properties), Eqs 1–2 (Gaussian beam), Eqs 5–6 (Mohr-Coulomb and tensile failure criteria).

• Kant, M. A., Rossi, E., Madonna, C., Höser, D., & von Rohr, P. R. (2017). A theory on thermal spalling of rocks with a focus on thermal spallation drilling. J. Geophys. Res.: Solid Earth, 122(7), 5345–5362. §2.2 (thermoelasticity), Eq 5 (uniform σ_xx), Eq 6 (stress intensity factor), Eq 8 (spallability), Eq 15 (onset temperature).

• Landau, L. D., & Lifshitz, E. M. (1986). Theory of Elasticity (Course of Theoretical Physics, Vol. 7), 3rd ed. Butterworth-Heinemann. §6: thermoelastic stress-strain relations; plane-strain formulation.

• Tikhonov, A. N., & Samarskii, A. A. (1963). Equations of Mathematical Physics. Dover. Heat equation in Cartesian and cylindrical coordinates; Green's function for the half-space.

MIT — see LICENSE.