Phase change doped sand storage for isothermal thermoelectric heat delivery

Composite thermal energy storage materials

Dopant selection rationale

Five candidate dopants were selected to span the principal thermophysical mechanisms by which solid-state additives modify packed-bed thermal performance: latent heat absorption, thermal conduction enhancement, and volumetric heat capacity augmentation. The selection drew on an initial screening of twelve candidate materials [Table 2], from which five were retained on the basis of thermal figure-of-merit, chemical compatibility with silica at temperatures up to 700 K, commercial availability, and cost^[30].

The composite index in Figure 3B is evaluated as a weighted sum of six normalized criteria. Thermal conductivity received the highest weight (0.25) because radial heat delivery to the TEG surface is governed by the effective conductivity of the annular composite layer, making it the rate-limiting property in this cylindrical geometry. Volumetric heat capacity was assigned a weight of 0.20 to reflect the sensible thermal energy accumulated over the 5 h charging window. Latent heat content, volumetric energy storage density, and chemical compatibility with the silica matrix were each weighted at 0.15. These three criteria are of comparable importance but operate through distinct mechanisms: latent heat governs isothermal buffering capacity at the phase-change transition, volumetric energy density integrates sensible and latent contributions into a single delivery metric relevant to the finite charging window, and chemical compatibility imposes a stability ceiling that no performance advantage can offset. Thermal diffusivity was assigned the lowest weight (0.10); as the ratio of conductivity to volumetric heat capacity, it carries information partly captured by the two highest-weighted terms, but its inclusion retains sensitivity to fast-transient thermal response that would otherwise be suppressed in a score dominated by steady-state properties.

Each criterion was normalized to [0, 1] by min-max scaling over the twelve candidates: $$ \tilde{x} $$ = (x - x_min)\/(x_max - x_min), with chemical compatibility treated as binary. A one-at-a-time weight perturbation of ±50%, with remaining weights rescaled to preserve unity, left Solar Salt, flake graphite, copper turnings, and SiC among the top five in all twelve variants, with Kendall rank correlations against the baseline exceeding 0.85 throughout, confirming that the selection is not sensitive to the precise weight assignments within physically reasonable bounds.

