Electro-chemo-mechanically coupled computational modelling of structural batteries

Structural batteries are multifunctional composites that combine load-bearing capacity with electro-chemical energy storage capability. The laminated architecture is considered in this paper, whereby restriction is made to a so called half-cell in order to focus on the main characteristics and provide a computational tool for future parameter studies. A thermodynamically consistent modelling approach is exploited for the relevant electro-chemo-mechanical system. We consider effects of lithium insertion in the carbon fibres, leading to insertion strains, while assuming transverse isotropy. Further, stress-assisted ionic transport is accounted for in addition to standard diffusion and migration. The relevant space-variational problems that result from time discretisation are established and evaluated in some detail. The proposed model framework is applied to a generic/idealized material representation to demonstrate its functionality and the importance of accounting for the electro-chemo-mechanical coupling effects. As a proof of concept, the numerical studies reveal that it is vital to account for two-way coupling in order to predict the multifunctional (i.e. combined electro-chemo-mechanical) performance of structural batteries.


Introduction
In recent years, there has been a growing interest in multifunctional materials which can enable pathways for energy efficient and sustainable transportation [1][2][3][4][5][6][7][8][9][10]. One class of such materials is the structural battery, which has the ability to simultaneously carry mechanical load and store electro-chemical energy. By combining these functionalities, the structural battery offers significant mass and volume savings for future electric vehicles and devices [11][12][13][14][15].
Several possible microstructural design principles can be envisioned. The laminated structural battery architecture was first proposed by Wetzel et al [13] and later demonstrated by Ekstedt et al [16] and Carlson [17]. The individual laminae correspond to different components in a battery cell (electrodes, separator, etc). When these laminae are stacked into a laminate, the resulting material provides mechanical properties similar to those of conventional fibre reinforced polymer composites that are commonly used in structural applications (see e.g. work by Johannisson et al [18]).
The internal structure of a conventional lithium ion battery and the laminated structural battery are shown schematically in figure 1 for comparison. The negative and positive electrodes of the laminated structural battery consist of structural battery electrolyte (SBE) [19,20] that is reinforced with carbon fibres. The SBE is a bi-continuous polymer network which contains liquid electrolyte. Hence, the SBE has the ability to both transfer mechanical load and transport ions. Since previous studies [21,22] have shown good specific mechanical and electrical properties of carbon fibres, they are well suited for multifunctional applications. In the positive electrode the carbon fibres are coated with lithium metal oxide or olivine based particles, e.g. LiFePO 4 , binder and conductive additives (see for example the work by Hagberg et al [23]). The two electrodes are separated by a porous polymer separator made from, e.g. a thin layer of SBE. The fibres in the negative electrode and the particles in the coating in the positive electrode act as active electrode materials (hosts for the lithium) in the structural battery cell.
In a conventional lithium ion battery, the electrodes consist of electrode particles and conductive additives (e.g. graphite particles and carbon black) adhered to a metal foil (current collector) using a polymer binder. The two electrodes are separated by a porous separator. This porous structure is soaked in a liquid electrolyte to allow for ion transport. Hence, the differences of the laminated structural battery cell, as compared to the conventional battery cell, are threefold: (i) The active electrode materials in the negative electrode are fibres instead of particles, (ii) SBE is used instead of liquid electrolyte and (iii) the fibres (in both electrodes) function as current collectors. The latter means that the upper and lower edges of the structural battery electrodes (see figure 1) do not need to be fully covered by metal foil like for the conventional battery cell. Furthermore, the ion conductivity of the SBE has been reported to be in the order of 10 −4 S cm −1 [19,20], whereas the conductivity of conventional liquid electrolytes is significantly higher, can be in the order of 10 −2 S cm −1 [19,24] at ambient temperature. Moreover, the electrical conductivity of carbon fibres [22] is about three orders of magnitude lower than for copper and aluminium, which are commonly used for current collectors in ordinary lithium ion batteries. This means that the geometric and topological characteristics of the cell, cf thickness of the different layers, volume fractions/fibre packaging, transportation properties of the constituents, etc are expected to have significant effects on the mechanical and electrical performance. There are currently no models available to assess these effects in structural batteries.
As to the theoretical prediction of the multifunctional properties and performance of the laminated structural battery, multiphysics modelling is needed. Apart from the seminal contributions by Newman and co-workers, e.g. [25][26][27][28], in the context of electro-chemo-mechanical modelling of conventional batteries with liquid and solid-state electrolytes, we note contributions by Purkayastha and McMeeking [29], Grazioli et al [30][31][32], Bucci et al [33], Wu and Lu [34], Ganser et al [35,36], Xu et al [37], Wan and Ciucci [38], Bower et al [39] and Wu and Lu [40], to mention a few. Moreover, Salvadori et al [41][42][43][44] have developed multi-scale and computational homogenization approaches for modelling conventional Li-ion battery cells. Examples of review articles are Xu and Zhao [45] and Zhao et al [46]. Indeed, one may favourably apply the very same conceptual model framework to structural batteries while keeping in mind that the main differences are that carbon fibres are used as the active material in the negative electrode and as current collectors in both electrodes.
In previous work by the authors, [47][48][49][50], implications on the mechanical performance from electro-chemical cycling have been studied. These studies show that the distribution of the lithium inside the active electrode materials and heat generation due to electro-chemical cycling significantly affect the effective elastic properties of the composite and its internal stress state. In these studies, one-way coupling between the electro-chemical and mechanical response was assumed (i.e. the electro-chemical analysis was used to provide input for the mechanical analysis) to simplify the analyses, and simplified geometries (representative for the so-called 3D-battery architecture, see Leijonmarck et al [51]) were used. In the context of theoretical modelling of SBEs, Tu et al [52] have studied their bifunctional performance while assuming linear constitutive relations for stiffness and ionic conductivity. Moreover, Dionisi et al [53] and Johannisson et al [54] considered a strongly simplified geometry and a simplified electro-chemical process in order to adopt an analytical model for predicting deformations and stresses in laminated structural batteries and an electrochemical actuator, respectively, due to cyclic volume change of the active materials. In summary, a computational framework to study the fully coupled electro-chemo-mechanical response of the laminated structural battery is lacking.
In this paper, we take a further step towards solving the complete multiphysics problem for the laminated structural battery cell by adopting a thermodynamically consistent theoretical framework. In particular, this means that the constitutive relations for (the energetic parts of) mechanical stresses, electrochemical potentials and the electric flux density are derivable from a volume-specific free energy density for the bulk materials. 1 Most importantly, the relevant electric, chemical and mechanical fields are resolved for a realistic microstructural design and the appropriate interface and boundary conditions. Moreover, the highly anisotropic behaviour of the fibres in the longitudinal and radial directions (transverse isotropy) is modelled. Employing the general computational framework for analysis of the laminated structural battery architecture, our main objective is to demonstrate the importance of the two-way coupling between the electro-chemical and the mechanical response, on the one hand, and on the other hand, as a proof of concept.
The paper is organized as follows: In section 2 we present the conceptual microstructural features of the laminated structural battery cell. In section 3, we restrict the analysis to the simplified problem of the negative half-cell (to obtain a manageable problem). Further, we present the governing equations for the individual domains of interest (fibre, electrolyte) as well as for the fibre-electrolyte interface and the external boundaries including the positive and negative collectors. The complete formulation of the potentiostatic and galvanostatic problems are presented in section 4 in the context of the time-incremental weak format. In section 5 the numerical implementation in the commercial finite element (FE) software COMSOL Multiphysics ® is described. Numerical results, as a proof of concept, are presented in section 6. Finally, concluding remarks and an outlook to future research are presented in section 7.

