Numerical method
Both the 1D and the 2D models are solved with the same strategy: an implicit control-volume discretisation of the heat equation (Patankar 1980), advanced in time until the daily solution becomes periodic. This page documents the discretisation, the solution algorithm, the convergence criteria, the differences with the reference paper (Barrios et al. 2016), and the validation record.
Discretisation
The domain is divided into control volumes: \(N_x = 200\) through the thickness in 1D (config.Nx), and an \(n_x \times n_y = 80 \times 160\) grid (width × thickness) in 2D (config2d). Integrating the heat equation over each volume with a fully implicit (backward-Euler) time scheme gives, per node, an algebraic equation of the form
\[ a_P\, T_P \;=\; a_E\, T_E + a_W\, T_W + a_N\, T_N + a_S\, T_S + b, \qquad a_P = \frac{\rho c\, \Delta x \Delta y}{\Delta t} + \sum a_{nb} \;(+\,h\,\Delta s), \tag{1}\]
where the neighbour coefficients \(a_{nb}\) are the face conductances (only \(E\)/\(W\) exist in 1D), \(b\) carries the previous-time-step term and the boundary sources, and the film coefficients \(h_o\), \(h_i\) enter \(a_P\) and \(b\) at the boundary nodes exactly as the boundary conditions prescribe.
Layer-to-mesh mapping (1D): the physical layers are mapped onto the uniform mesh by cumulative coordinates. Each cell takes for reference the material containing its centre, and its volumetric heat capacity is the thickness-weighted average of \(\rho c\) over the cell, so the total thermal mass \(\sum_j \rho c_j L_j\) is conserved exactly at any \(N_x\).
At interior faces, the conductance (per unit area) is built from the series resistance between the adjacent cell centres (Patankar 1980):
\[ G_f \;=\; \left[\, \int_{x_P}^{x_E} \frac{dx'}{k(x')} \,\right]^{-1} \;=\; \left(\, \sum_m \frac{\delta_m}{k_m} \,\right)^{-1}, \tag{2}\]
where the sum runs over the material segments crossed by the span. Material interfaces may therefore fall anywhere inside a cell, and layers thinner than \(\Delta x\) (metal sheets, membranes) contribute their exact resistance and mass at any position and any \(N_x\): a 0.5 mm metal sheet adds its true, negligible \(L/k\) instead of being deleted or smeared. When the interface lies exactly on the face with equal half-widths \(\delta_L = \delta_R = \Delta x / 2\) and perfect thermal contact, Equation 2 reduces to the familiar harmonic mean:
\[ k_f \;=\; \frac{2\, k_L\, k_R}{k_L + k_R}, \qquad G_f = \frac{k_f}{\Delta x}. \tag{3}\]
Solution algorithm
In 1D the discrete system is tridiagonal, and each time step is solved exactly in one pass of the Tridiagonal Matrix Algorithm (TDMA / Thomas algorithm). No inner iteration is needed.
In 2D each time step is solved line by line (Patankar 1980): a TDMA solves each line across the width implicitly, while the through-thickness neighbour terms and the non-linear cavity physics (the convective coefficient \(h_c(\Delta T)\) and the \(T^4\) radiation exchange) are evaluated with the previous sweep’s field (Jacobi lagging) and re-assembled at the start of every sweep.
Sweeps repeat until two non-cancellable criteria hold simultaneously (both in °C, config2d.tol_inner, \(10^{-8}\) by default):
\[ \max_P \left| T_P^{(r+1)} - T_P^{(r)} \right| \;\le\; \varepsilon_{in} \qquad\text{and}\qquad \max_P \frac{\left| a_P T_P - \textstyle\sum a_{nb} T_{nb} - b \right|}{a_P} \;\le\; \varepsilon_{in}, \tag{4}\]
i.e. the largest node update of the last sweep and the largest scaled residual of the discrete equations (Equation 1), the residual being evaluated with the accepted field’s own (re-assembled, non-linear) coefficients. This is the convergence monitor recommended by Patankar (1980). A per-step cap config2d.max_inner (\(10^4\) sweeps) guards against stalls; if it is hit, the solve completes but is flagged converged = False (see below). The default \(\varepsilon_{in} = 10^{-8}\) °C keeps the worst-case accumulated inner error over one day (8640 steps × \(10^{-8}\) ≈ \(10^{-4}\) °C) below the day-to-day tolerance tol_day.
Indoor-air coupling (free-running mode)
The conduction solver is implicit and therefore unconditionally stable; no time-step restriction comes from the wall itself. The coupling with the indoor air is, however, segregated: each step solves the wall using the indoor temperature of the previous step, and then advances the air node explicitly (forward Euler over the air balance):
\[ T_i^{new} \;=\; T_i \;+\; \mathrm{Fo}_a \left( \overline{T}_{surf} - T_i \right), \qquad \mathrm{Fo}_a = \frac{h_i\, \Delta t}{\rho_a\, c_a\, L_a}. \tag{5}\]
That scalar explicit step is the one piece of the scheme with a stability constraint: \(\mathrm{Fo}_a < 2\) for stability and \(\mathrm{Fo}_a \lesssim 1\) to avoid oscillations. With the fixed \(\Delta t = 10\) s, \(\mathrm{Fo}_a \approx 0.03\), far inside the stable region. This is why config.dt is locked: a user-enlarged \(\Delta t\) would push the air node toward non-physical behaviour while buying almost no speed (the inner iterations grow to compensate).
In the 2D solver the cavity-air nodes \(T_h\) are advanced the same way (forward Euler), and their stability parameter is not uniformly small:
\[ \lambda_h\,\Delta t \;=\; \frac{h_c\, P_{cav}\, \Delta t}{\rho_a\, c_a\, A_{cav}}, \tag{6}\]
with \(P_{cav}\) and \(A_{cav}\) the cavity perimeter and cross-section. It depends on the geometry and on the instantaneous \(h_c\): at \(\Delta t = 10\) s the documented example geometries give \(\lambda_h\Delta t \approx 0.2\) for the hollow-block wall cavity but up to \(\approx 1\) for the roof-slab cavities, stable (\(< 2\)) though not always in the monotone regime (\(\le 1\)).
In the air-conditioned mode the indoor node is a fixed setpoint, so Equation 5 does not apply to it. The scheme is still not purely implicit there: with Fill.AIR the cavity-air nodes \(T_h\) keep advancing explicitly (Equation 6 holds in AC mode too).
Convergence to the periodic regime
The boundary forcing (\(T_{sa}\)) is periodic with a 24 h period, so the physical solution of interest is the periodic one. Starting from a uniform field at the neutrality temperature \(\overline{T_n}\), the day is simulated repeatedly, each day starting from the previous day’s final state, until every persisted state stops changing between consecutive days: the solid field, the indoor-air node \(T_i\) (free-running mode) and each cavity-air node \(T_{h,c}\):
\[ \text{day\_error} \;=\; \max\!\left( \frac{1}{N}\sum_{\text{nodes}} \left| T^{day+1} - T^{day} \right|,\; \left|\Delta T_i\right|,\; \max_c \left|\Delta T_{h,c}\right| \right) \;\le\; C_s. \tag{7}\]
The tolerance is \(C_s = 5\times10^{-4}\) °C in both solvers (hard-coded in 1D; config2d.tol_day in 2D), and both cap the iteration: config2d.max_days (60) in 2D and the module constant MAX_DAYS = 60 in 1D. Both solvers report the outcome explicitly: after a solve, days holds the days used, day_error the final value of Equation 7, and converged is True only when the day cycle met the tolerance (in 2D, additionally, only if every inner step converged within max_inner; inner_iterations records the largest sweep count). A non-converged solve emits a RuntimeWarning instead of returning a result indistinguishable from a converged one.
As an independent physical check, the free-running solvers also expose energy_imbalance \(= |Q_{in} - Q_{out}| / \max(Q_{in}, Q_{out})\): in a perfectly periodic regime the gross energy entering the indoor node over the day equals the gross energy leaving it, so this ratio should be close to zero.
Differences with the 2016 paper
The reference paper describes the online tool Ener-Habitat (v2.2). This package solves the same physical problem but is not a copy of that implementation; the table lists the deliberate differences (each entry verified against the code).
| Aspect | Online tool (Barrios et al. 2016) | This package |
|---|---|---|
| Climate source | Built-in database, 80 Mexican cities (Meteonorm) | User-provided EPW file, any location |
| Average-day irradiance | Positive half-sine, amplitude = monthly max; diffuse varies linearly with tilt | Monthly hour-by-hour means interpolated in time; tilted-surface \(I_s\) by pvlib transposition |
| Long-wave factor \(RF\) | 3.9 °C horizontal → 0 vertical, linear in tilt | Same (linear in tilt) |
| Time step | 1 s | 10 s, fixed (explicit air-node coupling, Equation 5) |
| Mesh | 1D \(N_x = 100\); 2D \(160 \times 160\) | 1D \(N_x = 200\); 2D \(80 \times 160\) (config2d) |
| Indoor film coefficient \(h_i\) | 8.1 (wall); 9.4 / 6.6 (roof, by flux direction) | Same convention: 8.1 (wall), 9.4 / 6.6 selected per time step on roofs (hi_flow); all values configurable, and hi_flow = False restores a fixed 8.1 |
| Periodic tolerance | \(C_s = 10^{-5}\) °C | \(C_s = 5\times10^{-4}\) °C (+ max_days cap in 2D) |
| 2D geometries | Filled or hollow block (one cavity) | HollowBlock (walls) and Slab (roofs: N cavities, L-shaped rib), fill AIR/SOLID |
| AC setpoint | Neutrality temperature \(T_n\) | Same (\(T_n\) constant); DeltaTn provided as data only |
| Materials | Built-in + user database | User materials.ini only (no bundled data) |
Validation
The validation record combines the published validation of the method with the package’s own regression tests:
- 1D vs EnergyPlus (Barrios et al. 2016): six constructive systems (massive, insulated and mixed), non air-conditioned and air-conditioned, hot month. Maximum differences: indoor temperature 0.5 °C, decrement factor 0.03, lag time 0.4 h, cooling energy 3%, heating energy 4%; both programs rank the systems identically.
- 2D, filled block: setting the same material everywhere, the 2D solution equals the 1D one (Barrios et al. 2016). The package reproduces this 2D → 1D reduction as a test.
- 2D, hollow block vs experiment: steady-state heat flux through a concrete hollow-block wall compared with hot-box measurements (Borbón et al. 2010); maximum difference 6.5% (experimental uncertainty ±3.8%) (Barrios et al. 2016).
- Package vs legacy code: the 2D solver was ported from the C code of the online tool;
tests/golden/holds golden master dumps of the C solver and the suite (tests/test_eh2d_*.py) checks the port reproduces them, plus geometry, coefficient, single-step, full-day and AC tests.
Performance
Each solve is single-threaded, JIT-compiled with numba (the first call in a session pays the compilation cost). In-solver thread parallelism was measured (~1.06× on realistic meshes) and removed as not worthwhile. For parameter sweeps over many independent configurations, parallelise at the process level (multiprocessing, joblib): measured speed-ups are near-linear, ≈6× with 8 processes and ≈10× with 16 on an 18-core machine.