Solar Salt (60 wt% NaNO₃, 40 wt% KNO₃) was retained as the sole PCM, chosen for its well-characterized latent heat (L_d ≈ 100 kJ kg^-1 and a solid-liquid transition temperature (T_m ≈ 493 K) falling within the anticipated TEG operating window^[23,24]. The isothermal buffering effect of this plateau is evident in Figure 4A. Flake graphite (k_d ≈ 150 W m^-1 K^-1) and copper turnings (k_d ≈ 400 W m^-1 K^-1) occupy the high-conductivity end of the design space, targeting aggressive suppression of bed thermal resistance^[22,31]. The superior thermal delivery of these high-conductivity formulations relative to the pure sand baseline is also reflected in Figure 4A. Iron filings were included as a high-volumetric-heat-capacity additive (ρ_dC_p,d ≈ 3.5 MJ m^-3 K^-1) to evaluate whether enhanced sensible storage density can offset modest conductivity improvements^[32]. Silicon carbide (SiC, k_d ≈ 120 Wm^-1 K^-1) probes an intermediate regime, offering simultaneous elevation of conductivity and superior chemical inertness at elevated temperatures relative to metallic or carbonaceous fillers^[21,33]. The contrast in radial heat distribution between the pure sand baseline and the Solar Salt composite at t = 5 h is illustrated in Figure 4B and C, respectively.

Figure 4. Simulated TEG interface temperature during charging (A) and radial temperature distributions at t = 5 h for pure sand (B) and 30 vol% layered Solar Salt (C). TEG: Thermoelectric generator; PCM: phase change material.

Seven additional candidate materials were evaluated during the screening phase but ultimately excluded [Table 2]. Alumina (Al₂O₃) offers excellent chemical stability but its thermal conductivity at the grain scale (k ≈ 30 W m^-1 K^-1) provides insufficient enhancement over quartz sand to justify inclusion^[6]. Magnetite (Fe₃O₄) possesses favorable volumetric heat capacity (ρC_p ≈ 3.6 MJ m^-3 K^-1) but undergoes irreversible oxidation to hematite (Fe₂O₃) above 570 K in ambient atmosphere, compromising long-term property stability^[34]. Steatite, a magnesium silicate ceramic, was excluded for its low thermal conductivity (k ≈ 3 W m^-1 K^-1) despite high volumetric heat capacity^[35]. Aluminium turnings, although highly conductive (k ≈ 237 W m^-1 K^-1), melt at 933 K and suffer progressive surface oxidation that generates insulating Al₂O₃ layers and degrades interfacial contact conductance^[36]. Boron nitride (hexagonal, k_⊥ ≈ 30 W m^-1 K^-1, k_∥ ≈ 400 W m^-1 K^-1) was excluded because of the strong anisotropy in its thermal conductivity, which renders effective medium predictions unreliable without explicit microstructural resolution^[37]. Calcium chloride hexahydrate (CaCl₂·6H₂O) melts at 302 K, far below the target operating window^[36]. Sodium hydroxide (NaOH, T_m ≈ 596 K) presents severe corrosion risk to both quartz sand and metallic containment at the temperatures of interest^[38].

Interfacial thermal resistance and heat flux management

Transient thermal-phase change simulation

The interfacial thermal resistance between the composite storage medium and the TEG substrate was quantified through a dedicated 2D axisymmetric transient simulation incorporating phase-change dynamics, implemented in COMSOL Multiphysics 6.4. The simulation workflow is summarized in Figure 5.

Figure 5. Numerical workflow for the transient thermal simulation. (A) 3D geometric model and material stacking; (B) Applied boundary conditions with a pulsed heat flux; (C) Simulated cross-sectional temperature field. TIM: Thermal interface material.

The computational domain comprised four vertically stacked layers: the composite sand bed (incorporating the 30 vol% layered Solar Salt formulation identified in the preceding section^[39]), a compliant TIM layer of thickness d_TIM = 0.1 mm, an alumina (Al₂O₃) TEG substrate, and the TEG module body [Figure 5A]. The TIM layer was assigned a thermal conductivity k_TIM = 3.5 W m^-1 K^-1, representative of commercially available silicone-based gap pads used in power electronics packaging^[40]. The interfacial thermal resistance across this layer is expressed as:

(13) $$ R_{\mathrm{TIM}}=\frac{d_{\mathrm{TIM}}}{k_{\mathrm{TIM}}}+R_{\mathrm{c}, \text { upper }}+R_{\mathrm{c}, \text { lower }} $$

where R_c,upper and R_c,lower are the contact resistances at the upper (sand bed-TIM) and lower (TIM-substrate) interfaces, respectively. For the baseline simulation, these contact resistances were set to 5 × 10^-5 m² K W^-1 each, yielding a total interfacial resistance R_TIM ≈ 1.29 × 10^-4 m² K W^-1. The corresponding interfacial heat flux delivered to the TEG substrate under steady-state conditions is:

(14) $$ q_{\mathrm{TEG}}=\frac{T_{\mathrm{bed}}-T_{\mathrm{sub}}}{R_{\mathrm{TIM}}} $$

The contact resistances R_c,upper and R_c,lower in Equation (13) were assigned constant values of 5 × 10^-5 m² K W^-1, consistent with values reported for lightly clamped silicone-based gap pads against ceramic surfaces^[41]. In practice, contact resistance decreases with increasing contact pressure as interfacial asperities deform, with experimental correlations indicating a reduction of approximately 30%-50% as pressure increases from 0.5 MPa to 2.0 MPa for soft polymer-ceramic interfaces. The fixed-resistance assumption is therefore most representative of the reference condition at P_contact = 1.0 MPa. Across the parametric sweep, this introduces a bounded conservatism whose net effect on the identified optimum is modest, given that the sensitivity of η_del to R_c is second-order relative to its sensitivity to k_TIM in this regime. Incorporating a pressure-dependent contact resistance correlation into Equation (13) is recommended as a refinement prior to prototype testing.

The mesh strategy was designed to resolve the steep thermal gradients within and adjacent to the 0.1 mm TIM gap^[41]. A five-layer boundary layer mesh was deployed on both the upper and lower surfaces of the TIM domain, with a first-layer thickness of 2 μm and a geometric growth ratio of 1.5 (Figure 5A, inset). The maximum element size was constrained to 5 μm within the TIM layer itself, transitioning to a coarser mapped quadrilateral mesh in the bulk sand bed and substrate domains. A mesh independence study confirmed that further refinement of the boundary layer beyond five sub-layers altered the computed interfacial heat flux by less than 0.4%.

The thermal boundary conditions replicated the operating scenario of the sand battery during a single charging cycle [Figure 5B]. A time-dependent pulsed heat flux was applied at the bottom boundary to represent the resistive heating input:

(15) $$ q_{\text {in }}(t)=12000 \cdot\left(1-\exp \left(-\frac{t}{3600}\right)\right) \mathrm{W} \cdot \mathrm{~m}^{-2} $$

where the exponential ramp approximates the thermal inertia of the resistive heating element during start-up, reaching 95% of the steady-state flux of 12,000 Wm^-2 within the first hour. The top surface was subjected to a convective boundary condition with h = 150 W m^-2 K^-1, representing the combined natural convection and radiative exchange at the TEG cold side. The initial temperature was set uniformly to 300 K throughout all domains. The transient solver employed BDF, and the apparent heat capacity method [Equation (8)] was retained for the Solar Salt layer to capture phase-change nonlinearity^[42].

The computed 2D temperature field at t = 2.5 h is presented in Figure 5C. The thermal field reveals a pronounced temperature discontinuity across the TIM gap, with the sand bed side reaching approximately 520 K while the substrate side remains near 480 K, confirming that the 0.1 mm interfacial layer accounts for a temperature drop of approximately 40 K. This drop is comparable in magnitude to the temperature differential exploited by the TEG module itself, indicating that the interfacial resistance constitutes a non-negligible fraction of the total thermal pathway resistance. The fraction of stored thermal energy that is effectively deliverable to the TEG interface within the 5 h charging window was quantified by the thermal delivery efficiency:

(16) $$ \eta_{\mathrm{del}}=\int_{0}^{t_{f}} q_{\mathrm{TEG}}(t) d t / \int_{0}^{t_{f}} q_{\text {in }}(t) d t $$

For the baseline TIM configuration, η_del was computed to be 0.74, indicating that 26% of the input thermal energy remains trapped upstream of the interfacial barrier at the end of the charging cycle^[43]. A parametric sweep over k_TIM from 1.0 W m^-1 K^-1 to 10.0 W m^-1 K^-1 showed that increasing the TIM conductivity to 6.0 W m^-1 K^-1 raised η_del to 0.82, beyond which further gains were marginal.

Thermo-mechanical coupling and fatigue life assessment

The cyclic thermal loading imposed during repeated charge-discharge operation induces mechanical stresses at material interfaces where CTE values are mismatched. To evaluate the structural integrity and fatigue life of the TIM joint under service conditions, a one-way thermo-mechanical coupling analysis was performed in ANSYS Workbench 2024 R1.

A three-dimensional conformal mesh comprising approximately 500,000 hexahedral-dominant elements was generated over the full cylindrical assembly [Figure 6A]. Edge fillets of 0.5 mm radius were introduced at the perimeter corners of the TIM layer to suppress the geometric stress singularity that would otherwise arise at sharp bimaterial interfaces.

Figure 6. Thermo-mechanical and fatigue life analysis of the assembly. (A) Conformal mesh with corner fillets; (B) One-way coupling topology; (C) Imported transient temperature load; (D) Thermal stress concentration and fatigue life distribution.

The coupling topology is shown in Figure 6B: the transient temperature history computed in the Transient Thermal module was mapped as an imported body temperature load onto the Static Structural module, enabling sequential resolution of the thermal and mechanical fields on a shared mesh. This sequential strategy avoids the convergence difficulties associated with fully coupled thermomechanical solvers while preserving the fidelity of the thermal driving force at each structural time step.

The mesh was locally refined at the TIM boundaries to capture the steep stress gradients associated with the coefficient of thermal expansion (CTE) mismatch between the composite sand bed (α_TES ≈ 11.5 × 10^-6 K^-1) and the alumina substrate ($$ \alpha_{A l_{2} O_{3}} $$ ≈ 7.2 × 10^-6 K^-1) [Figure 6C]. The thermally induced strain at the interface scales with the product of the CTE mismatch and the applied temperature differential:

(17) $$ \varepsilon_{\mathrm{th}}=\Delta \alpha \cdot \Delta T=\left(\alpha_{\mathrm{TES}}-\alpha_{\mathrm{Al}_{2} \mathrm{O}_{3}}\right) \cdot\left(T(t)-T_{\mathrm{ref}}\right) $$

where T_ref = 300 K is the stress-free reference temperature. For a peak interfacial temperature of 520 K, this expression yields ε_th ≈ 9.46 × 10^-4, corresponding to a thermal mismatch stress that concentrates at the perimeter edge of the TIM joint where the radial constraint is most severe.

The transient temperature field T(t) was imported as a nodal body temperature load onto the structural model [Figure 6C], and a static structural solution was obtained at each time step corresponding to the peak and valley temperatures of a representative charge-discharge cycle (300 K to 520 K to 300 K). The equivalent (von Mises) stress field exhibited pronounced concentration at the perimeter edge of the TIM layer, where the displacement magnification factor reached approximately 100 under the applied thermal loading [Figure 6D]. The maximum von Mises stress at this location was 48.3 MPa, which remains below the tensile strength of the silicone-based TIM (typically 55-65 MPa) but exceeds 74% of the yield threshold, placing the joint in the high-cycle fatigue regime^[44].

Fatigue life estimation was performed using the ANSYS Fatigue Tool with the Goodman mean-stress correction. The alternating stress amplitude σ_a and mean stress σ_m were extracted from the computed stress cycle at the critical perimeter location. The fatigue life N_f was evaluated from the modified Basquin relation:

(18) $$ \sigma_{a}=\sigma_{f}^{\prime} \cdot\left(1-\frac{\sigma_{m}}{\sigma_{\mathrm{UTS}}}\right) \cdot\left(2 N_{f}\right)^{b} $$

where $$ \sigma_{f}^{\prime} $$ is the fatigue strength coefficient, σ_UTS is the ultimate tensile strength, and b is the fatigue strength exponent of the TIM material. The computed fatigue life at the most critically stressed location was N_f ≈ 8,400 cycles. Assuming one full charge-discharge cycle per day, this corresponds to a projected service life of approximately 23 years, which exceeds the 20-year design horizon specified for the HEV infrastructure.

The 20-year horizon corresponds to N = 7,300 cycles at one cycle per day and reflects the typical asset depreciation period for community-scale off-grid infrastructure, it serves as the structural design criterion for the containment vessel. The 23-year figure (N_f ≈ 8,400 cycles) is a computed output of the Goodman fatigue model at the baseline interface condition, not an independently specified target, and its role is to confirm that a sufficient margin exists before the interface is driven toward higher conductivity in the optimization study of Section 3.2.

The thermal ratcheting deformation visible at the outer perimeter in Figure 6D remained below 0.12 mm per cycle, confirming that progressive plastic accumulation does not compromise the geometric stability of the TIM joint within the anticipated service window^[45].

Thermal stress and material degradation

Thermomechanical simulation of the multi-layer storage assembly

To evaluate the global thermomechanical response of the storage vessel under sustained cyclic operation, a 2D axisymmetric simulation was constructed in COMSOL Multiphysics 6.4 using a one-way coupling between the Heat Transfer in Solids (ht) and Solid Mechanics (solid) modules. The computational domain comprised four concentric annular regions: the central heating element, the composite storage layer (inner pure-sand core and outer 30 vol% layered Solar Salt-sand composite), a ceramic insulation shell (alumina, 5 mm), and an outer metallic containment shell (stainless steel, 3 mm), as depicted in Figure 7A. All inter-domain boundaries were constructed with shared nodes to enforce displacement continuity^[46,47].

Figure 7. Computational framework for thermo-mechanical simulation. (A) 2D axisymmetric geometry and material distribution; (B) Finite element mesh with interfacial refinement; (C) One-way multiphysics coupling workflow; (D) Cyclic thermal loading profile.

Each domain was assigned a linear elastic constitutive model with isotropic thermal expansion referenced to T_ref = 300 K. The CTE values span nearly two orders of magnitude: α_sand ≈ 0.5 × 10^-6 K^-1, α_salt ≈ 40 × 10^-6 K^-1 (liquid phase), α_ceramic ≈ 8.0 × 10^-6 K^-1, and α_steel ≈ 17.3 × 10^-6 K^-1. This mismatch, identified as the principal driver of interfacial degradation in multi-material TES structures^[48,49], governs the radial thermal stress at any interface between adjacent layers i and j:

(19) $$ \sigma_{r}^{i j}(t)=\frac{E_{i} E_{j}}{E_{i}\left(1-v_{j}\right)+E_{j}\left(1-v_{i}\right)}\left(\alpha_{i}-\alpha_{j}\right)\left[T\left(r_{i j}, t\right)-T_{r e f}\right] $$

where i and j denote the indices of the two adjacent material layers; E_i(E_j), ν_i(ν_j), and α_i(α_j) denote the Young’s modulus, Poisson’s ratio, and coefficient of thermal expansion of layer i (layer j), respectively; and r_ij is the radial coordinate of their shared boundary. The finite-element mesh incorporated five-layer boundary layer elements at each material interface with a first-layer thickness of 0.05 mm [Figure 7B], and the total element count was maintained between 60,000 and 80,000, which a convergence study confirmed to be adequate for the 2D axisymmetric domain^[50,51]. This boundary-layer resolution at each bi-material junction is necessary to capture the steep stress gradients arising from CTE mismatches that span nearly two orders of magnitude across adjacent layers. The temperature field T(r,t) was obtained from the transient heat transfer solution and mapped onto the structural domain at each time step through the built-in variable ht.T, completing the one-way coupling illustrated in Figure 7C. The cyclic thermal loading followed the boundary condition profile shown in Figure 7D: each cycle comprised a 5 h charging phase during which the pulsed heat flux [Equation 15)] was applied at the inner boundary, followed by a 5 h cooling phase with zero heat input and natural convective dissipation (h = 5 - 10 W m^-2 K^-1, T_amb = 300 K) at the outer shell surface.