Microstructural design
The considered typical microstructural design of a structural battery composite, more specifically a so called laminated battery design, is shown in figure 2. The idea of the laminated design is to stack the laminae in such a fashion that the resulting mechanical properties become similar to those of conventional fibre reinforced polymer composites as well as providing an efficient battery function. The laminated micro-battery components are identified as follows: • The negative electrode (upper lamina in figure 2(a)): carbon fibres, with the ability to (de)insert (neutral) Li-ions, are embedded in a porous matrix (SBE). The typical radius of the fibres is 2.5 µm and the distance between the centre position of the fibres depends on the volume fraction of fibres. For a fibre volume fraction of 0.3 (assuming square fibre packing arrangement) this distance is 8.1 µm. • The positive electrode (lower lamina in figure 2(a)): carbon fibres, coated with a mixture containing Limetal-oxide or olivine based particles (such as LiFePO 4 ) and conductive additives (such as carbon black), are embedded in SBE. The typical thickness of the coating varies from around 1 µm to several micrometres (Hagberg et al [23]). • Separator/SBE (middle lamina in figure 2(a)): The separator is assumed to be made from SBE. The existence of a separator lamina assures that the active electrode materials in the two electrodes do not come in contact.
When the battery is charged, the current is brought about by Li-ions that are transported (migrate and diffuse) from the positive electrode (coating) through the separator and enters the negative electrode (fibre). Inside the fibres it is assumed that the Li-ion is neutralized and the current is conducted by electronic activity (whereby the current flows in the opposite direction to that of the electrons). When the battery is discharged the flow direction of ions and electrons is reversed.

Simplified architecture: The negative half-cell
It is possible to simplify the theoretical analysis and experimental investigation of the laminated structural battery cell by using the similarity of the two electrodes in the laminated architecture (both the negative and positive electrodes consist of fibres embedded in a SBE-based matrix material, as discussed above). It should however be noted that the main difference between the electrodes is that in the negative electrode the fibres act as the host for the lithium (i.e. expand/shrink) and current collectors simultaneously. In contrast, in the positive electrode the particles in the coating act as host for the lithium and the fibres only function as current collectors. Hence, in the negative electrode the carbon fibres will expand/shrink due to lithium insertion/extraction while in the positive electrode the particles in the coating will undergo such volume changes. Moreover, properties of the coating (e.g. particle size, distribution, etc) need to be incorporated in the model. Depending on the particle size (in comparison with dimensions of the fibres), the particles may need to be modelled explicitly or homogenized properties of the coating domain may be utilized. The provided framework can be used for modelling both electrodes but due to the relative maturity and the availability of experimental data, the negative electrode is selected for this study. We shall thus consider the negative half-cell in figure 3, which corresponds to the battery cell that was experimentally studied in previous works [18][19][20]. In principle, the positive electrode in the full-cell has been replaced by merely the collector of solid Li-metal. Moreover, the separator is excluded in the model (for simplicity).

