Sand Battery Under a Slab: Storage, Losses and Floor Temperature

The axisymmetric conduction equation is solved in the sand, concrete, insulation and surrounding soil:

Cell-centred finite volumes use annular volumes and shared-face conductances formed from the two half-cell resistances. Every interface aligns with a grid face. Backward Euler advances the coupled system through a 365-day year. The default main grid has 306 cells: 17 radial by 18 vertical, including two concrete layers and one layer across each thin insulation region.

The exterior soil surface follows a single annual harmonic. Initial soil and far radial/bottom boundaries use the Kusuda–Achenbach damped, phase-lagged solution

This follows Kusuda and Achenbach (1965), NBS Report 8972, PDF pp. 28–29, equations (2)–(5). The phase lag with depth is retained. The upper concrete boundary connects to the fixed room through the finish resistance and combined radiation/convection coefficient; the reported floor temperature is the room-facing finish surface.

Linearity allows an unheated baseline and the heater-induced temperature rise to be advanced separately with the same operator. Their sum gives actual temperatures. PV-attributable stored heat and boundary flows use the rise alone. Net base and edge fluxes are evaluated at the outer insulation faces, with the exposed insulation annulus counted as an edge loss. Implicit endpoint fluxes make heat input equal room transfer plus base transfer plus edge transfer plus the change in assembly energy. Soil outside the insulation is outside this accounting boundary.

January means span 1 January to 1 February; stored heat is read on 15 January. October–March spans 1 October to 1 April. The main time step is adjusted slightly, if necessary, to divide the year exactly. The two three-option sweeps repeat the same two-dimensional model on a coarser grid. Independent Q1 Galerkin finite-element / Crank–Nicolson solutions in SciPy verify two annual cases and a near-zero heat-flow regression; an NBS tabulated ground temperature checks the seasonal boundary.

• Equivalent circular footprint with axisymmetry; no wall, footing or service penetration bridges.
• Constant properties in each material and perfect thermal contact. Insulation volumetric heat capacity is an illustrative 0.06 MJ/(m³ K).
• Fixed room temperature and a prescribed effective surface coefficient; no air-property correlation is used.
• Sand and insulation initially share the entered initial temperature; concrete starts at room temperature. Surrounding soil starts at the June seasonal profile.
• Constant daily-average surplus during the selected summer window; the heater is distributed uniformly across the footprint in the 60–80% depth band. No active winter extraction or charge-temperature control.
• The heat-capacity result is $m c_p\Delta T$, independent of losses and of whether the charge actually reaches that swing.

Use for low-temperature conduction studies, with sand normally below 80 °C. Verify insulation, floor finishes, concrete, drainage and electrical controls for the proposed temperatures and loads. Element sheath temperatures and local hot spots are unresolved. Floor results are area means; the reported peak floor value is the maximum of that mean over the year.

Moisture, freezing/thawing and groundwater can change both storage and ground losses substantially. Fit the seasonal inputs to ground surface, including snow and shading effects, rather than directly to air temperature. The defaults are illustrative, not Wisconsin climate measurements. A circle has less perimeter than many real slabs, so elongated or irregular footprints need a three-dimensional or plan-specific model.

This is the first June-to-May cycle, not a repeatedly charged steady annual state. The room is a prescribed reservoir: positive incremental PV heat can mean avoided room-to-ground loss, not a floor warmer than the room. Annual room transfer includes unwanted summer release and must not be treated as annual useful heating or a commercial battery efficiency.

To keep the default evaluation near 3 s, the main solve uses 306 cells and daily steps; construction sweeps use 132 cells and the selected time step. Coarser sweep energies can differ by several percent: transfer an option to the main inputs and refine it. The daily default overpredicts January residual PV heat by about 5% against the independent fine FEM cases; floor means differ by less than 0.03 °C in those cases. Mesh and time-step checks quantify retained-energy convergence. Separate January heat-flow checks enforce absolute changes below 5 W on mesh refinement and 1 W on time-step refinement in the near-zero regression. The verification vector includes a fixed 1 MW normalization so the engine's relative convergence test imposes those exact watt limits; this is not a physical heat source. Short charge fluctuations and subdaily slab dynamics need a smaller time step; two concrete cells and one cell through each insulation region do not resolve rapid within-layer transients. Refine the outer soil extent as well as the grid when long-range diffusion matters.

The temperature picture uses 65 × 49 samples, PCHIP interpolation in radius and depth, and linear interpolation between approximately ten-day saved fields. This smooths the picture without adding solution accuracy. The view shows the near-ground region, while the solver extends to the entered far boundary. Installation design, structural checks and electrical protection are outside scope.

The heat-direction verdict is neutral when the magnitude of January net heat is inside the reported dead band. At the default or finer resolution, the band is ±0.35 W/m² times the slab area (±35 W for 100 m²). It is multiplied by $\max[1,(7/n)^2]$ for a coarser sand/soil mesh and by $\max[1,\Delta t/(24\,h)]$ for a longer time step. This screening allowance covers the demonstrated 22.66 W discrepancy in the 100 m² near-zero regression against an 11,058-node independent FEM. It is not a rigorous error bound for every input combination and excludes material and boundary-condition uncertainty. Small refinement changes alone do not eliminate the bias from fixed thin-layer counts and the slab-edge corner. The signed watt result remains available when direction is unresolved; schematic arrows are suppressed inside the same band.

The supported main mesh is limited to 4–9 cells through the sand depth (at most 462 total cells), with 6–48 h time steps. Unsupported requests are rejected before allocating or marching the dense system. The largest supported mesh at a 6 h step completed in 12.3 s in the build verification, below the 30 s engine timeout; defaults remain near 3 s. Higher-resolution continuum work requires an external FEM or a different sparse solver. This ceiling limits how closely the thin layers and edge can be resolved; it does not change the conduction equation.