Composite thermal energy storage materials

Simulations were conducted over a 5 h charging window. The TEG-side boundary temperature T_TEG (t) served as the primary performance metric. The undoped quartz sand baseline rose monotonically from 300 K to approximately 420 K over this period, establishing the reference trajectory against which all composite formulations were evaluated.

The Solar Salt configurations produced a qualitatively distinct thermal response not observed in any other dopant group [Figure 8A]. As T_TEG (t) approached the eutectic transition temperature of approximately 493 K, latent heat absorption during the solid-liquid phase change induced a well-defined isothermal plateau, visible as a pronounced inflection in the T_TEG (t) curve. The 30 vol% layered configuration reached this plateau earliest, at approximately t = 3 h, while the 10 vol% uniform blends did not stabilize until near t = 5 h. This isothermal window was entirely absent in every other dopant group and represents a qualitative, not merely quantitative, advantage for thermoelectric coupling: a stable hot-side temperature directly fixes the Seebeck voltage and suppresses efficiency oscillations in the TEG modules. Purely sensible-heat dopants cannot reproduce this buffering behavior regardless of loading fraction.

Figure 8. Simulated thermal performance of composite TES beds with five candidate dopants. (A-E) TEG-side boundary temperature as a function of charging time for uniform and layered configurations at 10%, 20%, and 30% dopant volume fractions; (F) Three-dimensional temperature field of the pure sand baseline at t = 5 h; (G) Three-dimensional temperature field of the 30 vol% layered configuration at t = 5 h. Comparison of panels (F) and (G) illustrates the redistribution of stored thermal energy toward the TEG boundary enabled by the phase-change front in the layered Solar Salt configuration. TES: Thermal energy storage; TEG: thermoelectric generator.