Time-continuous strong format: Individual domains, interfaces and boundaries
In this section the time-continuous strong format and modelling assumptions for the individual domains, interfaces and boundaries are presented. A general model assumption is isothermal condition, i.e. the absolute temperature θ(x, t) = θ 0 (x) is only a given parameter. As to the mechanical properties, it is assumed that all material phases are linear elastic (small deformations). However, material nonlinearities are present in the electro-chemical relations. Self-weight is ignored, i.e. no body load is present. Finally, any piezoelectric effect is ignored, which means (in particular) that stresses do not depend explicitly on the electric field.

Fibre domain(s)
The following special assumptions are introduced: (i) The material properties are characterized as transversely isotropic, whereby isotropy pertains to the cross-section (Cartesian coordinates x 1 , x 2 ); (ii) The single active species is Li, which is neutral and can move into the carbon fibre material; (iii) In the external connection (connecting Ω f to Γ + ) all current is carried by electron transport. However, all fibres (which serve as the negative electrode) are connected to a collector with the same potential Φ − (t). Due to the high electronic conductivity (as compared to the ionic conductivity), the potential is assumed uniform, As a consequence, the electric field vanishes and the current is not resolved in the fibres (figure 3(d)). We remark that the value Φ − (t) is either given a priori (potentiostatic problem) or is computed as part of the problem solution (galvanostatic problem).
We summarize the governing balance equations in the strong format: where σ is the (symmetric) stress tensor, j Li is the mass flux vector for lithium and c Li is the mass concentration (molarity) of lithium. The relevant constitutive relations are: (c) Schematic architecture of the negative half-cell that was experimentally studied in previous works [18][19][20]. (d) Illustration of the considered idealized half-cell. In the external circuit, Load represents electric Loading, e.g. in terms of a given resistance. The introduced notation is defined in section 3.
is the (small) strain tensor expressed as a linear operator of the displacement field u, and µ Li is the chemical potential for Li. Moreover,c Li is the normalized mass concentration of Li w.r.t. its maximum concentration, c Li,max , thus defined asc Li = c Li c Li,max . As to the material properties, we introduced the elasticity tensor that is pertinent to transverse isotropy as follows: where I = E 1 + E 2 + E 3 is the 2nd order identity tensor (E i := e i ⊗ e i is the ith base dyad), I sym := 1 2 [I⊗I + I⊗I] is the (symmetric) 4th order identity tensor, 2 whereas 3 ] is a 4th order symmetric tensor. Lamé's first parameter is denoted L and the shear modulus 3 is denoted G. 3 The explicit representation of E in Voigt matrix notation is presented in appendix A, for clarity. Further, we introduced the lithium insertion strain ϵ ch (c Li ) and the mobility tensor M Li (c Li ) as follows: where the coefficient ζ = c −1 Li,max is introduced for dimensionality and α ch is a second order tensor containing the transversely isotropic coefficients of the insertion induced expansion/shrinkage of the fibres. Moreover, the reference value c Li,ref defines the state at which no chemical strains are present in the material. For simplicity, it is set equal to 0. Furthermore, the mobilities in the fibres in the transverse and longitudinal directions are assumed equal and are defined as: The constitutive relation for the chemical potential µ Li consists of two parts: The first part ζα ch : σ(ϵ, c Li ) represents the stress-driven diffusion. The second part represents standard concentration-driven diffusion. As to the model parameters, µ 0 Li is a reference/standard value for the pure species, R is the universal gas constant and θ 0 is a reference temperature. The activity coefficient f Li is chosen in accordance with the definition for an ideal solid solution of non-interacting particles on a lattice in Bazant [55]:

