On the coupled thermo–electro–chemo–mechanical performance of structural batteries with emphasis on thermal effects

Carbon fibre (CF) based structural batteries is a type of battery designed to sustain mechanical loads. In this paper, a fully coupled thermo–electro–chemo–mechanical computational modelling framework for CF-based structural batteries is presented. We consider the combined effects of lithium insertion in the carbon fibres leading to insertion strains, and thermal expansion/shrinkage of the constituents leading to thermal (free) strains, while assuming transverse isotropy. The numerical studies show that the developed framework is able to capture the coupled thermo–electro–chemo–mechanical behaviour. Moreover, it is found that the dominating source for heat generation during galvanostatic cycling is associated with discontinuities in the electrical and chemical potentials at the fibre/electrolyte interface. Further, a limited parameter study shows that the temperature change during electrochemical cycling is significantly influenced by the applied current, thermal properties of the constituents and heat exchange with the surroundings. Finally, for large temperature variations, e.g. as identified during relevant (dis)charge conditions, the magnitude of the thermal strains in the structural battery electrolyte (SBE) are found to be similar to the insertion induced strains.


Introduction
Lithium-ion (Li-ion) batteries have been used for energy storage in electric vehicles and portable devices since the mid-'90s (Cano et al., 2018). However, the performance of conventional Li-ion batteries is limited by their energy storage to weight ratio (Cano et al., 2018;Armand et al., 2020). For example, to realize fully electric regional aircraft, this ratio needs to be increased compared with today's standard (Schäfer et al., 2019). A potential route to improve this ratio is to develop energy storage solutions with the ability to sustain mechanical loads Federico et al., 2021;Moyer et al., 2020;Johannisson et al., 2018;Yin et al., 2020;Wang et al., 2019;Lutkenhaus and Flouda, 2020;Asp et al., 2019;Ladpli et al., 2019). One such solution exploits electrodes consisting of carbon fibres (CFs) embedded in structural battery electrolyte (SBE). This type of material can be used in so-called ''structural batteries'', i.e. batteries with mechanical load-bearing capability Moyer et al., 2020;Johannisson et al., 2018;Yin et al., 2020) (see schematic illustration of the conceptual design of the laminated structural battery in Fig. 1a), in structural electrochemical actuators (Johannisson et al., 2020), or in structural energy harvesting and strain sensing materials (Jacques et al., 2013b;Harnden et al., 2022). * Corresponding author.
The SBE in the studied material is a bi-continuous bi-phasic composite which consists of two phases: (i) A solid phase, corresponding to a mechanically robust porous polymer network with an open pore system; (ii) A liquid phase which contains liquid electrolyte with Lisalt (cf. Ihrner et al. (2017) and Schneider et al. (2019)). The liquid phase in the porous polymer network enables ion transport between the electrodes, while the solid phase makes it possible to distribute mechanical loads. A schematic illustration of the SBE is presented in Fig. 1b, which corresponds to a numerically generated nano-structure, from Tu et al. (2020), representing an idealized fine-scale geometry of the SBE. The characteristic pore size of the polymer network is roughly in the order of 50-200 nm (Schneider et al., 2019). This can be compared with the typical diameter of a carbon fibre, which is approximately 5 μm.
It is well known that the electrochemical performance of conventional batteries is highly influenced by the temperature, see e.g. Al Hallaj et al. (1999), Thomas and Newman (2003) and Forgez et al. (2010). Moreover, the mechanical properties of conventional CF reinforced polymers (as construction material) are known to be highly affected by temperature (Gigliotti et al., 2011). Due to the fact that electrochemical cycling of a battery generates (or absorbs) heat, associated with e.g. transport of charge species and redox reactions, the  Hagberg et al. (2018) and Sanchez et al. (2021). (b) Numerically generated porous nano-structure, from Tu et al. (2020), representing an idealized fine-scale geometry of the SBE that consists of a polymer skeleton (solid phase) saturated with a liquid electrolyte phase. (c) Illustration of the developed thermo-electro-chemo-mechanical framework, and comparison with previous work (Carlstedt et al., , 2022b. temperature inside a battery will change during operation/service. This will alter the temperature dependent properties and induce thermal expansion/shrinkage of the constituents. Hence, it is crucial to be able to predict the distribution and evolution of temperature during operation as part of evaluating the multifunctional (i.e. coupled thermal, electrochemical and mechanical) performance of the structural battery material.
For conventional Li-ion batteries, a variety of approaches for modelling the thermo-electro-chemo-mechanical problem have been proposed, see e.g. Hu et al. (2017), Zhang et al. (2020), Duan et al. (2018) and Wang et al. (2021). However, only a few are developed based on a rigorous thermodynamics setting. In this context, we note an important contribution from Salvadori et al. (2018), who developed a fully coupled modelling framework for mass and heat transport, mechanics, and chemical reactions with trapping (which refers to the fact that only a fraction of the total available mass of mobile species is transported, whereas a significant counterpart remains immobile), based on non-equilibrium rational thermodynamics and small strain kinematics.
With respect to computational modelling of structural batteries, studies of various interactions between the electrochemical and mechanical fields have been carried out by the authors  and by Xu et al. (2018a,b). Moreover, Tu et al. (2020) studied the bifunctional performance of the SBE on the nano-scale numerically. Goudarzi et al. (2022) have developed an efficient computational approach for modelling and simulation of electrochemical phenomena taking place in (three-dimensional) fibrous battery electrodes. Further, Pejman et al. (2021) developed a multi-physics design optimization framework for structural batteries and Yin et al. (2020) compared experimental and numerical results for modified CF electrodes in liquid electrolyte. In Yin et al. (2020), and in the accompanied work by Hong et al. (2021), the computational framework is based on the battery modelling scheme originally developed by Newman and coworkers (Newman and Tiedemann, 1975;Doyle et al., 1993;Newman and Thomas-Alyea, 2004) that presumes one-way coupling between the electrochemical and mechanical processes. Recently, the authors (Carlstedt et al., 2020) developed a thermodynamically consistent computational modelling framework to study the electro-chemo-mechanical properties of structural batteries while allowing for two-way coupling between the electrochemical and mechanical fields (Fig. 1c). In addition, the contribution to the ion transport from convection, i.e. from seepage of the liquid phase of the SBE, was accounted for in Carlstedt et al. (2022b), and numerical predictions were compared favourably with experimental data in Carlstedt et al. (2022a). However, in these studies isothermal conditions were assumed. With respect to the thermal interaction on the performance of structural batteries, only limited evaluation has been conducted. For example, Schutzeichel et al. (2019Schutzeichel et al. ( , 2021 studied, both experimentally and numerically, the thermal behaviour of multifunctional composite materials made from polymer electrolyte coated carbon fibres. Moreover, Carlstedt and Asp ) developed a semi-analytical framework to predict the mechanical consequences (thermal and diffusion induced stresses) from electrochemical cycling. In , the electrochemical analysis was highly simplified and the analysis was limited to one-way coupling between the electrochemical and thermal analysis and the mechanical analysis. Hence, a complete modelling framework which accounts for the fully coupled thermal, electrochemical and mechanical processes in structural batteries is currently lacking.
In this paper we extend a previously developed (thermodynamically consistent) computational framework for structural batteries, , to account for thermal interaction effects (cf. Fig. 1c). The energy balance then becomes an additional governing equation, and the fully coupled thermo-electro-chemo-mechanical problem is solved for a simplified geometry with the appropriate constitutive relations and interface/boundary conditions. Numerical studies are performed to demonstrate the capability of the developed framework and to evaluate the effects of different coupling parameters. In particular, the significance of the individual terms linked to heat generation (or absorption) during galvanostatic (dis)charge conditions are evaluated, and a parameter study is carried out in order to assess the (additional) effect of the thermal properties on the combined thermal-electro-chemical-mechanical performance.

Studied half-cell and coupled analysis
To simplify the analysis we study the CF-SBE electrode half-cell illustrated in Fig. 2a. This simplification is motivated by the similarities of the two electrodes in the laminated structural battery architecture (both electrodes consist of fibres embedded in SBE), cf. Carlstedt et al. ( , 2022b and Fig. 1a. The basic design principle of the laminated structural battery is to create a laminate by stacking laminae with different functionalities. The design presented in Fig. 2a has the following laminae/components: (1) The working electrode, made from CFs embedded in SBE; (2) A separator layer, e.g. made from a thin layer of SBE; (3) The counter (and reference) electrode made from Li-metal.
The flow of charged species within the studied battery cell is schematically illustrated in Fig. C.10a-b in Appendix C. We note that the electrons in the fibres are assumed to flow in the fibre direction only. Further, we assume a linear distribution of the electronic current along the fibre (illustrated in Fig. C.10c). This leads to the idealized conditions where the current input is the same at all cross-sections along the fibre direction, which allows us to simplify the analysis and only study one cross-section (2D-analysis). Moreover, the electron flow in the Li-metal as well as the current collections, (electron) conductive adhesive and external circuit are neglected (for simplicity).
The idealized material representation of the studied half-cell is presented in Fig. 2b and corresponds to a repeatable unit in the 1direction of the CF-SBE electrode lamina. It should be noted that the separator is excluded (for the sake of simplicity) but can easily be added as an additional electrolyte phase, as e.g. done is Carlstedt et al. (2022a). Moreover, we note that the SBE is treated as a continuum at the studied length-scale and effective properties are utilized (cf. Carlstedt et al. (2022b)). The introduced notation is defined in Section 3.
Finally, a schematic illustration of the coupled thermo-electrochemo-mechanical problem is presented in Fig. 2c. The pertinent balance equations are linear momentum balance, charge balance (Gauss law), mass balance, and energy balance. The basic independent fields are the displacement field , the electrical potential , the ion concentrations and chemical potentials of species =Li,X, and the (absolute) temperature .

Time-continuous strong format and modelling assumptions
The time-continuous strong format and modelling assumptions for the individual domains, interfaces and boundaries are presented in this section. We thereby extend the theory, as compared with previous work by the authors , to account for thermal effects associated with the electrochemical and mechanical processes. Moreover, a detailed derivation of the thermodynamics for a generic electro-chemical-thermal-mechanical system (not necessarily the actual battery) is presented in Appendix A.

Fibre domain(s)
The fibre domain(s) are denoted f = ∪ fibres f, (cf. Fig. 2b). The following special assumptions are introduced: (i) Material properties are characterized as transversely isotropic (isotropy pertains to the cross-section, defined by Cartesian coordinates 1 , 2 ); (ii) The single active species is Li, which moves into (and out from) the fibre; (iii) The fibres ( f ) and the positive connector ( + ) are connected via an external circuit. It is assumed that the electric potential is uniform (and equal) in all fibres, ( , ) = − ( ), for ∈ f = ∪ fibres f, ; consequently, the electric field is zero within the fibre domain. Although the analysis is two-dimensional (restricted to the 1 , 2 -plane), so-called resistive (Joule) heat generation along the fibres is accounted for in an approximate fashion.

Balance equations
The governing balance equations in the fibre domain(s) ( f ) are expressed in the strong format as follows: where is the (symmetric) stress tensor, Li is the ion flux vector for Li, is the heat flux vector, Li is the ion concentration of Li, 1 Li is the chemical potential of Li, and is the (absolute) temperature. The generated heat due to the electric current in the fibre is denoted elec (see Section 3.1.3). Moreover, is the bulk density and v is the volume-specific heat capacity.
Remark 1. The contributions to the energy balance with respect to and Li , denoted and Li , are derivable theoretically in a thermodynamically consistent fashion, cf. Appendix A. However, we have simplified the energy equation from the outset by setting = and Li = 0. To ignore is standard for quasi-static mechanical loading. Further, Li is often neglected in the battery modelling literature (see e.g. Hu et al. (2017)).

Constitutive relations
In order to obtain explicit expressions of the purely energetic quantities and Li , we introduce the following partial free energy densities me-th-ch ( , , Li ) = 1 2 where 0 Li is a reference/standard value of the chemical potential and is the universal gas constant. As to the introduced variables, [ ] is the (small) strain tensor expressed as a linear operator of the displacement field , and̃L i ∶= Li Li,max is the normalized mass concentration of Li (with respect to the assumed maximum Li-concentration Li,max ). It should be noted that we have assumed that the specific heat ( v ) and the elastic tensor ( ) are temperature independent (for simplicity). Further, Eq. (2b) is expressed such that the reference state for the generic/idealized model representation) of the studied half-cell. In the external circuit, Load represents electric loading. (c) A schematic illustration of the coupled thermo-electro-chemo-mechanical problem. The introduced notation is defined in Section 3. chemical potential is defined at = , Li = Li,0 and = 0 . Moreover, ch and th are the insertion and thermal (free) strains, defined as where the reference/initial values Li,0 and 0 denote the state of a pristine material, i.e. the state at which no chemical or thermal strains are present in the material, respectively. Further, ch and th are second order tensors containing the coefficients of the lithium insertion and thermal induced expansion of the fibres, respectively, and they are defined as 2 where ch ⟂ , th ⟂ and ch ∥ , th ∥ are the transverse and longitudinal expansions coefficients associated with lithium insertion and thermal expansion, respectively. Finally, as to the elasticity tensor that is pertinent to transverse isotropy, it is defined by five independent parameters 2 ∶= ⊗ is the i:th base dyad. as follows: where = 1 + 2 + 3 is the 2nd order identity tensor, sym ∶= 1 2 [ ⊗ + ⊗ ] is the (symmetric) 4th order identity tensor, 3 whereas ∶= 1 4 [ 3 ⊗ + 3 ⊗ + ⊗ 3 + ⊗ 3 ] is a 4th order symmetric tensor. Lamé's first parameter is denoted , the shear modulus is denoted , whereas ∥ is the uniaxial strain modulus. Since measurements have shown that it is mainly the transverse elastic modulus of the carbon fibre that is significantly influenced by lithiation , ⟂ ( Li ) is the only elastic modulus that is chosen to depend on the Li-concentration.
We are now in the position to compute It appears that the initial values of and Li are en ( , Li,0 , 0 ) = and en Li ( , Li,0 , 0 ) = 0 Li . Next, we consider the ''dissipative'' variables Li and , whose constitutive relations are chosen as where we introduced the mobility tensor and thermal conductivity tensor as follows: The mobilities in the transverse and longitudinal directions are denoted Li,⟂ and Li,∥ , respectively. Similarly, the thermal conductivity coefficients in the transverse and longitudinal directions are denoted ⟂ and ∥ , respectively.

Generated heat in the carbon fibre
The (only) source term of heat supply in Eq. (1c), denoted elec , stems from the power that is dissipated due to electric current along the fibres, whereby we recall that the fibres work as both electrodes for the ionic transport (diffusion) as well as current collectors. The question is how to compute this term in a two-dimensional analysis? In a proper three-dimensional analysis the situation is more obvious, since the effect of any current in the energy equation can be written as elec = − e ⋅ , cf. Eq. (A.20), where − e is the (part of the) current that is due to electronic transport. To motivate the proper expression of elec in the present two-dimensional analysis, we thus argue as follows: First, we recall the assumption that the electric potential is constant ( = − ) within each fibre cross-section, collectively occupying the domain f . This means that there is no electric current in the crosssectional plane. However, in reality must vary along the fibres, giving rise to an electric field component f ∶= ( ) 3 = − 3 and, thus, a resulting current f ∶= ( − e ) 3 . Further, we assume that the constitutive equation for f is f = f f , where f is a (constant) electric conductivity. Hence, we obtain the simple expression It only remains to determine the value of either f or f . In the particular case of galvanostatic control with a prescribed supplied current − via the fibre perimeter per unit length of the fibre [A/m] (which is the case studied in this paper), we assume that f is equal in each fibre and that the current varies linearly along the fibres. Further, we approximate the generated heat at the evaluated cross-section as the average value of the total elec (with respect to the total fibre length ). This corresponds to assuming an average current of where f is the total cross-sectional area of all fibres and is the fibre length (cf. Fig. 2a). This gives the approximated expression elec = We note that the total current extracted from the cell then becomes − . In a potentiostatic problem on the other hand, the value − is part of the solution, i.e. Eq. (11) has to be evaluated as part of the problem formulation. It should be noted that we assume a finite electric conductivity of the fibres when estimating the generated heat in the fibres, while we assume a constant electric potential along the fibre (as discussed in Section 2) for the cross-sectional problem. We note that this is not consistent with the energy conservation of the system. However, this assumption is made to simplify the analysis in order to allow for utilizing a 2D-model representation while estimating the temperature increase during operation. Without this assumption, the potential at the fibre would be different at different cross-sections.

Electrolyte (or SBE) domain
The electrolyte (or SBE) domain is denoted e (cf. Fig. 2b). The following special assumptions are introduced: (i) Material properties are characterized as isotropic; (ii) The Li-ions are positively charged (cation) and the companion X-ions (anion) are negatively charged; (iii) The current density is carried both by Li + and the companion anion X − .

Balance equations
The governing balance equations in the electrolyte domain ( e ) are expressed in the strong format as follows: In addition to those variables already introduced for the fibre domain, we introduced the following notation: is the electric flux density vector (dielectric displacement), Li and X are the ion concentrations of Li and X, respectively (defined as ion mass, in moles, per unit mass of fluid, in kg). Moreover, X is the chemical potential and X is the ion flux vector for the species X. Finally, F is the intrinsic density of the fluid (electrolyte) in the SBE and is Faraday's constant.
Remark 3. In analogy with the situation for the fibres, we have simplified the energy equation from the outset by setting = , Li = 0 and X = 0, cf. Appendix A.

Constitutive relations
In order to obtain explicit expressions of the purely energetic quantities , , Li , and X , we introduce the following partial free energy densities where ∶= − is the electric field. 4 The only free strain in the electrolyte is the thermal induced strain, which is given as where th ∶= th is the isotropic thermal expansion tensor. The standard isotropic elasticity tensor is expressed in terms of Lamé's parameters as where is Lamé's first parameter and is the shear modulus. Moreover, the isotropic permittivity tensor is given as: where 0 and are the vacuum and relative permittivity, respectively. In analogy with the situation for the fibre domain, we define the normalized ion concentration of species Li and X as̃L i = Li We are now in the position to compute It appears that the initial values of , Li and X are en ( , 0 ) = , en Li ( Li,0 , 0 ) = 0 Li and en X ( X,0 , 0 ) = 0 X . The dissipative variables Li , X and are chosen as where the isotropic mobility and thermal conductivity tensors are defined as: In these expressions, is the mobility coefficient of the species for = Li,X. Obviously, we have introduced the simplification that there is no (constitutive) coupling between the diffusion of Li + and X − . Moreover, Eqs. (19a)-(19b) differ from Eq. (8a) principally in the sense that the mobilities vanish when = 0 or = ,sat , cf. Kaessmair et al. (2021). Further, is the isotropic thermal conductivity.
Finally, using Faraday's law of electrolysis (while exploiting (18a) and (18b)), we derive the constitutive relation for the current density: where we introduced the ionic conductivity  ∶= 2 [ Li + X ].

Interior boundaries
The fibre/electrolyte interfaces, denoted fe , represent interior boundaries. We assume that the displacement field is continuous across fe (i.e. perfectly bonded). However, due to the redox reactions in the close vicinity of the interface, Li , and may be discontinuous across fe , whereas Li, is continuous across fe . The discontinuities in Li and are modelled as linear interface relations involving the ''electric resistance'' which turns to have a model structure that is similar to a linearized Butler-Volmer relation. With respect to the ion flux we assume that Li, (= e Li, ) ∶= Li ⋅ , where is the normal on fe pointing from the electrolyte domain e into the fibre domain f , is governed constitutively by an interface mobilitȳsuch that In Eq. (21), we introduced the jump operator (commonly referred to as the electrochemical potential 6 ), where e ∶= ( e ) is evaluated in the electrolyte at the fibre-electrolyte interface. The current density flux ∶= ⋅ across the interface fe is defined as where we introduced the assumption X, ( e ) = 0, ∈ fe , i.e. the transport of X − is blocked at the fibre-electrolyte interface. Further, we introduced the interface ionic conductivity ∶= 2̄. Next, we introduce the constitutive assumption for (= e ) ∶= ⋅ where is the interface permittivity. With respect to the thermal conditions, we evaluate the thermal energy balance at the interface (cf. Eqs. (1c) and (12e)). Firstly, we introduce the normal components of the heat flux vector as f ∶= f ⋅ f and e ∶= e ⋅ e , where (once again) we note that e = − f ∶= . The balance equation including the interface may then be expressed as where we used that e Li, = − f Li, ∶= Li, and e = − f ∶= , i.e. the normal component of Li , as well as of , are continuous across the interface. Now, inserting the ansatz into (24), we obtain where we used the identity = Li, . The quantity represents the heat absorbed at the interface (i.e. ''heat sink''), and it is noted that the expression in Eq. (26) follows directly from the heat balance and does not require any additional constitutive assumption. This is in line with conventional battery modelling frameworks, see e.g. Newman and Thomas-Alyea (2004).
Finally, we adopt a standard expression for the heat transfer across the interface (from the fibre to the electrolyte) in terms of the temperature jump where th is the interface transfer resistance coefficient (or effective layer conductance).
Remark 4. The interface's contribution to the dissipation inequality may be included in the formulation, e.g. in the case of mechanical damage (cf. Singh and Pal (2020) and Rezaei et al. (2021)).

Exterior boundary
The entire external boundary comprises ext = ∑ 3 =1 ∪ ext, and the collector boundary + . With respect to the mechanical conditions, related to displacements ( ) and tractions ( = ⋅ ), we assume that: These conditions are motivated by the fact that the studied (part of the) lamina will in practice constitute a layered plate structure (cf. ). Further, for chemical conditions related to ion flux ( Li, = Li ⋅ , X, = X ⋅ ), we set: Li, = 0 on ext,2 ∪ ext,3 (29b) The assumed chemical conditions are motivated by the fact that the ion flux occurs between the Li-metal and fibres (i.e. mainly in the 2direction) and that the height of the studied unit corresponds to the thickness of the electrode lamina, cf. Fig. 2b. With respect to the electrical conditions, related to the electric flux density ( ), we set: where ext denotes all exterior boundaries except + . It is noted that + is a spatially constant value in the collector (Li-metal) just outside + and e is the electrolyte potential along + . Moreover, we assume that + is prescribed at 0 V.
Finally, the thermal conditions related to the heat flux ( = ⋅ ) are defined as: In Eq. (31b) we assume a convective heat flux where th ext is the heat transfer coefficient and ext is the external temperature. These conditions are motivated from experimental studies, Asp et al. (2021) and Johannisson et al. (2018). Finally, it should be noted that we neglect heat that is generated or absorbed along + .

FE-model geometry and loading conditions
The model geometry and (triangular) finite element (FE) mesh are presented in Fig. 3a. An SEM-image of the corresponding cross-section (from Carlstedt et al., 2022a) is provided in Fig. 3b, for comparison. The fibre volume fraction ( f ) is set to f = 0.45, and the dimensions and fibre volume fraction are motivated by experimental measurements in Johannisson et al. (2018) and Carlstedt et al. (2022a).
The numerical implementation is done in the commercial FE software COMSOL Multiphysics version 5.4. The time-incremental weak format of the governing equations (presented in Appendix B) are set up and solved in a monolithic fashion, i.e. without any staggering between the different physical mechanisms. The polynomial order of the triangular Lagrange elements (Fig. 3a) used for the various primary fields in the fibres and the SBE domains are presented in Appendix D.1. It should be noted that the dependent variables in the respective domains (fibre/SBE) are defined independently. Hence, no interface elements are used at the interior boundaries.
To illustrate the utilized convention for charge/discharge conditions, we provide a schematic illustration of the studied CF-SBE electrode half-cell in Fig. 3c. The corresponding electric potential and applied mass specific current versus time are shown in Fig. 3d. We then assume that the fibres are lithiated during discharge, i.e. Li-ions move from the Li-metal to the fibres. On the other hand, the fibres are delithiated during the charge phase, i.e. Li-ions move from the fibres to the Li-metal. Only galvanostatic conditions are considered in this study.
With respect to the mechanical loading condition in the 3 -direction, we considered the case of generalized plane stress ( 33 ( 1 , 2 , ) =̄3 3 ( )), defined by the condition: Here we note that defines a surface in 2D and that an extra condition is needed to computē3 3 ( ) as part of the FE-problem. With respect to the electrochemical loading conditions we apply a constant (dis)charge current, defined as where − is the total applied current per unit length, f is the applied mass-specific current and f is the fibre mass of the unit cell.

Parameter values
The parameter values used in the analysis are listed in Table D.2 (in Appendix D). The particular choice of relevant values is discussed in Carlstedt et al. ( , 2022b. The reference/standard value of 0 Li for Li in the fibres is based on measurements by Kjell et al. (2013), while 0 Li = 0 in the electrolyte. As to the Li-concentration dependent elastic properties of the fibres, we set ⟂ ( Li ) = 0 ⟂ [1 + 1.07̃L i ] in accordance with measurements by Duan et al. (2021). Moreover, the effective values of Lamé's parameters ( , ) and ion mobilities ( ) for the SBE are based on estimated values for a porosity of = 0.4 (cf. Carlstedt et al. (2022b) and Fig. 1b). It is also noted that the total applied current ( − ) is defined based on the mass-specific current f (cf. (33)). As reference state we utilize f = 168 A kg −1 of fibres, which corresponds to the electric current needed to discharge the carbon fibre electrode half-cell in approximately 1 h (i.e. C-rate = 1), based on measurements by Jacques et al. (2013a). Further, Li,max 7 and the insertion induced expansion coefficients ( ch ⟂ , ch ∥ ) of the fibre are based on measurements for this particular current, cf. Jacques et al. (2013a). For simplicity, these values are not altered for different values of f . Finally, we set the saturation (maximum) concentrations in the electrolyte to ,sat = 3 mol kg −1 (Nyman et al., 2008).
The interface mobility is set as:̄= 0 ∕[ 0 ]. The exchange current density ( 0 ) is assumed constant for simplicity and is based on measurements by Kjell et al. (2013). Further, the interface permittivity is chosen as = ∕ , where is the assumed thickness of the electric double layer. The relative permittivity is set to = 10, cf. Fontanella and Wintersgill (1988) and Ganser et al. (2019), and the thickness of the electric double layer is set to 0.5 nm, cf. Ganser et al. (2019) and Braun et al. (2015).
The thermal expansion coefficients of the fibres ( th ⟂ , th ∥ ) and SBE ( th ), are based on data from Bowles and Tompkins (1989). Moreover, the thermal conductivity ( ∥ , ⟂ ), volume specific heat capacity ( v ) and electronic conductivity ( f ) of the carbon fibres are based on measurements from the supplier (Torayca, 2018). The thermal conductivity ( ) and volume specific heat capacity ( e v ) of the SBE are based on data from Zantout and Zhupanska (2010). It should be noted that the thermal properties of the SBE are based on conventional D. Carlstedt et al. polymer systems for carbon fibre reinforced polymers (due to lack of experimental data for the developed SBE systems (Ihrner et al., 2017;Schneider et al., 2019)). Further, the assumed fibre-SBE interface heat transfer coefficient ( th ) is based on data from Macedo and Ferreira (2003) for polymer-based carbon fibre composites. Moreover, the fibre length is allowed to vary between 0.1-0.2 m in the parametric study performed in Section 5.4. For the remaining studies we set to = 0.1 m. The heat exchange coefficient with the surroundings ( th ext ) and the volume specific heat capacity of the SBE ( e v ) are allowed to vary between 1-10 W m −2 K −1 and 10 5 -10 7 J m −3 K −1 , respectively, in the parametric study performed in Section 5.5. For the remaining studies (i.e. excluding the parametric study), we set th ext = 1 W m −2 K −1 (cf. Gigliotti et al. (2011) and Xu et al. (2022)) and e v = 2 ⋅ 10 6 J m −3 K −1 (cf. Zantout and Zhupanska (2010)).
Finally, we note that as initial state for the simulations we set Li,0 = 0.01 in the fibres,̃, 0 = 1 in the SBE and 0 = 293.15 K.

A complete discharge/charge cycle under galvanostatic control
To demonstrate that the developed framework can predict the coupled thermo-electro-chemo-mechanical behaviour, we perform a complete discharge/charge cycle under galvanostatic control with a current of f m = 168 A kg −1 of fibres. Results from this (dis)charge process are presented in Fig. 4. The applied current corresponds to a (dis)charge time of approximately 1 h, cf. Jacques et al. (2013a). After the discharge and charge process respectively, the current is set to zero for 500 s to allow the cell to rest. During discharge, the Li-ions move from the Li-metal to the fibres via the electrolyte and insert into the fibres while the electrons move via the external circuit, see Fig. 4a. This alters the electric potential ( − ) as well as the Li-concentration in the electrolyte (̃e Li ) and the fibres (̃f Li ), as shown in Fig. 4b.
The Li-ions enter the fibres via the fibre/electrolyte interface and cause the fibres to expand, generating internal stresses and strains (Fig. 4c). Due to the fact that the free expansion of the fibres is hindered by the surrounding SBE, in addition to the assumed mechanical boundary conditions, this results in a continuous increase of the compressive stress state in the fibre during discharge (i.e. lithiation). In the SBE on the other hand, the stress state depend significantly on the assumed fibre arrangement and boundary conditions, and are both of tensile and compressive nature. During charge, the mass and charge transport processes are reversed (compared with the discharge process). This causes the stresses to decrease (based on the assumption of zero stress/strain at̃f Li = 0) as the Li-concentration in the fibres reduces.
The transport of charged species and electrochemical reactions generate (or absorb) heat, which alters the temperature (Fig. 4d). Heat is generated due to transport of charged species during both the discharge and charge processes. How the average temperature evolves within the half-cell is presented in Fig. 4d. It is noted that the temperature increases initially before it reaches a ''plateau''. Further, the temperature distribution in the through thickness direction ( 2 -direction) is found to be negligible throughout the complete (dis)charge process. This is considered reasonable given the dimensions of the studied geometry (thickness of electrode 25 μm). Moreover, it is noted that the temperature inside the material quickly returns to the temperature of the surroundings when the current is set to zero (Fig. 4d). Finally, we observe that the average temperature during charge is slightly higher as compared with the case of discharge.  Finally, we quantify the error associated with the simplification made: assuming a constant electric potential in the fibres while approximating a heat generation in the fibres due the electric current (cf. Section 3.1.3). For the given applied current the electric field component f is estimated (using Eq. (9)) to be approximately 0.26 V per metre. Hence, over the fibre length = 0.1 m this corresponds to a voltage drop of 0.026 V. This can be compared with the electric potential presented in Fig. 4b, which is ranging from approximately 0.25 to 0.55 V. It is evident that this contribution is small in comparison with the approximated electric potential. Further, in the case of a so-called full-cell, i.e. when the Li-metal is replaced with a positive electrode (see e.g. Asp et al. (2021)), this effect will be even smaller as the operating voltage is further increased (roughly 2 -3.5 V). However, we note that this effect should be considered to accurately estimate the power delivered by the battery, or to quantify the energy efficiency during operation.

Effect of (dis)charge rate/applied electric current
The influence of (dis)charge rate, i.e. applied electric current, under galvanostatic discharge is studied. The computational results for a 200 s discharge pulse at currents: f m = 84, 168, 504 A kg −1 of fibres (which corresponds to an assumed discharge time of 2, 1 and 0.33 h, respectively) are presented in Fig. 5. The initial Li-concentration in the fibres is set tõL i = 0.01. In Fig. 5b, the electric potential and applied current density versus time are shown for the three cases. It is evident that the electric potential is highly influence by the current density, as expected. Further, in Fig. 5c the temperature evaluation for the studied discharge currents are presented. Since the generated heat associated with the electronic current in the fibres depends on the applied current f m (cf. Eq. (9)), the temperature at which an equilibrium between heat generated and heat exchange is reached will increase for larger currents. For example, at the simulated discharge current f m = 504 A kg −1 (discharge time 0.33 h, or C-rate = 3) the temperature increase is approximately 9 degrees and stabilizes after approximately 200 s. In Fig. 5d, the normalized Li-concentrations and temperature at time 1 = 200 s for f m = 504 A kg −1 of fibres are presented. The temperature distribution inside the material is found to be negligible while the variation in concentrations is significant both in the SBE and fibres. The large variations in concentrations is due to the transport properties of the constituents and the applied current.

Coupling terms associated with the thermal interaction
With respect to coupling terms linked to the thermal interaction with the electrochemical and mechanical processes, we identify the following relevant couplings: (i) The electric potential − is correlated with the chemical potential of Li in the fibres, which depends on the mechanical strain, Li-concentration and temperature: Li ( , Li , ) in f and ( , ) in e , cf. (6b) and (17c)-(17d); (ii) The mechanical stress field depends on the mechanical strain, Li-concentration and temperature, i.e. ( , Li , ) in f and ( , ) in e , cf. Eqs. (6a) and (17a). To study these coupling terms we repeat the galvanostatic discharge pulse at f m = 168 and 504 A kg −1 (for 200 s) for the following two cases: • Case 1: Temperature dependence of the chemical potentials ( ) included (i.e. problem is solved as described in Section 3) • Case 2: Temperature dependence of the chemical potentials neglected. This is achieved by defining the chemical potentials (Eqs. (6b), (17c) and (17d)) in terms of a reference temperature corresponding to the starting temperature of the simulation (i.e. 293.15 K) In Fig. 6a-b, the electric potential ( − ) and the average temperature in the half-cell, versus time, are presented for the two cases presented above. It is clear that temperature dependence of the chemical potentials only has minor effect on both electric potential and temperature evolution. The effect becomes more pronounced when the temperature change increases due to the increased electric loading (Fig. 6b), as expected. It should be noted that idealized definitions of the chemical potentials are used.
In Fig. 6c-d, the in-plane strain components in the 1 -direction (cf. Eqs. (3a), (14) and (3b)) at time 1 = 200 s are presented for Case 1 and the two studied electric loading scenarios. The thermal strains ( th 11 ) in the fibres are found to be orders of magnitude lower compared with the other strain components. However, the magnitude of the thermal strains in the SBE are found to be comparable with the insertion induced strains, in particular for large temperature increase. Hence, in this case of large temperature increase (or decrease), e.g. due to large applied current and/or large variations in the surrounding temperature, the thermal strains in the SBE may have noticeable effect on the resulting strain level. Finally, we emphasize that the stiffness D. Carlstedt et al. of conventional polymers, as well as the ion conductivity of the electrolyte, are known to be highly temperature dependent (Gigliotti et al., 2011;Ihrner et al., 2017). However, these dependencies are currently not accounted for.

Significance of the individual terms in the thermal energy balance
To assess the influence of the individual contributions to the thermal energy balance, we identify and denote the respective terms associated with heat generation and/or absorption in Eqs. (1c), (11), (26) and (12e) as follows: It should be noted that 2 and 3 are relevant only in the SBE, whereas 4 relates to the heat generated at the fibre/electrolyte interface and 5 to the electron transport along the fibres. The results are extracted by running six simulations: one simulation where all terms are included ( ∑ ) and five simulations where each individual term is accounted for ( ) while the remaining terms are set to zero.
In Fig. 7, the individual contributions are shown for the case of galvanostatic discharge during 200 s at f m = 168 and 504 A kg −1 (and fibre length = 0.1 m and 0.2 m). It appears the heat generation associated with 4 is the dominating contribution to the thermal energy balance for low current (Fig. 7a). Further, the resistive (Joule) heating of the fibres ( 5 ) provides a relevant contribution, in particular for large D. Carlstedt et al.  (Fig. 7c-d). The remaining contributions are small compared with the dominating terms (cf. Fig. 7b). Further, we note that the contribution associated with 5 is proportional to 2 , where is the fibre length, cf. Eq. (11). The computed temperature evolution is confirmed qualitatively by results of Schutzeichel et al. (2021), who considered a discharge current corresponding to approximately C-rate = 12 (when setting the fibre length = 1 m) and th ext = 5 W m −2 K −1 . It should be noted that in the case of f m = 504 A kg −1 and = 0.1 m the electric field component f is approximately 0.78 V per metre, which corresponds to a drop in electric potential over the fibre length of 0.078 V. This can be compared with the operating voltage of the half-cell being around 0.2 V. Further, in the case of f m = 168 A kg −1 and = 0.2 m the corresponding drop is approximately 0.1 V while the operating voltage is around 0.4 V.

Parametric studies: th ext and v for the SBE
To assess the influence of the boundary condition associated with the convective heat flux (cf. Eq. (31b)), we repeated the galvanostatic discharge process with duration 200 s at f m = 168 A kg −1 of fibres for different heat transfer coefficients th ext (Fig. 8a). The resulting temperature evolutions are presented in Fig. 8b. It is evident that the temperature increase is highly influenced by the heat transfer coefficient. Further, to illustrate the influence of material properties on the temperature evolution, we re-ran the galvanostatic discharge pulse with different values of the volume specific heat capacity of the SBE ( e v ) for th ext = 1 W m −2 K −1 . The resulting temperature evolutions are presented in Fig. 8c. As expected, the temperature increase within the material is faster for lower values of v . It should be noted that experimental data for the developed SBE systems is largely lacking.

Conclusions and outlook to future work
In this work, a fully coupled thermo-electro-chemo-mechanical computational modelling framework for carbon fibre based structural batteries is presented. We thereby extend a previously developed computational framework to account for thermal interaction effects, which means that the energy balance now becomes an additional governing equation. We have limited the numerical studies to galvanostatic conditions.
The numerical results reveal that the dominating source for heat generation during galvanostatic cycling is associated the discontinuities in the electrical and chemical potentials at the fibre/electrolyte interface. For large temperature variations, e.g. due to large (dis)charge current or variations in external temperature, the magnitude of the thermal strains in the SBE are found to be similar to the insertion induced strains. Moreover, it is shown that the temperature change during electrochemical cycling is significantly influenced by the applied current, thermal properties of the constituents and heat exchange with the surroundings.
As to future work, we note that reliable material data for the studied material system (in particular linked to thermal dependencies) are lacking. For example, the effective thermal conductivity of the SBE is (to the authors knowledge) not available in the literature. Moreover, how the stiffness of the SBE depends on the temperature is unknown; hence, this dependence is currently not accounted for. In addition, the effective ion conductivity of the SBE is known to be highly influenced by temperature (Ihrner et al., 2017). Adding such material dependencies is part of future work. Moreover, realistic material geometries are needed when proceeding with more detailed studies, e.g. to study the mechanical damage tolerance or longevity of the material. Further, a more rigorous approach for dealing with the three dimensional problem is required in order to e.g. estimate the power delivered by the battery or quantify the energy efficiency during operation. In particular, such an approach must account for the different length scales (thickness vs. length/width) involved, e.g. by combining the micro-scale (cross-section) problem with the macro-scale structural problem. Such extension is planned for future work.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Chalmers Centre for Computational Science and Engineering (C3SE) provided by the Swedish National Infrastructure for Computing (SNIC).

A.1. Balance of momentum, charge and mass
We consider a (part of a) body that occupies the (arbitrarily chosen) domain  ⊂ R 3 with boundary  ⊂ R 2 . The balance equations of momentum, charge (Gauss' law), and mass read where is the (symmetric) stress tensor, is the electric flux density (dielectric displacement) vector and is the bulk density. Further, for any given species , is the ion concentration, is the ion flux vector, whereas ′ = is the specific charge ( is the valence number and is Faraday's constant). To simplify matters, we ignored any volume-specific reaction term in (A.1c).

A.2. Balance of energy
The body may exchange energy of different nature (electrical, chemical, thermal and mechanical) with the environment.
The global format of the energy balance reads where (the parametrization of which is discussed later) is the volumespecific internal energy. Referring to Fig. A.9, we define the power supply to the body of mechanical, thermal, electrical and chemical nature as follows: where is the displacement field, is the volume-specific body force, is the volume-specific rate of heat generation source, and is the heat flux. Further, Max ∶= × is the so-called Maxwell current, where is the magnetizing field, [ ] =A m −1 and is the electric field. Moreover, is the chemical potential for any given species, [ ] =J mol −1 . Introducing the identity = ∑ ′ , transforming surface to volume integrals and using the appropriate balance equations, we obtain where is the strain tensor. In order to derive (A.4a), we used the momentum balance (A.1a); to derive (A.4c), we used Ampere's law and the expression for the current density ; while we used the continuity Eq. (A.1c) to derive (A.4d). We may thus localize (A.2) as where ′ ∶= + ′ is the chemical potential including the effect of the electric potential (often referred to as the electrochemical potential, cf. Newman and Thomas-Alyea (2004) and Salvadori et al. (2015)) associated with the :th species.

A.3. Entropy inequality -constitutive equations
The standard global format of the entropy inequality reads where is the volume-specific internal entropy and  is the external entropy supply In (A.7), is the absolute temperature. Eliminating the expression − ⋅ upon using the localized energy balance (A.5), we may obtain the localized form of (A.6) as follows: Separating terms conveniently, we may express this inequality as where the different terms are defined as A stronger version of the entropy inequality is -Clausius-Planck Inequality (CPI):  ener ≥ 0 -Fourier Inequality (FI):  heat ≥ 0 -Fick-Ohm Inequality 8 :  el−chem ≥ 0 We shall consider these inequalities in turn:

A.3.1. Clausisus-Planck inequality
We first note that the natural parametrization of is ( , , , ), 9 where the components of are the individual concentrations . The next step is to introduce some kind of free energy by replacing with via the appropriate Legendre transformation. We shall also replace with . However, as to the role of as the independent thermodynamic variable, the choice is less obvious. The two main possibilities (which both have advantages and disadvantages) are -is kept as the independent variable (which is the option chosen here) -is replaced by as the independent variable A volume-specific free energy density ( , , , ) is thus obtained via the double Legendre transformation Upon evaluating the min in (A.11), we obtain the corresponding stationarity conditions − = 0, and − = , from which it is assumed that it is possible to solve for , in terms of , , , . Upon inserting into (A.10a), we obtain the equivalent inequality which can be rewritten as We thus obtain the energetic constitutive relations

A.3.2. Fourier inequality
It follows that  heat ≥ 0, if the following (classical) sufficient condition is chosen: where is the symmetric and positive (semi)definite thermal conductivity tensor, [ ] = W m −1 K −1 that may represent anisotropy. In the simplest case is state-independent.

A.3.3. Fick-Ohm inequality
It follows that  el−chem ≥ 0 if the following sufficient condition is chosen: where are symmetric and positive (semi)definite electrochemical mobility-diffusion tensors, [ ] = mol 2 m −1 s −1 J −1 that may represent anisotropy. In the simplest case is state-independent. Moreover, are defined as ∶= ′ . The mass flux is thus composed of two physical mechanisms: Diffusion (driven by chemical potential gradients) and migration (driven by the electric field). As will be shown later, the chemical potential consists of a combination of concentration and (mean) stress.

A.4. The energy equation revisited
In standard fashion, we rephrase the energy balance in (A.5) in terms of the free energy density (rather than the internal energy ) using the Legendre transformation (A.11). We then obtain the equation It is possible to expand to obtain the linearized energy equation where we introduced the quantities Remark 5. If electronic current e − is accounted for, in addition to the ionic current, then (A.17) is generalized as Upon introducing the identities ′ = + ′ and = ∑ ′ + e − , we arrive at the alternative formulation which is the generic format that holds whether electronic current is accounted for or not. □ 9 No dissipation is associated with deformation.
As to the further specification, we shall (in this paper) assume the following additive decomposition of the free energy density : ( , , , ) = me-th-ch ( , , ) + el ( ) + ∑ ch ( , ) + th ( ), (A.22) representing various couplings of mechanical, electrical, chemical and thermal contributions. These terms are defined explicitly in the main text. We note that, as a consequence of the specific parametrization (A.22), = .

Appendix B. Time-incremental weak format
The time-incremental weak format of the galvanostatic problem is derived as follows. First, the time intervals = ( −1 , ), whose length is = − −1 , are introduced. Next we employ the Backward Euler method for time integration; however, we deviate from the fully implicit rule by replacing the constitutive mobility tensor ( ) by −1 ∶= ( −1 ) for = Li,X, which infers forward differencing. Hence, we evaluate ∶= at = as where  Li = Li and  X = − X . The relevant solution (and test) spaces for solutions at the updated time are defined as: For the galvanostatic problem, the potential value + ( ) is prescribed along the boundary + , whereas the total current − ( ) from/to the negative collector (fibre) domains f ∶= ∪ f, is assumed to be a known function. The entire problem of solving for the updated fields at = can now be posed as follows: Find ∈Û, ∈F, Li ∈M Li , X ∈M X , Li ∈ L 2 ( f ∪ e ), X ∈ L 2 ( e ), ∈T, and − ∈ R, that solve the set of equations (B.3a)-(B.3h) given in Box I, where the pertinent constitutive relations are presented in Sections Section 3.1.2 ( f ), Section 3.2.2 ( e ), Section 3.3 ( fe ) and Section 3.4 ( ext ∪ + ).

Appendix C. Flow of charged species during electro-chemical cycling
The flow of charged species within the studied battery cell during operation is schematically illustrated in Figs. C.10a-b. Further, we assume a linear distribution of the electronic current along the fibres as illustrated in Fig. C.10c.

Appendix D. Symbols and parameters
Symbols and parameters used in the analysis presented in this paper are listed in Tables D.1 and D.2, respectively.

D.1. FE-approximations
The polynomial order of the triangular Lagrange elements used for the various primary fields in the fibres and the SBE domains are as follows: cubic (fibre, SBE), quadratic (SBE), Li quadratic (fibre, SBE), X quadratic (SBE), Li quadratic (fibre, SBE), X quadratic (SBE),    (2018) quadratic (fibre, SBE). The selected polynomial orders (for a given triangulation) are based on a convergence study, which shows that the results are reliable and not flawed by discretization errors.