Graphite delivered the most aggressive early-stage heating among the five candidates [Figure 8B]. The 30 vol% layered variant reached approximately 660 K at t = 5 h, while the 30 vol% uniform blends reached around 650 K, both substantially above the pure sand baseline. The initial heating rate for the 30% uniform configuration approached 110 K h^-1 during the first hour, a direct consequence of the elevated effective thermal diffusivity. The marginal gains between 20% and 30% loading were small relative to the step from 10% to 20%, consistent with a percolation-mediated conductivity transition in which a continuous high-conductivity network spanning the bed provides diminishing returns in conductivity with further addition while progressively reducing volumetric heat capacity.

Iron filings and SiC occupied an intermediate thermal performance regime [Figure 8C and D]. Iron filings exhibited steady, near-linear temperature evolution governed by their high volumetric heat capacity, reaching approximately 620 K under 30 vol% layered doping at t = 5 h. SiC achieved a marginally higher terminal temperature of approximately 645 K under the same loading conditions, attributable to its higher intrinsic conductivity relative to iron. Neither material displayed plateau behavior; temperature profiles remained monotonically ascending throughout the 5 h window. The absence of an isothermal window means that the TEG hot-side temperature drifts continuously during operation, introducing time-varying mismatch between the storage subsystem and downstream power conditioning electronics.

Copper turnings produced the highest absolute temperatures across all configurations tested [Figure 8E], reaching approximately 660 K at t = 5 h under 30 vol% layered doping. The early-phase rise from 300 K to around 550 K within the first 2 h is consistent with copper’s outstanding thermal conductivity of approximately 400 W m^-1 K^-1. Practical deployment of copper, however, raises non-trivial engineering concerns. The high material density (8,960 kg m^-3) introduces a substantial gravimetric penalty, and progressive surface oxidation at sustained temperatures above 500 K may degrade interfacial thermal contact conductance over repeated charge-discharge cycles. These considerations collectively preclude copper turnings from further consideration despite its conductivity advantage.

Across all five dopant systems, layered configurations consistently outperformed their uniform counterparts at equivalent dopant fraction, with the performance disparity most evident at lower loadings. The dominant thermal resistance in the radial direction is concentrated in the outer annular shell, where the temperature gradient must drive outward heat flux to the TEG surface. By concentrating the high-conductivity dopant in this rate-limiting zone, the layered architecture reduces the thermal bottleneck while preserving thermal inertia of the inner core. Uniform mixing, by contrast, distributes the conductivity enhancement across the interior where its contribution to TEG-side heat delivery is attenuated by the 1/r dependence of radial heat flux in cylindrical geometry.