Remark 1. Due to the transverese isotropy of fibres, the contribution from stress to µ Li involves not only the mean stress; indeed, it is only in the case of isotropy [when
, see Xu et al [37].

Electrolyte domain Ω e
The following special assumptions are introduced: (i) The material properties are characterized as isotropic; (ii) The Li-ions are positively charged (cation, Li + ), whereas the companion X-ions (anion, e.g. PF − 6 ) are negatively charged. Insertion does not occur in the electrolyte; (iii) The current density is carried both by Li + and the companion X − . No electronic activity is permitted, i.e. no current due to motion of electrons i e − = 0; (iv) The electric potential φ may be discontinuous along each fibre-matrix interface Γ fe,i , i = 1, 2, . . . , N fibres . This discontinuity is modelled via a linearized Butler-Volmer type of 'electric resistance' relation. This relation, which involves the potential φ in the electrolyte (along the interface) and the value Φ − in the fibres, will serve as a Robin type boundary condition along each Γ fe,i (see below).
We summarize the governing balance equations in the strong format: where i MAX := ∂ t d + i ions is the Maxwell current density that is composed of two contributions: d is the electric flux density vector (dielectric displacement), whereas i ions = z ′ Li j Li + z ′ X j X is the current density due to motion of ions (known as Faraday's law of electrolysis). Here we introduced the specific charge z ′ α = Fz α for α = Li,X, where z α is the valence number and F is Faraday's constant. Moreover, we define the normalized mass concentration of species Li and X in the electrolyte asc Li  The relevant constitutive relations are: where we introduced the notation L α := z ′ α M α for α = Li,X. We introduced the standard isotropic elasticity tensor and Lame's parameters can be expressed in terms of Young's modulus (E) and Poisson's ratio (ν) as follows: . We also introduced the isotropic permittivity tensor E := εI, where ε = ε 0 ε r is the permittivity (i.e. the material's ability to transmit an electric field). The permittivity in vacuum is denoted ε 0 and the relative permittivity is denoted ε r . Moreover, φ is the electrical potential and µ X is the chemical potential for species X. The activity coefficients are assumed constant and equal to 1 in the electrolyte (f Li (c Li ) = f X (c X ) = 1). Moreover, the isotropic mobility tensors M Li (c Li ), M X (c X ), are defined as follows: Here, we made the simplification that there is no coupling between the diffusion of Li + and X − . Moreover, the mobilities are defined as: M Li (c Li ) = η Li c Li and M X (c X ) = η X c X , where η α is the mobility coefficient of species α for α = Li, X in the electrolyte.
Finally, combining the expressions above, we derive the constitutive relation for the Maxwell current: where we introduced the 'effective' electric conductivity

The fibre/electrolyte interface Γ fe
It is assumed that the redox reactions and load transfer occur along the entire fibre/electrolyte interface while assuming uniform properties. As to the mechanical integrity, it is assumed that the interface is perfectly bonded such that the displacement field u as well as the traction (σ n = σ · n, where n is the normal on Γ fe pointing out from the fibre domain Ω f ) are continuous across Γ fe . Next, we consider the electrochemical conditions. Firstly, due to redox reactions, the electrically charged Li-ion will pick up an electron and become neutral when passing through Γ fe from Ω e to Ω f , whereafter it becomes inserted in the carbon fibre. Next, we assume that µ Li , as well as φ, may be discontinuous across Γ fe , whereas j Li,n is continuous across Γ fe . A simple, yet very useful (as will be shown below), assumption is that j Li,n is governed constitutively by an interface mobilityM such that where we defineL := z ′ LiM . Here, we introduced the jump operator Li is assumed continuous across Γ fe and where φ e := φ(x e ) is the electric potential in the electrolyte at the fibre-electrolyte interface. Clearly, this model carries over directly to the current density flux i n and, consequently, to i MAX,n across the interface Γ fe . To this end, we first obtain where we introduced the assumption j X,n (x e ) = 0, x ∈ Γ fe , i.e. the transport of X − is blocked at the fibre-electrolyte interface. Further, we introduced the 'effective conductivity'K eff := [z ′ Li ] 2M = z ′ LiL . Next, we introduce the constitutive assumption (12) Now, noting that i Max,n = ∂ t d n + i n , we are in the position to give the constitutive relation for i Max,n as follows: Remark 2. It is of considerable interest to look into the nature of the relation for i n in equation (11). Recalling the explicit expressions for the chemical potential µ Li , as defined in equation (2c), we may rewrite [[µ Li ]] as follows: wherec e Li andc f Li are the normalized Li-concentration in the electrolyte and fibre, respectively, and U eq is the so called equilibrium potential that is completely consistent with the choice of the chemical potential µ in the fibre and the electrolyte in equations (2c) and (6d), respectively. This corresponds to the Nernst equation, see Newman and Thomas-Alyea [28]. To capture the features of the equilibrium potential of the carbon fibres in a more realistic fashion, it is possible to refine the chosen simple format of the activity coefficient f f Li (e.g. as done in Doyle et al [26]). Hence, equation (11) can be written more explicitly as When stress-driven diffusion is ignored (which is the classical situation in the electro-chemistry literature), it is possible to identify the expression in equation (16) as a linearization of the classical Butler-Volmer expression, which reads whereη is defined asη(x) := − φ e − Φ − + U eq , x ∈ Γ fe , and is the so called 'surface overpotential'. Further,S is the resistance coefficient,S = F Rθ0 . For small values ofη it is possible to linearize this expression, whereby equation (16) (without the stress-diffusion term) is obtained with the identityK eff = i 0S .
The current I − (t) corresponds to the total current in all the fibres that are connected to the negative connector (e.g. a strip of metal foil) that is assumed to be attached at the fibre ends. Clearly, this current has to be transported/conducted via electronic conduction along the fibres; however, it cannot be modelled in a 2D-setting. Its value is known (prescribed) in the case of a galvanostatic problem, whereas it is a 'reaction' in the case of a potentiostatic problem, when Φ − is prescribed. In the latter case the corresponding test function is set to zero.
In either case we have the conditionˆΓ which can be rephrased in the weak form as part of the galvanostatic problem.

Exterior boundaries Γ ext ∪ Γ +
Conditions at the exterior boundaries are imposed (typically) as follows: Mechanical conditions: The motivation for the traction-free condition on Γ ext,1 ∪ Γ + in equation (19b) is that the studied (part of the) lamina will in practice constitute a layered plate structure, whereby the assumption about small magnitude of the normal stress across the plate thickness is well taken.
Chemical conditions: The assumed chemical conditions are motivated by the fact that the mass flux occurs between the Li-metal and fibres (i.e. mainly in the x 2 -direction in figure 3(d)) and that the height of the studied unit corresponds to the thickness of the electrode lamina.
Electrical conditions: where 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 henceforth prescribed at 0 V (as a given reference potential).

System characteristics for discharging and (re)charging phases
The basic problem formulations are the potentiostatic and galvanostatic problems, whereby the potential Φ − (t) and the current I − (t), respectively, is the controlled (prescribed) quantity. A more general situation is that the battery cell is connected to an external circuit that contains an 'electric loading device' , see figure 3(d).
In such a situation neither Φ − (t), nor I − (t), is known; rather, they are both part of the problem solution.
Here we use the convention of a carbon fibre electrode vs. Li-metal half-cell, i.e. lithiation of the fibres is referred to as discharge and delithiation as charge. We are now in position to characterize the distinct phases of discharging and charging (which can be repeated in a cyclic fashion): Charging (delithiation): During the charging phase, the electric loading device is disconnected and either a potentiostatic or galvanostatic problem is solved. This phase ends when the battery is fully charged, which is the case when either Φ − (t) or I − (t) reaches a predefined threshold value. Given the selected convention, charging corresponds to delithiation of the fibres.
Discharging (lithiation): During the discharging phase, the electric loading device is connected. This phase ends when either Φ − (t) or I − (t) falls below a predefined threshold value. Given the selected convention, discharging corresponds to lithiation of the fibres.

Preliminaries
Let us introduce time intervals I n = (t n−1 , t n ), whose length is ∆t = t n − t n−1 . We employ the Backward Euler method for time integration; however, we deviate from the fully implicit rule by replacing the constitutive mobility tensor M α ( n c α ) by n−1 M α := M α ( n−1 c α ) for α = Li,X, which infers forward differencing. Hence, we evaluate j α := n j α at t = t n as As a direct consequence, we evaluate The relevant solution (and test) spaces for solutions at the updated time t n are defined as:

Potentiostatic problem: Controlling the electric potential Φ − (t)
The most straightforward situation is when the potential value Φ − (t) (in addition to Φ + (t) = 0) is a prescribed function in time within the negative collector (fibre) domains Ω f := ∪ i Ω f,i . The entire problem of solving for the updated fields at t = t n can now be posed as follows: Find u ∈Û, φ ∈F, µ Li ∈M Li , µ X ∈M X , c Li ∈ L 2 (Ω f ∪ Ω e ), and c X ∈ L 2 (Ω e ), that solve the set of equationŝ where the pertinent constitutive relations were given as follows: σ is defined in equation (2a) on Ω f and in equation (6a) on Ω e ; d is defined in equation (6b) on Ω e ; j Li is defined in equation (2b) on Ω f and in equation (6c) on Ω e ; j X is defined in equation (6e) on Ω e ; µ en Li and µ en X are the energetic constitutive expressions of equations (2c), (6d) and (6f ), which are given here for completeness: Moreover, j Li,n is defined in equation (10) on Γ fe and in equation (20a) on Γ + . Finally, i MAX is given in equation (23), whereas i Max,n is defined in equation (24) on Γ fe and Γ + , respectively. Remark 3. We employ a mixed method, since µ Li and µ X are treated as independent fields in addition to c Li and c X . This choice requires the additional constraint conditions in equations (26e) and (26f).

Galvanostatic problem: Controlling the electric current I − (t)
We consider the problem when the potential value Φ + (t) is prescribed (function in time) along the (positive) collector boundary Γ + , whereas the total current I − (t) from/to the negative collector (fibre) domains Ω f := ∪ i Ω f,i is assumed to be a known function. The entire problem of solving for the updated fields at t = t n can now be posed as follows: Find u ∈Û, φ ∈F, µ Li ∈M Li , µ X ∈M X , c Li ∈ L 2 (Ω f ∪ Ω e ), c X ∈ L 2 (Ω e ), and Φ − ∈ R, that solve the set of equationŝ Upon comparing with the formulation of the potentiostatic problem, we note that the value of Φ − is now part of the solution; hence, the additional equation (28g) is required.

FE-approximation and implementation in COMSOL Multiphysics ®
The numerical implementation is done in the commercial FE software COMSOL Multiphysics version 5.4. The Weak form PDE module is used to set-up the time-incremental weak format of the half-cell problem presented in section 4 and the built-in solver MUMPS (MUltifrontal Massively Parallel sparse direct Solver [56]) is used to solve the system of equations. The partial differential equations for the mass and charge balance in the SBE phase (electro-chemical analysis) are discretized with quadratic triangular Lagrange elements. The mass balance in the fibre domains is discretized with cubic triangular Lagrange elements. Finally, the displacement field (linked to the mechanical analysis), in both fibre and SBE domains, are discretized using quartic triangular Lagrange elements. The Fully Coupled Approach available in the COMSOL suite is used to solve the coupled problem which means that the complete system of equations is solved in a monolithic fashion, i.e. without using any staggering between the different physical mechanisms. Moreover, all boundary conditions are applied as Weak Contributions.

Model geometry and loading conditions
The geometry and element mesh for the chosen two-dimensional FE-model with generic/idealized geometry are illustrated in figure 4. The height (h ext ) and width (w ext ) are set equal and defined such that the fibre volume fraction is 30% (V f = 0.3). The fibre volume fraction is kept constant (for simplicity) and 30% is selected as a rough estimate based on cross-section images of carbon fibre-SBE electrodes (see figure 3(b)). The fibre distribution in the idealized model is selected arbitrarily. A biased mesh ( figure 4(b)) is used with mesh size ranging from 0.1 to 0.6 µm. The mesh size is determined such that the overall solution has converged. It should be noted that local phenomena, such as the deviation from electroneutrality in the immediate vicinity the fibre-SBE interface (see section 6.4), is highly affected by the mesh size. Although the FE-analysis is two-dimensional, it must be kept in mind that the stress state is three-dimensional due to the assumed stress/strain condition in the x 3 -direction (along the fibres), see specification below.
The mechanical boundary conditions were given in equations (19a) and (19b). As to the stress/strain conditions in the x 3 -direction, i.e. in the fibre direction (henceforth denoted loading conditions), several options are possible. Here, we opt for conditions that simulate those which are typical for beam (and plate) kinematics, whereby the x 3 -direction is the beam axis. We then consider the following alternatives: • Load(i): Standard plane strain, i.e. ϵ 33 (x 1 , x 2 , t) =ε 33 = 0. Postprocessing then gives the field σ 33 (x 1 , x 2 , t) and the normal force N 33 (t). • Load(ii): Generalized plane strain, i.e. ϵ 33 (x 1 , x 2 , t) =ε 33 (t) with a prescribed time-variation ofε 33 (t). In order to allow for comparing with experimental data (e.g. from Jacques et al [57]), we choose 'ramp-loading' of the formε where a is a constant associated with the magnitude of the applied strain and t 1 denotes the time at which no additional strain is applied. Even in this case postprocessing will give the field σ 33 (x 1 , x 2 , t) and the normal force N 33 (t). • Load(iii): Generalized plane stress, defined by ϵ 33 (x 1 , x 2 , t) =ε 33 (t) and the condition where we note that Ω defines a surface in 2D. This is the extra condition that is needed to computeε 33 (t) as part of the FE-problem. Clearly, postprocessing will provide the field σ 33 (x 1 , x 2 , t).

Material parameters
All parameters are homogeneous within the fibres and the SBE, respectively. The complete set of parameter values used in the analysis is listed in table A2. The mobility of Li in the fibres is estimated based on the longitudinal diffusion coefficient for sized IMS65 carbon fibres atc Li = 0.05 reported by Kjell et al [58] according to the relation η = D Rθ , where η is the mobility coefficient and D the measured diffusion coefficient, see Salvadori et al [43]. Moreover, the mobility coefficients for Li and X in the SBE are chosen equal and the values are approximated as η Li = η X = [59]), where K eff is the measured ion conductivity of the SBE reported by Ihrner et al [19], and c Li,ref = c X,ref = 1000 mol m −3 corresponds to the initial salt concentration in the SBE.
The reference/standard value of µ for Li in the fibres is based on measurements by Kjell et al [58] as µ 0 Li = 1 2 FU eq (c Li = 0.05), while µ 0 Li = 0 in the electrolyte. The specific capacity of the carbon fibres (C f ) is defined based on measurements on negative half-cells (carbon fibres in liquid electrolyte) by Jacques et al [60]. The selected capacity is reported for lithiation with a current rate corresponding to approximately one hour discharge time. The maximum Li-concentration in the fibres is defined based on the specific capacity according to: c Li,max = C f ρ f 3600/F. The coefficients of the insertion tensor α ch are based on measurements by Jacques et al [60] for the corresponding current. The mobility resistances linked to the mass flux and current flow at the fibre-electrolyte and Li-metal-electrolyte interfaces are approximated as:M = i 0 /[Rθ 0 F] (assuming small overpotential). For simplicity, we assume that the exchange current density (denoted i 0 ) is constant. We also assume that the interface resistance associated with permittivity can be expressed as E = ε/δ, where δ is the assumed thickness of the electric double layer. The relative permittivity is assumed equal to ε r = 10 based on previous work on polyether electrolytes (Fontanella and Wintersgill [61]) in accordance with Ganser et al [35]. The thickness of the electric double layer is set to 0.5 nm, see [35,62].

