Description

This PR introduces a MultiLayer 3D Thermal Architecture to PyBaMM, enabling true thermal resolution through the stack thickness ($x$-direction) of pouch cells. Previously, thermal domains were shared or averaged across the stack. This update assigns independent 3D temperature fields to each electrochemical layer within the stack, allowing the model to capture critical thermal stratification—such as localized hotspots generated by individual aged or high-resistance layers.Key features include:

• Layer-Specific Thermal Feedback: Each layer's electrochemistry (SPM, SPMe, or DFN) responds to its own local volume-averaged temperature.
• Thermal Interface Physics: Configurable boundary conditions between layers to model either perfect thermal contact or realistic thermal contact resistance (e.g., from adhesives or air gaps).
• Flexible Electrical Coupling: Built-in algebraic constraints to natively support both parallel and series electrical connections across the stack layers.

Motivation

Background

In modern battery packs, pouch cells are often stacked together to form modules. Similarly, each cell may be composed of many layers of parallell (or series) connected unit stack layers. When there are multiple stack layers in a cell, or multiple cells stacked:

• Heat is generated differentially across layers due to varying local conditions
• Temperature gradients develop through the stack thickness (perpendicular to the pouch faces)
• Thermal contact resistance between layers creates additional thermal impedance
• Asymmetric cooling (e.g., cold plate on one side) creates heterogeneous conditions within stack layers.

Limitation of Existing Models

PyBaMM's existing Basic3DThermalSPM model treats the entire cell as a single 3D thermal domain. While it captures in-plane temperature gradients (within x-y-z coordinates of one cell), it cannot resolve:

1. Independent electrochemical states in different physical layers
2. Through-stack temperature gradients when layers are thermally isolated or poorly coupled
3. Layer-specific current distribution in parallel-connected stacks
4. Inter-layer thermal contact resistance effects

For example, in a stack of 120 unit cells with asymmetric cooling (cold plate on one side, insulated on the other), the cells closest to the cold plate will be significantly cooler than those far away. This affects:

• Local reaction kinetics (temperature-dependent exchange current density)
• Current sharing in parallel connections (hotter layers draw less current)
• Capacity fade and aging (non-uniform temperature accelerates degradation in hot zones)

Possible Implementation

This feature takes the thermal model formulation in the src/pybamm/models/full_battery_models/lithium_ion/basic_spm_with_3d_thermal.py and applies a different formulation for the temperature in a case where it is a multilayer pouch cell with many layers.

The math is outlined as follows.

For each layer $i$, the temperature field $T_i(x, y, z, t)$ satisfies the heat equation:

$$ ρc_p ∂T_i/∂t = ∇·(λ∇T_i) + Q_{vol}^i $$

where $Q_{vol}^i$ is broadcast from the electrochemical model to the 3D mesh.

Volume-averaged temperature (couples 3D thermal to 0D electrochemistry):

$$ T_av^i = (1/V_i) ∫∫∫_{V_i} T_i(x,y,z) dx dy dz $$

This is enforced as an algebraic constraint linking the SPM variable $T_{av}^i$ to the spatial average of $T_i$.

Adjacent layers $i$ and $i+1$ exchange heat through their shared interface. The heat flux $q$ [W.m-2] is driven by the temperature difference across the thermal contact resistance $R_{th}$:

$$ q_{interface}^{i→i+1} = \frac{(T_i|_{x=x_{max}} - T_{i+1}|_{x=x_{min}})}{R_{th}} $$

This is implemented as matching Neumann boundary conditions:

$$ -λ_i (∂T_i/∂x)|_{x=x_{max}} = q_{interface}^{i→i+1} $$

$$ -λ_{i+1} (∂T_{i+1}/∂x)|_{x=x_{min}} = -q_{interface}^{i→i+1} $$

The outermost layers have convective cooling on their exposed faces:

$$ -\lambda \left. \frac{\partial T}{\partial x} \right|_{x=0} = h_{\text{edge}, x_{\text{min}}} (T - T_\infty) \quad \text{[layer 0, left face]} $$

$$ -\lambda \left. \frac{\partial T}{\partial x} \right|_{x=L} = h_{\text{edge}, x_{\text{max}}} (T - T_\infty) \quad \text{[layer N-1, right face]} $$

Similarly for y and z faces (all layers):

$$ -\lambda \left. \frac{\partial T}{\partial y} \right|_{y=0, W} = h_{\text{edge}, y} (T - T_\infty) $$

$$ -\lambda \left. \frac{\partial T}{\partial z} \right|_{z=0, H} = h_{\text{edge}, z} (T - T_\infty) $$

where $h_edge$ is the heat transfer coefficient [W.m-2.K-1].

Asymmetric cooling can be modeled by setting different h values on different faces.

Electrical Coupling

Parallel Connection
In parallel, all layers share the same voltage but carry different currents:

$$V_0 = V_1 = \dots = V_{N-1} = V_{\text{terminal}}$$

$$I_{\text{total}} = I_0 + I_1 + \dots + I_{N-1}$$

Introducing current fraction $f_i$:

$$I_i = f_i \cdot \frac{I_{\text{total}}}{n} \quad (\text{where } n = \text{layers per zone})$$

$$\sum_{i=0}^{N-1} f_i = 1$$

The voltage equality is enforced via algebraic constraints:

$$V_i - V_0 = 0 \quad \text{for } i = 1, 2, \dots, N-1$$

Series Connection
In series, all layers carry the same current but have different voltages:

$$I_0 = I_1 = \dots = I_{N-1} = \frac{I_{\text{total}}}{n}$$

$$V_{\text{terminal}} = V_0 + V_1 + \dots + V_{N-1}$$

No algebraic constraints are needed; current is directly imposed.

Model Formulation

When initiating the model, there are a few new parameter to control how it is setup.

Parameter	Symbol	Meaning
num_physical_layers	N_phys	Total physical unit cells in the real stack
num_subdivisions	N_zones	Number of thermal/electrochemical zones modeled
layers_per_zone	n = N_phys / N_zones	Physical cells lumped per zone
mesh_h	h	FEM mesh characteristic length (controls resolution)
connection	—	Electrical connection: "parallel" or "series"
R_contact	R_th	Inter-layer thermal contact resistance [K.m2.W-1]

The model solves N_zones independent SPMs/SPMes/DFNs, each representing n physical cells. Electrode thicknesses and per-unit-cell overpotentials are preserved; only the thermal extent and stack capacity are rescaled.

Scaling for Physical Layers

When layers_per_zone = n > 1, each modeled zone represents n physical cells:

1. Thermal extent: Zone thickness = n × (single cell thickness)
2. Capacity: Stack capacity = N_phys × (single cell capacity)
3. Current per SPM: I_spm = I_zone / n (preserves per-cell current density)
4. Electrode thicknesses: L_n, L_p unchanged (per-cell geometry preserved)

The apply_stack_scaling() method automatically multiplies the nominal cell capacity by N_phys to match the physical stack.

Variables

Variables for different stack layers are denoted by prefixing Layer i before the variable name.

Additional context

No response