To quantify the extent of this architectural advantage, the total radial thermal resistance is decomposed into its inner-core and outer-annulus contributions using the cylindrical-shell formula R = ln(r_out/r_in)/(2πkL). For the 30 vol% Solar Salt formulation, the volumetric mixing rule [Equation (3)] gives k_eff = ϕ_sk_s + ϕ_dk_d = 0.70 × 1.4 + 0.30 × 0.55 = 1.145 W m^-1 K^-1, which falls below the pure-sand value of 1.4 W m^-1 K^-1. In the uniform configuration, this reduced conductivity governs the entire radial domain; in the layered configuration, the inner core retains k_s = 1.4 W m^-1 K^-1 while the composite is confined to the outer annulus. Using the vessel geometry (R_i = 0.025 m, R_c = 0.20 m, R_o = 0.30 m), the total radial resistance per unit length for the uniform case is R_unif = ln(R_o/R_i)/(2πk_eff) = 0.345 m K W^-1, while the layered case yields R_layer = ln(R_c/R_i)/(2πk_s) + ln(R_o/R_c)/(2πk_eff) = 0.236 + 0.056 = 0.292 m K W^-1, a reduction of approximately 15%. This outcome reflects a property specific to Solar Salt: because k_d < k_s, distributing the PCM uniformly penalizes conductive transport in the inner domain without any compensating latent-heat benefit there. Excluding Solar Salt from the inner core preserves the higher-conductance sand pathway over the greater part of the radial domain, while the outer-annulus resistance remains identical in both architectures. Beyond the conductive pathway, approximately 44% of the total Solar Salt mass in the uniform configuration resides in the inner domain (r < R_c), where phase-change absorption at 493 K proceeds remote from the TEG surface and cannot directly regulate T (R_o,t). In the layered configuration, all PCM volume is confined to the outer annulus, ensuring that the complete latent heat reservoir acts on the zone governing the TEG hot-side temperature, consistent with the formation of a near-isothermal band between 485 K and 500 K within the outer annulus. The 15% lower total radial resistance and the full spatial concentration of latent heat capacity at the TEG boundary together constitute the quantitative basis for the performance advantage of the layered architecture over the uniform configuration at equivalent dopant loading.

The COMSOL thermal field maps in Figure 8F and G provides spatial evidence for the mechanisms inferred from the T_TEG (t) curves, comparing the pure sand baseline against the 30 vol% layered Solar Salt configuration at hourly intervals from t = 1 h through t = 5 h. In the pure sand baseline [Figure 8F], the temperature field at t = 1 h shows a steep radial gradient concentrated in the inner third of the domain, with the central heating element above 550 K while the outer annulus remains below 340 K. The low thermal diffusivity of quartz sand confines the thermal penetration front to a narrow region around the heat source. By t = 5 h, the TEG-side temperature reaches only approximately 420 K, confirming that internal conduction resistance, rather than external convective coupling, limits the rate of useful heat delivery.

The 30 vol% layered Solar Salt configuration [Figure 8G] shows a markedly different spatial evolution. At t = 1 h, the thermal front has already penetrated further into the outer annulus relative to the pure sand case, because the composite shell provides a lower-resistance pathway for radial heat transport. By t = 2 h, the temperature within the outer composite annulus has risen above 450 K with a notably shallower radial gradient than in the undoped baseline. Between t = 3 h and t = 4 h, the outer annulus traverses the Solar Salt melting range near 493 K, and the latent heat sink absorbs the incoming thermal flux without a corresponding temperature rise. The field map at t = 3 h reveals a nearly isothermal band within the outer shell, with temperature contours compressed into a narrow interval near 485 K to 500 K. By t = 5 h, phase change is largely complete, and the TEG interface receives heat through a composite layer whose temperature varies by less than 10 K across its radial extent, compared with more than 80 K variation across the same region in the pure sand baseline. The interior hot-spot temperature in the pure sand case exceeds 700 K at t = 5 h while the TEG surface reaches only 420 K, a 280 K internal mismatch representing stored energy that cannot be extracted within the charging window. In the layered Solar Salt configuration, the maximum interior temperature is lower (approximately 640 K) because enhanced radial transport draws heat outward more effectively, and the TEG surface temperature stabilizes near 495 K, reducing the internal mismatch to roughly 145 K.

The comparative analysis identifies 30 vol% layered Solar Salt as the preferred composite formulation for this platform. While copper turnings achieved the highest absolute T_TEG (t), the operational objective of a sand battery coupled to thermoelectric generators is not to maximize peak temperature but to deliver a temporally stable thermal driving force over the full charging cycle. For a monotonically rising T_TEG (t), as exhibited by copper, graphite, SiC, and iron filings, the continuous drift of the operating point generates a time-varying mismatch between the thermally optimal load resistance and any fixed downstream power conditioning circuit, reducing cycle-averaged electrical output below the theoretical maximum. The 493 K isothermal plateau furnished by Solar Salt directly resolves this mismatch: the latent heat of fusion absorbs input variability and re-emits stored energy at a near-constant temperature. The spatial field maps in Figure 8G confirm that this temporal plateau originates from a spatially coherent phase-change front propagating through the outer annulus, not from a cancellation of opposing gradients.

The choice of the layered over the uniform architecture at 30 vol% is further supported on thermomechanical grounds: confining the PCM to the outer annulus ensures that the inner core retains the mechanical integrity and dimensional stability of pure quartz sand, reducing the risk of thermally induced cracking near the heating element during repeated thermal cycling. Among PCM candidates that could provide a latent heat plateau within the target temperature window, Solar Salt is distinguished by its established industrial track record in CSP, its cost advantage relative to metallic or advanced ceramic fillers, and its demonstrated chemical compatibility with silica-based matrices at the temperatures of interest.

Interfacial thermal performance and thermo-mechanical reliability

The preceding sections established that the 30 vol% layered Solar Salt composite maximizes TEG-side thermal delivery among the candidate formulations investigated. However, the practical realization of this thermal advantage depends critically on the efficiency with which heat traverses the interface between the granular storage medium and the solid TEG substrate. This section presents a parametric analysis of the coupled effects of TIM thermal conductivity (k_TIM = 1.0~10.0 W m^-1 K^-1) and mechanical contact pressure (P_contact = 0.5~3.5 MPa) on thermal delivery efficiency, transient heat flux, interfacial temperature drop, thermomechanical stress distribution, and cyclic degradation behavior, with the objective of identifying the design envelope that simultaneously satisfies the thermal performance and structural durability requirements of the HEV sand battery module.