Electro-chemical cycling: Galvanostatic vs. potentiostatic control
To demonstrate the characteristics of the two basic controls of electro-chemical cycling, the studied structural battery half-cell is discharged (i.e. the fibres are being lithiated) under galvanostatic conditions (figure 5) and potentiostatic conditions (figure 6), while assuming the loading condition Load(iii) (generalized plane stress). The results are derived by solving the weak form of the governing equations described in section 4. Hence, two-way coupling between the electro-chemical and mechanical fields is utilized (i.e. the system is solved with full interaction). The cell potential (Φ − ) and current (I − ) during the two different discharge processes are presented in figures 5(b) and 6(b).
For the galvanostatic control (figure 5), the battery cell is discharged at a constant current I − (t) = I pre , estimated as: I pre = C f m f , where C f is the assumed specific capacity for the fibre and m f is the mass per unit length of the fibres. During the discharge process the current remains constant while the cell potential drops. A stopping criterion is defined as max x∈Ω fc Li (x) = 1, which corresponds to Φ − ≈0.25 V. At this moment it is assumed that the cell is fully discharged, i.e. that the fibres are assumed to be fully lithiated. It should be noted that the stopping condition is based onc Li (in the fibres), rather than on Φ − . This is due to the non-linear characteristics of the variation in equilibrium potential (see equation (15)) asc Li in the fibres approaches 1. The time for the battery cell to be fully discharged under the given conditions turns out to be approximately 3580 s, which is less than the nominally expected 1 h.
Compared with experimental data reported by Kjell et al [58] and Johannisson et al [18], the overall behaviour of the potential Φ − (t) during galvanostic discharge, as shown in figure 5(b), is in agreement (given the expected deviation linked to the simplified expression for the activity coefficient of the fibres). For the potentiostatic control (figure 6), the potential remains constant while the current drops. In this case the potential is predefined as Φ − (t) = Φ pre = 0.05 V. For simplicity the same stopping criterion (as for the galvanostatic case) is used. The current density at which this limit is reached is approximately 0.55 A m −2 . Under these conditions the discharge time is 450 s. The reason for the relatively short discharge time (compared with the galvanostatic discharge process) is the fact that the current density is controlled by the applied potential and the assumed variation of the equilibrium potential (see equation (15)). Hence, under the given conditions the current density is larger during the entire discharge process compared with the galvanostatic case.
The normalized Li-concentrationc Li and the strain field ϵ 11 (x 1 , x 2 , t) at the two time instances t 1 = 100 s and t 2 = 300 s are presented in figures 6(c)-(f). In comparison with the galvanostatic control, a more pronounced variation in the Li-concentration in the fibres and the SBE is observed. Moreover, in the beginning of the discharge process (i.e. for low c Li in the fibres) the current density is much larger compared with later in the process. This is problematic from the viewpoint of electric power losses, mass transport limitations, etc. The mechanical strains are found to develop similarly during galvanostatic and potentiostatic discharge processes. Note, however, that the strains develop at rates associated with the set electric loading conditions.
It should be noted that the utilized electrical interface condition (equation (13)) can be expressed as the linearized form of the classical Butler-Volmer equation (see equation (17)). The exponential term Fη Rθ0 in equation (17) is equal to 0.39 in the galvanostatic case and varies (with the current density) between 4.5 and 2.6 for the potentiostatic case. This means that the validity of the linearization (associated with the assumption of small overpotentials) of the Butler-Volmer equation [28] is violated for the large current densities present during the beginning of the potentiostatic discharge process. For the galvanostatic case, the introduced error/deviation (compared with the standard form of the Butler-Volmer equation in equation (17)) is on the third decimal (i.e. is considered negligible).