The seven panels of Figure 9 collectively demonstrate that thermal delivery efficiency and mechanical durability are governed by the interaction between TIM conductivity and contact pressure, and the following discussion traces the physical mechanisms and quantitative bounds of this coupling. The coupled influence of TIM thermal conductivity and mechanical contact pressure on the thermal delivery efficiency of the sand battery assembly is mapped in Figure 9A. The contour field reveals that η_del increases monotonically with both k_TIM and P_contact, but the sensitivity is asymmetric: below k_TIM ≈ 3.5 W m^-1 K^-1, efficiency gains are dominated by conductivity improvements, whereas above this threshold the marginal return from further conductivity enhancement diminishes and contact pressure becomes the controlling parameter. The baseline configuration (k_TIM = 3.5 W m^-1 K^-1, P_contact = 1.0 MPa) yields η_del = 0.74, indicating that 26% of the input thermal energy remains inaccessible to the TEG interface at the end of the 5 h charging window. Increasing k_TIM to 6.0 W m^-1 K^-1 at P_contact = 2.0 MPa raises η_del to 0.82, while the most aggressive configuration (k_TIM = 10.0 W m^-1 K^-1, P_contact = 3.5 MPa) achieves 0.89. The Pareto front extracted from the η_del - N_f trade-off space traces a trajectory from the lower-left to the upper-right of the design map, confirming that no single parameter can be optimized in isolation without compromising either thermal performance or structural durability. The inset marginal profile at P_contact = 2.0 MPa shows that η_del saturates above k_TIM ≈ 6.0 W m^-1 K^-1, suggesting limited practical benefit from ultra-high-conductivity interface materials. This asymmetry is also reflected in the contour spacing of Figure 9A, where the iso-efficiency lines are closely spaced along the conductivity axis below 3.5 W m^-1 K^-1 and widen appreciably above it, while their spacing along the pressure axis remains comparatively uniform. The weak pressure dependence in this regime is consistent with the bounded influence of the constant contact-resistance assumption adopted in the model, and it identifies interface conductivity, rather than clamping pressure, as the parameter most worth optimizing once a moderate contact load has been secured.

Figure 9. Parametric analysis of interfacial thermal performance and thermo-mechanical reliability. (A) Coupled effects of thermal conductivity and contact pressure on thermal delivery efficiency; (B) TEG-side temperature and instantaneous conversion efficiency under different contact pressures. Solid curves denote T_TEG on the left axis, and dashed curves denote η_TEG on the right axis; the blue and yellow dashed curves correspond to P_contact = 0.5 and 3.5 MPa, respectively; (C) Heat flux and cumulative delivered thermal energy under different TIM conductivities. Solid curves denote q_TEG on the left axis, and dashed curves denote cumulative delivered thermal energy on the right axis; the blue and yellow dashed curves correspond to k_TIM = 1.0 W m^-1 K^-1 and 10.0 W m^-1 K^-1, respectively; (D) Interfacial temperature drops as a function of TIM conductivity and contact pressure; (E) Radial stress distribution indicating edge stress concentration; (F and G) Cyclic degradation trajectories and comprehensive comparison balancing thermal efficiency and fatigue life. Panel (G) identifies the design point at k_TIM = 6.0 W m^-1 K^-1 and P_contact = 2.0 MPa as providing the optimal balance between thermal delivery efficiency and projected service life. TEG: Thermoelectric generator; TIM: thermal interface material.

The feasible design region for the structural containment is bounded by η_del ≥ 0.74 and N_f ≥ 7,300 cycles (corresponding to the 20-year service horizon at one cycle per day). For the TIM joint, the fatigue criterion is applied differently: as a compliant gap pad at the accessible outer vessel surface, the TIM is replaceable during a planned maintenance shutdown, so its fatigue life governs the length of the servicing interval rather than the structural design life of the system. Superimposed isothermal contours of ΔT_TIM provide a complementary perspective: configurations falling below the ΔT_TIM = 20 K isoline maintain interfacial temperature drops comparable to or smaller than the TEG operating differential, whereas those exceeding 40 K suffer parasitic thermal losses that substantially erode conversion efficiency.

The transient TEG-side temperature histories at four contact pressures [Figure 9B] reveal that elevated contact pressure accelerates the thermal response, with the 3.5 MPa case reaching the Solar Salt phase-change plateau near 493 K approximately 0.8 h earlier than the 0.5 MPa case. The instantaneous TEG conversion efficiency, plotted on the secondary axis, increases with temperature differential and stabilizes during the isothermal plateau, confirming that the PCM buffering effect identified in the composite optimization study translates through the interfacial barrier to the TEG operating point. Figure 9C quantifies the corresponding heat flux delivered to the TEG substrate for four k_TIM values at fixed P_contact = 2.0 MPa. The k_TIM = 10.0 W m^-1 K^-1 configuration sustains a terminal flux approximately 2.4 times that of the k_TIM = 1.0 W m^-1 K^-1 case, and cumulative energy delivery over the charging window scales accordingly. The parametric sweep of ΔT_TIM against k_TIM [Figure 9D] confirms that only configurations with k_TIM above approximately 5 W m^-1 K^-1 at P_contact ≥ 2.0 MPa satisfy the 20 K target threshold across the full pressure range investigated.

The radial distribution of von Mises stress across the TIM joint [Figure 9E] exhibits a pronounced concentration at the outer perimeter (r/R_o > 0.85), where the CTE mismatch between the composite sand bed and the alumina substrate is mechanically constrained. The peak stress at P_contact = 0.5 MPa reaches 62 MPa, exceeding the TIM ultimate tensile strength of 55 MPa, indicating that insufficient clamping pressure permits localized edge debonding. At P_contact = 2.0 MPa, the peak stress reduces to 48.3 MPa (safety factor of 1.14), consistent with the fatigue life estimate of 8,400 cycles reported in the thermo-mechanical analysis. Higher contact pressures further suppress edge stresses but introduce compressive creep concerns not captured by the elastic model.

The long-term degradation trajectories of η_del under cyclic thermal loading [Figure 9F] show that all four configurations exhibit an initial exponential decay followed by stabilization beyond approximately 5,000 cycles. The baseline configuration (k_TIM = 3.5 W m^-1 K^-1, P_contact =1.0 MPa) reaches 80% retention of its initial η_del at approximately 8,400 cycles, consistent with the fatigue life prediction. The low-conductivity, low-pressure case (k_TIM = 1.0 W m^-1 K^-1, P_contact = 0.5 MPa) degrades below the 80% retention threshold before 5,000 cycles, rendering it unsuitable for the 20-year design window. The summary comparison across five representative configurations [Figure 9G] confirms that the optimized design point (k_TIM = 6.0 W m^-1 K^-1, P_contact = 2.0 MPa) provides the most favorable balance between thermal delivery efficiency and fatigue life, achieving η_del = 0.82 with N_f ≈ 6,200 cycles. This 6,200-cycle life is lower than the 8,400 cycles obtained at the baseline interface (k_TIM = 3.5 W m^-1 K^-1), and the difference is a direct consequence of the higher interface conductivity. Raising k_TIM increases the heat flux delivered across the joint and therefore the absolute interfacial temperature sustained during each charge-discharge cycle, which enlarges the cyclic thermal-stress amplitude at the TIM perimeter and correspondingly reduces the fatigue margin. The gain in delivery efficiency and the reduction in fatigue life are thus two outcomes of the same increase in interfacial heat transfer, representing an inherent efficiency-durability trade-off. Although the k_TIM = 10.0 W m^-1 K^-1 configuration attains the highest η_del of 0.89, its fatigue life of 3,200 cycles falls below the 20-year threshold, excluding it from the feasible design region.

Taken together, the parametric results establish that interfacial thermal management constitutes a first-order design constraint for compact sand-based TES modules, comparable in importance to the bulk composite formulation. The thermal delivery efficiency of the assembly is governed not by a single interfacial property but by the coupled interaction between TIM conductivity and mechanical contact pressure, with diminishing returns in k_TIM above approximately 6 W m^-1 K^-1 shifting the optimization burden to contact pressure control. The identification of a bounded feasible region in the k_TIM - P_contact space, constrained simultaneously by η_del ≥ 0.74 and N_f ≥ 7,300 cycles, provides a quantitative basis for material selection and mechanical design of the interface joint. The optimized configuration (k_TIM = 6.0 W m^-1 K^-1, P_contact = 2.0 MPa) delivers 82% of the stored thermal energy to the TEG interface with a projected TIM fatigue life of approximately 6,200 cycles (~17 years). This value is below the 7,300-cycle structural threshold but reflects the distinct serviceability of the TIM layer: unlike the welded containment vessel, the gap pad is accessible at the outer surface and can be replaced in a single planned maintenance event, after which the assembly returns to its original thermal-mechanical performance state. The structural containment independently satisfies N_f ≥ 7,300 cycles at a safety factor of 2.18 (Section 3.3), so the 20-year system design horizon is maintained through one scheduled TIM replacement within the operating window. This configuration therefore represents the preferred operating point for subsequent prototype development within the HEV platform.

The thermal delivery efficiency of 0.82 achieved at the optimized interface condition can be compared with the experimentally measured performance of analogous packed-bed TES configurations. Trevisan et al. (2022) reported a maximum thermal efficiency of approximately 72% in a first-of-kind radial-flow packed-bed TES prototype tested at temperatures up to 700 °C, observing that non-uniform bed porosity adversely affected the thermal performance in the outer region of the bed^[54]. The additional 10 percentage points recovered in the present study arise from the systematic co-optimization of TIM thermal conductivity and mechanical contact pressure, a coupled interfacial design dimension that has not been resolved in conventional packed-bed configurations and that the parametric mapping in Figure 9A identifies as the binding constraint above k_TIM ≈ 6.0 W m^-1 K^-1.

Cyclic stress response and material-pair selection

Having established the thermomechanical modeling framework and the candidate material matrix in the preceding section, attention now turns to the quantitative response of the assembly under sustained cyclic operation and to the identification of a material pairing that reconciles mechanical reliability with thermal-shock resistance and capital cost. Figure 10 presents six complementary analyses - corresponding to the six panels (A) through (F) - applied across the nine insulation-shell material combinations to support the final material-pair selection.

Figure 10. Cyclic stress response and material-pair selection. (A) Peak von Mises stress over 500 cycles; (B) radial stress and temperature profiles for Al₂O₃ + 304SS; (C) CTE design map with iso-lifetime contours; (D) Goodman-corrected fatigue margins; (E) multidisciplinary composite ranking across nine material combinations; (F) 2D stress field for the optimal ZrO₂ + Inconel pairing. Panel (E) consolidates the multidisciplinary ranking across nine material combinations, confirming ZrO₂ paired with Inconel 718 as the preferred selection on the basis of fatigue safety factor, CTE compatibility, and thermal shock resistance. CTE: Coefficient of thermal expansion.