Assessment of coupling between electro-chemical and mechanical fields
The electro-chemo-mechanical coupling effects will be assessed in the proposed framework as follows: • Coup(i): One-way coupling: The electro-chemical problem solved independently and the result (concentration field) is used as input data for the solution of the mechanical problem. Achieved by setting Λ = 0, where the parameter Λ := ζα ch : σ is a measure of the coupling strength. • Coup(ii): Two-way coupling: The system is solved with full interaction. The computational results are obtained while considering a single discharge cycle, assuming the loading condition Load(iii) (generalized plane stress condition), and adopting galvanostatic control. How the electric potential Φ − varies with time is shown in figure 7 for the two cases of (de)coupling as defined above. It appears that the difference in Φ − (t) between Coup(i) and Coup(ii) increases with time due to the change in stress state as fibres expand with increasing c Li .
Due to the fact that the coefficients of α ch increase with reduced (dis)charge current (as reported by Jacques et al [60]), this coupling effect will be more pronounced for lower (dis)charge currents.

Assessment of out-of-plane loading conditions
The out-of-plane loading conditions are assessed by studying the three loading cases presented in section 5.2 (illustrated in figure 8(a)). The computational results in figures 8 are obtained while considering a single discharge cycle, assuming full electro-chemo-mechanical coupling and adopting galvanostatic control. For Load(ii), the prescribed strain is set to ε 33 = 0.01 and is applied at the initial time step. The stress field σ 33 (x 1 , x 2 , t) is shown in figures 8(d)-(f). Clearly, the axial stress field depends significantly on the chosen loading condition. As expected, the magnitude of σ 33 is much smaller for Load(iii), for which the axial load vanishes. The difference in σ 33 is translated to the value of Λ, which affects the transport properties. This fact is manifested in the shift of the transient part of Φ − (t), as shown in figures 8(b)-(c). Hence, the shift in potential depends on the assumed parameters associated with Λ. This means that the elastic properties of the constituents and mechanical loading/boundary conditions strongly affect the electro-chemical performance.  Compared with experimental data reported by Jacques et al [57], the shift in Φ − (t) due to applied mechanical strain is in the same order of magnitude. Jacques et al [57] reported a shift in electric potential for a carbon fibre half-cell (in liquid electrolyte) of +4.5 mV for a mechanical strain of 0.6% applied in the fibre direction. The estimated shift in Φ − (t) using the developed framework is approximately +5 mV for a mechanical strain of 1% as presented in figure 8(c).
It should be noted that the studied half-cell corresponds to a regular battery cell (i.e. a full-cell), but where one of the electrodes is replaced with Li-metal (as reference). Hence, the predicted results for the half-cell are expected to correlate with the effects on the complete structural battery (i.e. the full-cell illustrated in figure 1(b)). For example, the predicted shift in the electric potential (figures 8(b)-(c)) will also occur in the complete structural battery as the cell potential is simply the potential difference between the electrodes.
In all components made from structural batteries the insertion induced expansion and deformation of the cell will be constrained. Loading conditions 1 and 3 (i.e. Load(i) and Load(iii)) represent the extreme cases of constrained vs. unconstrained conditions. In order to determine which of these loading conditions that occur in service it is necessary to analyse the complete macrostructure of the laminated battery cells.