Figure 10A tracks the peak von Mises stress at the insulation-shell interface over 500 charge-discharge cycles for the three retained pairings. All three curves exhibit a two-stage behavior: an initial transient softening during the first 20-40 cycles, followed by a quasi-stationary plateau. Al₂O₃ + 304SS settles near 240 MPa, driven by the 9.2 × 10^-6 K^-1 CTE mismatch, while Mullite+carbon steel (C-steel) and ZrO₂ + Inconel stabilize at 165 MPa and 120 MPa, respectively. The shaded envelopes around each curve denote the ±3% cycle-to-cycle variation in peak stress recorded over the stabilized portion of the 500-cycle simulation, and therefore quantify the numerical scatter of the ratcheting solution at the plateau. The radial stress profile for the Al₂O₃ + 304SS baseline [Figure 10B] shows that σ_θ changes sign across the PCM-ceramic interface and reaches 260 MPa at the inner shell surface, where ceramic constraint and a 220 K temperature drop act in concert. The von Mises envelope peaks within the ceramic layer rather than at either interface, indicating that bulk ceramic failure, not interfacial debonding, is the limiting mechanism for this pair.

The CTE design map [Figure 10C] isolates the role of thermal expansion mismatch. Iso-lifetime contours at N = 3,650, 7,300, and 14,600 cycles delineate a narrow feasible corridor centered on α_ceramic ≈ 10 - 12 × 10^-6 K^-1 and α_shell ≈ 11 - 14 × 10^-6 K^-1). Al₂O₃ + 304SS lies outside this corridor, whereas Mullite + C-steel and ZrO₂ + Inconel fall within it, with the latter closest to the iso-stress minimum. The map assumes the reference geometry and excludes creep relaxation effects that would marginally expand the feasible region at elevated temperatures. Figure 10D quantifies the fatigue margin via the Goodman-corrected stress amplitude. The operating S_a for Al₂O₃ + 304SS (120 MPa) approaches the endurance limit of 230 MPa, yielding a safety factor of 0.89 that is unacceptable for a 7,300-cycle horizon. Mullite + C-steel returns 1.34, while ZrO₂ + Inconel attains 2.18.

The composite ranking in Figure 10E integrates fatigue life, peak stress, capital cost, and ceramic thermal-shock resistance under the weighting (0.40, 0.30, 0.15, 0.15). ZrO₂ + Inconel emerges as the preferred pairing despite its highest relative cost, owing to its low CTE mismatch, high ceramic fracture toughness (K_Ic = 8.0 MPa·m^1/2), and the favorable cyclic-hardening exponent of Inconel 718. The 2D stress field for this optimal pairing [Figure 10F] confirms that the peak stress of 142 MPa localises at the top and bottom corners of the metallic shell near r₄, where axial constraint intensifies the radial expansion mismatch, while no hot-spot is predicted at the ceramic-metal interface itself.

Taken together, these results establish that CTE compatibility, rather than absolute strength, governs the long-term viability of the multi-layer storage assembly under the imposed thermal cycle. The conventional Al₂O₃ + 304SS pairing, while adequate for short-duration prototypes, fails to meet the 20-year design horizon once realistic mean-stress corrections are applied. ZrO₂ + Inconel satisfies all four selection criteria with a fatigue safety factor above 2, and its adoption shifts the governing failure mode from material fatigue to geometric stress concentration at the shell end-caps, a consideration that should guide the subsequent detailed design of closure welds and edge fillets.

The Goodman-corrected safety factor of 2.18 reported for the ZrO₂-Inconel 718 configuration can be placed in perspective by reference to thermomechanical analyses of conventional TES containment published in the recent literature. Mishra et al. (2025) performed a validated coupled thermal-stress simulation of a commercial-scale molten salt storage tank constructed from standard 347H stainless steel and found that a single charge-discharge cycle was sufficient to drive the tank floor into the elastic-plastic regime, producing a peak von Mises stress of 107 MPa at the floor-wall junction with measurable plastic strain accumulation after just one thermal event; low-cycle fatigue and its interaction with creep were explicitly identified as unresolved concerns for long-duration service^[55]. A parallel set of observations was reported by Xu et al. (2025), who showed in a molten salt packed-bed TES model that operational conditions favouring improved thermal discharge performance simultaneously push tank wall stresses toward structural failure thresholds, illustrating the inherent conflict between thermal and mechanical performance that arises when containment materials are selected without explicit consideration of CTE compatibility^[52]. The present results suggest that this conflict can be substantially mitigated through material pairing rather than geometric reinforcement. The CTE mismatch for ZrO₂ and Inconel 718 (Δα ≈ 2.5 × 10^-6 K^-1) is approximately 3.7 times smaller than that of the Al₂O₃ + 304SS combination evaluated and rejected in this study, and this difference in mismatch propagates directly into the threefold improvement in Goodman safety factor (2.18 vs. 0.89). The quasi-stationary peak von Mises stress of 120 MPa in the preferred pairing remains comfortably within Inconel 718’s high-cycle endurance limit over the full 7,300-cycle design life, a condition that the existing literature on both steel TES tanks and rock-bed packed systems has thus far been unable to satisfy without significantly restricting operational temperature range^[29,52,55].

Table 5 consolidates the key design parameters, material selections, and quantitative performance outcomes across the three parallel investigations, providing a concise reference for practitioners adapting this framework to other household-scale TES configurations.

Several limitations should be noted. The effective medium treatment homogenizes grain-scale heterogeneities and does not capture discrete particle rearrangement or local contact network evolution, phenomena that discrete element analyses have shown to contribute additional hoop stress on containment walls. The apparent heat capacity formulation smooths the phase-change front over a 10 K interval, which may underestimate the sharpness of the thermal plateau observed in physical systems. All thermomechanical results are derived from linear elastic constitutive models; creep relaxation at sustained elevated temperatures and progressive oxidation of metallic interfaces remain unresolved. The present work was conducted entirely through numerical simulation because the target operating conditions, involving temperatures up to 700 K sustained over thousands of thermal cycles within a sealed multi-material assembly, preclude meaningful characterization through bench-scale experiments without purpose-built instrumentation for in-situ interfacial temperature and stress measurement. The simulation-first approach enables rapid parametric exploration of a design space that would be prohibitively expensive to survey experimentally, while providing quantitative targets for subsequent prototype validation.