Assessment of electroneutrality
The total (free) charge in the electrolyte is where we used that the valence numbers are z Li = +1 and z X = −1. A very frequent assumption in the electrochemical literature on conventional Li-batteries, used as an a priori constraint, is that γ = 0. Such an electroneutrality assumption would thus infer that c Li = c X . In order to have an objective measure of the deviation from electroneutrality, we introduced the normalized (with respect to Li) chargẽ Figure 9(a) shows the fieldγ in the electrolyte domain at the time t = 1000 s for the loading case Load(iii) and galvanostatic discharging. The individual fields c Li and c X are presented in figure 9(b). Since it can be concluded from the computational result that |γ| ≤ 6 · 10 −4 , the (pre)assumption of electroneutrality would have been justified; however, close to the interfaces (e.g. fibre-electrolyte) the discrepancy between c Li and c X is more visible, which has also been observed in the literature, e.g. Ganser et al [35].

Conclusions and Outlook to future work
In this paper we present a thermodynamically consistent modelling approach for studying the electro-chemo-mechanical properties of structural batteries. The laminated architecture is studied, and restriction is made to a so called half-cell. While the SBE is considered isotropic on the studied geometric scale, the strongly anisotropic (transverse isotropy) character of the fibres is taken into account. The proposed modelling framework accounts for stress-assisted transport in addition to standard diffusion and migration. We demonstrate that the framework can be used to simulate the galvanostatic and potentiostatic charge/discharge conditions of structural batteries. Further, the numerical studies reveal that it is vital to account for two-way coupling between the mechanical and electro-chemical processes. In the case of generalized plane stress conditions (when the magnitude of the out-of-plane stresses is small), the coupling effects have minor influence on the electro-chemical performance. However, in the case of severely constrained deformation (such as plane strain) or applied mechanical loading, then the coupling effects become more pronounced. Hence, for structural batteries that are intended to carry mechanical load it becomes crucial to account for the coupled effects.
As to future development, it is desirable to refine some of the model assumptions. For example, a more refined expression for the activity coefficient of the fibres is desirable. In this work a simplified expression is used which results in some discrepancy between the numerical prediction and experimental data [18,58] for the time variation of the electric potential of the studied half-cell. Moreover, previous studies on conventional lithium ion batteries have shown that the thickness of the electrodes has a significant effect on the electrical performance [63][64][65]. Due to the inferior transport properties of the constituents in the structural battery [19,20,58], as compared with a conventional battery, the performance is expected to be highly affected by the electrode dimensions. Evaluation of such effects will be the subject of future work. Finally, the provided framework can be used for evaluating the performance of the complete structural battery by adding the separator phase and the positive electrode. To realize this, procedures for dealing with the electrode coating in the positive electrode need to be developed. This will also be the scope of future work.

Appendix A. Symbols and parameters
Symbols and parameters used in the analysis presented in this paper are listed in tables A1 and A2. Moreover, the explicit representation of the elasticity tensor E (for transverse isotropy) in Voigt matrix notation is defined as