Modelling of local mechanical failures in solid oxide cell stacks

Abstract Solid oxide cells can deliver highly efficient energy conversions between electricity and fuels/chemicals. A central challenge of upscaling solid oxide cells is the probability of failure of the brittle ceramic components. The failures of the ceramic components may lead to significant degradation or eventual failure of a stack. To predict mechanical failures in a stack, a full stack model is needed, together with a local assessment of stresses at the vicinity of failing regions, e.g. the contact points between the cells and interconnects. A conventional three-dimensional model requires a very fine discretization of the mesh to capture stress intensities. Computational resources needed for such a model are therefore immense, and it is highly unlikely to compute at stack scale, as well describe the evolution over time. In this work, the homogenization modelling framework for solid oxide cell stacks is extended to identify local mechanical failures. Thus, the fracturing within a local failing point is examined by using a localization approach, where stresses in the stack model are linked to the local stresses and the energy release rate at the crack tip of the relevant interface. This is done in a general manner, such that the local stresses and the energy release rate can be evaluated at every point in the stack at every instant of time without loss of computational efficiency. A 100-cell stack can be modelled in three dimensions with all coupled multiphysics in steady state within 3 min on a current workstation computer.


Introduction
Solid oxide cells (SOCs) are able to convert chemical energy to electrical energy and vice versa with high efficiency. They can be operated to convert the chemical energy of a fuel into electricity as solid oxide fuel cell (SOFC), or run reversely, converting electrical energy into chemical energy as solid oxide electrolysis cell (SOEC). Comparing to other energy conversion technologies, SOC technology has many attractive features such as high efficiency, high fuel flexibility, good scalability for commercial purposes, etc., cf. [1].
Mechanical failure is one of the key challenges for SOC technology to withstand the desired long operation times and changes in operation. Stresses induced by temperature gradients across the stack generated during operation can cause failure in the components or interfaces between the layers [2][3][4]. Fractures in the cells and sealings may lead to cross-over of reactants, resulting in hot-spots, and ultimately cause failure of an SOC stack. Observations and measurements on the interface between the electrode and interconnect indicate that the interfaces around the oxygen electrode-contact layer are the weakest link. Loss of contact between the cell and interconnect will result in a redirection of current, higher resistance to overall current transport, as well as local Ohmic heating, cf. [3,5,6]. This can in turn accelerate other degradation phenomena, such as corrosion or degradation of the * Corresponding author.
E-mail address: ximi@dtu.dk (X.-Y. Miao). electrodes [7]. The contact layer between the electrode and metallic interconnect has gained some attentions due to the low fracture energy, and the degradation effects to the SOC stacks [8][9][10]. To estimate the mechanical response in the local regions, efforts in thermo-mechanical modelling of the metallic interconnect have been reported in [11] and [12].
Experiments provide the ultimate proof of a certain SOC stack concept and operation modes, as well as performance evaluation over the designed lifetime. This is however by the expense of time and money, in particular when considering the very high number of variables in the design and operation of SOC stacks. Computational simulations are therefore a cost and time saving complementary method to estimate the overall behaviour of a full stack from the interactions of multiple physical phenomena over the entire lifetime [13,14]. Particularly, computational modelling can evaluate if there are possible function degradations or failures with relevant local quantities that cannot be measured properly in experiments.
SOC stacks have many geometric features and are subjected to multiple physical coupled processes (electrochemical, current, mass and heat transport, mechanical, etc.). Models that contain more physical processes, material, and geometric details imply higher computational cost. The state-of-the-art SOC stack models require large computational resources (computer clusters) and still spend around 40 h for steady state simulations, cf. [15] and [16]. Simulations of stack performance at various operating conditions for a particular stack design will thus be a very time-consuming task, and transient simulations over the full lifetime of a stack seem unfeasible.
To greatly reduce the computational time of stack scale modelling, a new modelling concept using the homogenization approach has been developed [16], which allows simulations of the relevant physical phenomena, i.e. electrochemical reactions, fluid flow, heat, mass and electron transfer, as well as thermo-mechanical stresses for a full stack, e.g. a 100-cell stack in three dimensions within 15 min. In the homogenization framework, the SOC stack is considered as a porous chemical reactor (Fig. 1), where the reaction rates are determined from the equivalent electrical resistances from the current running through the cells. This provides an overall response of a stack model, which is similar to that of the actual stack, but without directly including the detailed information of e.g. local stresses between the cell and interconnect. Furthermore, for presenting the modelling concept in a comprehensible manner, various analytical approximations to the effective transport properties were employed. However, because the geometrical details of the stack are implicitly described in the homogenization modelling, the local response of the sub-structures cannot be examined.
In this work, we extend the previously presented homogenization modelling framework for SOC stacks [16] with the ability to examine for local mechanical failures using fracture mechanics, exemplified by the failure in the contact layer between the cell and interconnect. This approach makes it possible to investigate for local fractures in components on full-scale stack for the first time. With such a model the stack design operating conditions can be optimized towards better stack reliability. The modelling framework thus provides a very efficient tool for fast screening of materials and design parameters in relation to the fuel uniformity to reduce the risk of local mechanical fracturing. Finally, we also explain and demonstrate how the effective transport properties can be determined using numerical modelling on a realistic interconnect geometry, rather than a simple crude non-industrially relevant geometry. This paper is organized as follows. In Section 2, homogenization model concept is briefly reviewed. In Section 3, effective properties for the homogenization treatment to the demonstrated stack are obtained by numerical modelling of a repeating volume element where a realistic interconnect geometry is considered. The main contribution of the localization approach for the identification of the energy release rate at the wedge front of the contact layer is presented in Section 4. Validation of the proposed approach is performed and discussed in Section 5. A series of case studies for an SOC stack associated with different external loadings, manifold dimensions, number of cells, design of headers, and flow configurations are conducted in Section 6. Further discussions and main conclusions of this paper are stated in Sections 7 and 8, respectively.

Model concept
The homogenization approach aims to describe complex geometric features of a composite structure, as well as physical processes within this structure, when a fully detailed model heavily challenges the computational resources. By this approach the structural response can be determined from in an efficient manner, cf. [14,15,17] and [16].
In this work, we apply the homogenization concept to model the SOC stack shown in Fig. 1(a) under operation. To do so, an equivalent homogeneous material phase with anisotropic properties is introduced to represent all components associated with different materials and geometric details integrated in this stack, as demonstrated in Fig. 1(b). The relevant equivalent properties can be obtained by applying the volume weighted method [16,18], see Eq. (1). The detailed description of the homogenization modelling for SOC stacks has been presented in our previous work [16]. Insert on cracking between the oxygen electrode and interconnect after [6]. where represents the domain of the representative volume element (RVE), is the volume of the RVE. is an arbitrary physical quantity, and̄is the volume averaged quantity. This can be done for many physics, and the equations used for each of the domains, named in Fig. 2, are shown in Table 1. The definitions of relevant parameters covered in this work are given in Nomenclature.
The main challenge is however to describe electrochemistry, which is localized to distinct parts of the representative volume. This has been described in detail in [16]. Here it is shown that the driving forces from the electrochemical conversion could be added to the charge transfer equation as an external current density, i.e. and Where is the electrical conductivity tensor, the voltage, and ext the external current density vector. is the cell open circuit voltage. ℎ RVE is the effective height of the representative volume. The area specific resistance ( ) is for the cell according to the work of Leonide et al. [19].
Since the stack is considered as a porous chemical reactor in homogenization modelling, the gases distribute over the entire active volume. The reaction rates of the reacting species are determined by the volumetric current density. The diffusion of species are thus also described in an overall manner in the entire representative volume, as well the transport of charges. The diffusion of species within the active porous electrodes are addressed using a concentration overpotential included in the expression. All the overpotentials indicate losses in the conversion between chemical and electrical energy, i.e. heat and are thus included as heat sources in the homogeneous media, see [16]. The local overpotentials can provide local quantities for e.g. describing the degradation. In this work, we provide a method for the mechanical localization by a mapping between macroscopic quantities, i.e. average stresses and a local failure criterion. The novel mapping approach thus establishes a link between the local energy release rate at the failing zone and the average stresses computed in the homogenized stack model. Importantly is that the mapping is through an analytical expression such that the computational efficiency of the homogenized model is not compromised, see Fig. 2.
Another important point is that the material parameters for this type of model can be obtained using the actual materials and components used in SOC stacks, see e.g. [9,10,20].

Numerical determination of homogenized parameters
The determination of effective physical properties for homogenization modelling can be performed with analytical or numerical approaches in terms of the structural complexity of a composite [21][22][23]. An analytical method was employed in [16], because a simple straight rectangular channel was adopted in their stack model. Here we present how the effective properties can be determined with numerical simulation for a realistic interconnect geometry, e.g. an interconnect obtained with the shaping of a thin plate. The stack dimension and interconnect channel geometry chosen to exemplify this procedure is shown in Fig. 1, and the geometric parameters of the representative volume is provided in Table 2.
The evaluation of the effective flow parameters, thermal conductivity, stiffness matrix and thermal expansion are presented in the following sections.

Homogenization for flow in the headers
The flow in the channels and headers can also be homogenized. In the headers the gas streams flow between two parallel plates (see Fig. 1(b)), and the flow resistance is thus equivalent to that of a very  wide channel, and this can be solved analytically. The flow in the channels must be solved numerically for the chosen geometry. With relevant flow velocities, the flows in both domains are laminar. Thus for the homogeneous medium analogy, the flow can be described with advantage using Darcy's law by finding the equivalent permeability as in [16], as this is far less computationally demanding than solving the Navier-Stokes equations for this region. The equivalent permeability can be deducted by equating the description of the laminar shearing flow in the channel with Darcy's law. The flow velocity distribution in the channel can be described by solving the Hagen-Poiseuille differential equation where is the fluid velocity, the pressure, and the dynamic viscosity.
The average velocity can be determined from the resulting velocity profile solving the equation above.
where * is a solution specific constant, and ( , ) is a function introduced to describe the velocity profile over the cross section. Darcy's law states Combining Eqs. (6) and (7), it can be seen that the permeability can be expressed as The effective permeability of the channel for the flow can consequently be calculated by volume averaginḡ where ch is the domain of the channel region. The out-of-plane component of the permeability is zero, and the permeability tensor components are set to equal for the in-plane ( -plane, see Fig. 1(b)). Using this approach the flow distribution in the interconnect channels can be modelled as shown in Fig. 1(b). The effective permeability of the homogeneous porous medium representing the channels is calculated with Eq. (9), and̄c hannel = 5.2 × 10 −8 m 2 . The effective permeability for modelling the flow in the headers can be obtained in a similar manner, however by using the analytical solution for the velocity variation between the two parallel plates (̄h eader = 3.9 × 10 −8 m 2 ).

Average thermal conductivity tensor
The diagonal thermal conductivity tensor components̄are obtained by in turn imposing a temperature difference on the two opposing sets of boundaries in the direction of each of the three Cartesian coordinates, see Fig. 3 as an example. Specifically, we first solve the heat flux based on the applied temperature difference using numerical simulation. Then the total flux through a boundary can be obtained by integrating the heat flux over the oriented boundary, i.e.
where represents the boundary of the representative volume, and is the total flux through the boundary. With the known temperature difference and heat flux, the effective thermal conductivity can be calculated using Eq. (11).
where is the cross-sectional surface area, and is the distance between the two opposing boundaries.
This numerical approach can be used to determine other similar effective transport properties associated with the operating processes.
The thermal properties of the materials in the stack are shown in Table 3. Mechanically and thermally the cell can to a great extend be represented by the fuel electrode support layer, as this is comparatively thicker than the remaining. Of course a weighted average could be chosen for an actual stack, but in this methodology paper we refrain from doing so.

Effective stiffness matrix
As we use a homogeneous volume to represent the three-layer representative volume composite, the linear elasticity of such an equivalent homogeneous material phase can be represented by the stiffness matrix of an orthotropic material imposed on the two boundaries of representative volume perpendicular to direction.

Table 3
Thermal properties of the stack components.
Material layers Thermal expansion coefficient Thermal conductivity   The partial sealing area allows the gas streams from the manifolds flow to the headers or vice versa, see Fig. 1(a).
wherēare the components of the effective stiffness matrix.̄,ā re the components of the volume averaged stress and strain tensor. Eq. (12) represents a system of linear equations. The effective stiffness components can be determined by imposing different loading cases (uniaxial, shear) on the RVE sub-model as illustrated in Fig. 4 [12,22,31]. To ensure a single solution for the linear equations, specific boundary constraints are imposed on each single RVE sub-model, where volume averaged stresses and strains are calculated by applying the Gauss's theorem to the corresponding boundary values. A detailed description can be found in [31]. The RVE sub-models are performed in Abaqus/CAE 2018.
The mechanical properties of the materials in the stack shown in Table 4 are used as an example.

Average thermal expansion tensor
The average thermal expansion of the representative volume is used to estimate the linear elastic response of the homogenized stack model induced by the temperature fields, the boundary conditions (e.g. boundary constraints, external loads), as well as the thermo-mechanical coupling in the homogenized stack model, see Eq.
wherēis the effective/volume averaged thermal expansion coefficient tensor, which is used to assemble the global constitutive equation over the equivalent homogeneous volume [16] =̄− 1 ∶ The integration in Eq. (18) is thus done over the representative volume with the material parameters listed in Tables 5 and 6. The detailed calculation of other effective physical properties associated with the operating processes of an SOC stack can be found in [16].

Modelling localized fractures
The homogenization modelling estimates the average behaviour of an SOC stack. The stresses solved in this model are thus the average stresses̄in the equivalent homogeneous medium. The variation of local stresses in the representative volume is not directly available from the homogenized macro-model. To identify the occurrence of cracks in the local critical regions, the local stresses in the RVE submodel must be linked to the average stresses in the homogenized stack model. The correspondence between the average stresses and local stresses is explained in Section 4.1, which could be used to determine the probability of failure in e.g. the cells, using Weibull statistics and measured strength variations. For fractures at e.g. the contact layer between the electrode and interconnect, a novel method linking the local energy release rate and average strain energy is presented in Section 4.2. The latter approach is validated for three simple cases in Section 5.

Localization of stresses
We calculate the local stresses by the use of sub-modelling of the representative volume, on which the average stresses in the macromodel are applied as the boundary conditions (i.e. tractions), see Fig. 5.  7.6 7.6 0 11.9 11.9 11.9 7.9 7.9 (air) 0 7.6 7.6 (fuel) 0 By doing so we are able to link the homogenized mechanical response to the local stress state at any point in the representative volume. However, running a sub-model for every representative volume in the stack is a computationally intensive task and should be avoided. Thus we establish a link between the local stresses and the average stresses once and for all. The local stresses, , at a single point in the representative volume are determined by building a system of linear equations, where the variables are the average stresses applied on the representative volume,̄: = ∶̄ (19) where the components of the coefficient matrix can be obtained by solving the RVE sub-model for a set of simple stress states where only one of the stress components is nonzero, i.e. uniaxial tension and pure shear stress. Once is determined for the representative volume, the local stresses in any representative volume can be determined from the present average stresses̄using Eqs. (20)-(21) based on the linear superposition principle for linear elasticity. This concept is not new. But it is important to show here, as this leads to the next novel step.
We exemplify the calculation of for a two-dimensional case. In that case Eq. (19) is simplified from six to three linear equations: The coefficients 2 for a specific point with the local stresses { , , } are determined by loading the representative volume with a uniaxial tension in the direction shown in Fig. 4(a). Becauseā nd̄are zero, Eqs.

Localization using fracture mechanics
For fractures in the contact layer, knowing the location of the critical region where cracks are most likely to occur, the above process can be refined to obtain the critical local stresses for any combination of the average stresses. However, such critical stresses are strongly influenced by the geometry at the origin of the crack. For example, stress concentration at a sharp notch is higher than at a blunt geometry, and the effective strength of the joint will vary according to this, cf. [32]. If stresses were used as a failure criterion, the statistical distribution (e.g. Weibull distribution) of the strength must be measured for specific geometries with varying geometric features/flaws. However, manufacturing samples for testing the strength of a single contact point or even a group of contact points are challenging, and imposing a robust loading on such a sample would be even more challenging. If the geometry of the interconnect, or the processing method was changed, a different variation of geometry should be expected, and the strength variation has to be remeasured.
An alternative approach is to use linear elastic fracture mechanics, i.e. using fracturing parameters such as stress intensity factor, strain energy release rate, or the -integral to characterize the crack tip properties. In this work, the strain energy release rate is used as an indicator for the identification of crack propagation, i.e. investigate if the energy release rate is larger than the critical energy release rate c . By this approach it is examined if a pre-initiated crack will grow further. This is a conservative approach as the energy release rate needed to initiate a crack is often higher than that needed to propagate an already formed crack. However, the advantage is that the measurement of fracture energy is not dependent on the geometry of the interconnect, or the local geometry at the origin of the crack. Furthermore, the measurement of fracture energy of interfaces can be done in a robust manner at both room temperature and elevated temperature [9,10,20]. For the first realization of the proposed approach, we consider a very long interconnect rip as those used in many stack designs, which can be approximated with a plane strain model as shown in Fig. 4. We hereby neglect out-of-plane variations in the direction.
The energy release rate can be obtained by evaluating the -integral in the framework of linear elastic fracture mechanics ( = , see Table 7 Three mapping factors for linking the -integral to the average strain energy for a plane strain RVE sub-model. 11 33 13 / mm 2.0 7.9 0.7 Fig. 2). This can be done by solving the RVE sub-model with a preinitiated crack at the location of interest (e.g. at the contact interface) and evaluating the -integral along a path around the crack tip, see Fig. 5. Again, solving -integral for every contact point of the stack will be a highly computational expensive task. As in Section 4.1, there is a possibility of linking the local mechanical state to the average stresses. Because the -integral can be linked to the local strain energy around the crack tip (cf. [33]), it is possible to determine the correspondence between the -integral and the average strain energy for the simple stress states (i.e. uniaxial and pure shear stress). Again this can be done by introducing a mapping factor, e.g.
= 33̄̄f or the uniaxial tension shown in Fig. 4(a). The -integral for an arbitrary stress state can thus be obtained by linearly combining the contribution of strain energy from each simple stress state. Because the total strain energy is an addition of the contributions from each simple stress state, e.g. for a two-dimensional casē=  (Table 7). These are obtained by separately computing the -integral and average strain energy for each simple stress state in the RVE sub-model, leaving as unknowns to be determined. The -integral and average strain energy are calculated using Abaqus/CAE 2018 (e.g. Fig. 5), where the average (volume averaged) strains can be extracted by integration using a python script in the post processing of Abaqus.

Validation of the localization approach (J-integral expression)
To confirm the validity of the linear superposition principle, theintegral at the wedge front between the cell and metallic interconnect for the different specific stress states listed in Table 8 are solved using Eq. (23) and directly in Abaqus.
As seen in Table 8 the results from the analytical expression Eq. (23) based on the linear superposition principle are very close to the direct numerical simulation of the arbitrary stress state. A small deviation is present, which stems from the slight distortion of the representative volume boundaries under the fundamental loading cases. This would not be present in a repeating geometry, and therefore not a completely correct representation of the repeating geometry. Maybe these distortions can be mitigated in future work by imposing other loading conditions with variable tractions on the boundaries of the representative volume.

Results
In this section, we demonstrate the use of the homogenization stack model with the developed addition of the evaluation of the local energy release rate for a failure in the contact layer between the cell and interconnect. The proposed model is applied to explore the occurrence of local fracturing in a specific SOC stack (fuel cell mode), and how the different external loadings, number of cells, and flow configurations affect the energy release rate and thus the risk of failure.

Base case
We first simulate a base case of a 100-cell stack with manifolds to examine the effects of the multiphysics couplings. The cell model used here is from [19] and the implementation of this is further described in [16]. The dimensions of the stack chosen for demonstration are listed in Table 2, corresponding to Fig. 1(a). The geometric sizes of each component in the representative volume are also given in this table.

Stack boundary conditions
The boundary conditions used for the base case are presented in Table 9. A distributed compressive load of 1 MPa (cf. [34]), corresponding to a total load of 1019.7 kg, is applied on the top surface of the active area. The bottom surface is constrained in the vertical direction, and two boundary vertices are also constrained to prevent a rigid body motion and rotations. Due to the symmetrical characteristic of the chosen stack geometry, only half the stack is modelled, and symmetrical boundary conditions are imposed.
The boundary conditions specified for operating the stack correspond to the typical boundary conditions for SOC stacks running in fuel cell mode. The two gas streams (air and fuel) are operated in co-flow configuration and treated as ideal gases. Note that the fluid properties are variables in the model and they are a function of the gas composition, e.g. the viscosity is a function of the current gas mole fractions [35]. The flow directions are shown in Fig. 1.
Reynolds numbers are estimated based on the definition for flow in a pipe to determine the flow regime. The air flow in the manifolds is modelled as turbulent flow as the Reynolds number for the chosen flow conditions reaches 5000. This is done by a k-turbulence model using Reynolds-averaged Navier-Stokes equations. The fuel flow in the manifolds is modelled as laminar flow as the Reynolds number is below 100. The Reynolds numbers for the air flow in the active domain is 16 ∼ 18, and for the fuel flow is 0.1 ∼ 0.5, which reflect a laminar flow regime. The gas flows are evaluated as fully-developed flows. The reason is that we only investigate the stack under a steady operation, and that the flows in the channels are laminar. Thus the velocity profiles do not change along the flow direction. Therefore the convective heat transfer between the gas flows and the solid walls is modelled using Newton's laws of cooling, where the heat transfer coefficient is evaluated from Nusselt number. The overall removal of the heat is thus limited by the flow rate (heat capacity) of the gases.
The assembly temperature of 800 • C is set to be the reference temperature (zero stress condition) of the solid parts. It is a uniform temperature condition which is applied on all of the solid domains. This is a good approximation as the stresses in cells relax at reduction [36,37] and the glass sealing crystallizes at this temperature. Inlet temperatures are set to 700 • C for both air and fuel flow at the bottom of the inflow manifolds, see Fig. 1. In the current model, we consider the stack is placed alone and insulated. Radiation through the insulation material is neglected, as conduction is more dominant, cf. [35]. At the outer surfaces of the solid parts, insulation is modelled by a heat flux proportional to the temperature difference between the solid parts and the surrounding environment multiplied by the insulating properties (thermal conductivity, thickness), see Eq. (50) in [16]. The surface-to-surface radiation was comprehensively investigated in [26,38] and [14]. The effects of the radiation heat transfer on the overall temperature are usually small [26,38]. However, the radiation heat exchange inside the flow channels may affect the temperature field, cf. [38,39] and [40]. As the focus of the present work is modelling of local mechanical failures, a detailed investigation of radiation heat transfer is beyond the scope of this work. Thus we does not consider the surface-to-surface radiation as well.

Results on multiphysics modelling
For the envisaged stack design and operating conditions, different physical quantities distribute across the stack are shown in Fig. 6. Different basic expected responses from the stack model can be seen: · Significant temperature gradients induced by Joule heating and electrochemical reactions of the passing gases ( Fig. 1(b)) are observed in Fig. 6(a).
· The gases are thus heated as they flow through the stack. This is the primary factor playing into the temperature distribution in the stack, i.e. high temperatures are present at the outlet side of the gas flows.
· Due to the heat conduction in the solids and out through the insulation around the stack, small temperature gradients perpendicular to the flow direction are also formed, both horizontally and vertically.
· Species concentration shown in Fig. 6(b) indicates how the hydrogen is consumed during the fuel cell operation as expected, but the flow of hydrogen is rather non-uniform.  · The average stresses induced by the temperature gradients are shown in Fig. 6(f) as the first principal stress distribution in the active domain. The variation of the stress components which contribute to the fracturing of the contact layer between the cell and interconnect is illustrated in Fig. 7. The fracture toughness values of the commonly used contact layers are between 1 ∼ 2 J m −2 (cf. [20] and [9]), which indicates that the -integral value calculated with the present model settings is large enough to crack the contact layer. Specifically, we extract a profile of -integral value across the active area where the location of the highest -integral is covered, see Fig. 7(d). It is seen that the outermost 12 mm channel parts in both inlet and outlet side will probably delaminate in this stack, as the -integral values in these areas are higher than the maximum limit of the contact layer, i.e. 2 J m −2 .
In reality the stresses act on a representative volume of 4 mm width. A more precise effect of the very non-linear stress variation at boundaries could also be determined by imposing the stresses as traction boundary conditions of a representative volume to determine the -integral. Another, simpler approach is to assume that the average stresses over the outermost 4 mm acts on a single RVE and determine the -integral from that using Eq. (23). The average stresses acting over the representative volume are determined to 16   ; (d) current density; (e) pressure on the air side; (f) first principal stress.

Case studies
To explore the influence of stack design on the risk of fracturing of the contact layer, we present five case studies: case (A) with varying external loads on a 100-cell stack; case (B) with different manifold dimensions for redistribution of flow over the height of the stack; case (C) with different number of cells; case (D) with additional header with high flow resistance for the fuel flow; case (E) with counter-flow configuration.

Case A: variation of external loadings
In case A, we impose a series of distributed loads ranging from 0.01 to 2 MPa to the 100-cell base case stack to investigate how the external compressive loading affects the -integral.   which is the highest one among the present loading cases. The maximum reduction of the -integral to 30.4 J m −2 is reached at 2 MPa. This is because the applied compressive load counteracts the internal vertical expansion due to the temperature field in the stack, reducing the contribution of the normal tensile stresses in the direction, which is the primary contributor to the -integral.
Such high compressive loadings are however complicated to realize in reality and the interconnects may creep over time due to the high external loads, cf. [12,13] and [42]. Modelling the impact of creep in this stack modelling framework will be addressed in future works.

Case B: influence of manifold dimensions
In the base case model, we observe high stresses formed in the outermost parts of the active area due to the differences of thermal expansions across the domain induced by the non-uniformly distributed temperature, which lead to high -integral values. In this case study, we study whether these high stresses can be reduced by optimizing the thermal profile in the stack through enlarging the dimensions of the manifold for the air flow and hereby obtain a more uniform flow field and temperature distribution. The size of the manifold is increased by extending the trapezoidal interconnect in the flow direction and maintain the outermost width perpendicular to the flow direction.
It can be seen from Fig. 8(b) that the highest -integral value reduces from 35 J m −2 to 32 J m −2 when the area of the manifold is enlarged from 11 cm 2 to 16 cm 2 . The -integral value decreases monotonically with the increase of the area of the manifold. It is also worth noting that the impact of increasing the dimension of the manifold decreases with the increase of the manifold area. A further increase upon 16 cm 2 in the manifold area has no notable effect on the -integral as the flow has almost been uniformly distributed towards the corners under the current operating conditions.

Case C: variation of number of cells
In case C, the influence of the number of cells on the -integral is examined, having 20, 50, 70, 100 cells. Note that we conduct the comparison on the condition that the fuel utilization is the same in all simulated stack models. The air flow is also adjusted proportionally to the number of cells. The air flow in the manifolds for the stacks with less than 50 cells is modelled as laminar flow, as the Reynolds number is below 2300. Fig. 8(c) shows the maximum -integral value in the active area changing with the number of cells. We can see that this -integral value is smaller in the stacks with fewer cells, because lower stresses are found in these stack models. This is because more uniform temperature gradients are formed in the stack models with fewer cells due to the more uniform fuel flow distributions from the manifolds to the active area.
Comparing the -integral values calculated in case A -case C, we can conclude that the number of cells of the stack has a greater influence on the -integral value than the external loading and manifold dimension. From a mechanical perspective, a stack with less cells should be more stable.

Case D: header with high flow resistance
As revealed in case C, the uniformity of the fuel flow distribution has a significant impact on the temperature distribution across the active domain. A more non-uniform distributed flow will lead to higher temperature gradients not only along the flow direction but also along the height and width of the stack. As a consequence, higher stresses will occur.
In this case study, we thus add an additional header adjacent to the active domain to the base stack model, see Fig. 9(a) to remedy this. This new header has a thickness of 2.5 mm and has high flow resistance that ensures uniform pressure in the outer header and thus regulates the fuel flow to be more evenly distributed before it flows into the active area.
In the homogenization modelling, the high flow resistance of header is realized by lowering the permeability. The permeability of this extra header is set 15 times smaller than the original headers. We can see from Fig. 9(c)

Case E: counter-flow configuration
In case E, we examine whether flow configuration has a notable effect on the -integral. In this case we thus switch the boundary conditions for the fuel flow to redirect the flow to be opposite to the air flow.
The -integral for a counter-flow stack under the base case conditions is shown in Fig. 10. The sub-figures are plotted to highlight the critical regions with high -integral values. Particularly, a new critical region is formed within the inner part of the active domain in this model configuration as shown in Fig. 10. This is mainly due to the change in the temperature gradients over the stack under the counterflow configuration, compare Figs. 10(a) and 6(a). Besides, there are also some stress concentrations at the outermost corners of the active area, where the highest -integral (10 J m −2 ) is found.

Numerical settings and computational times
The presented homogenized stack model is realized with the commercial software COMSOL Multiphysics 5.5. We use either the segregated solution or the fully coupled approach to solve the models. The fully coupled approach requires relatively more memory and computational time to solve all the unknowns simultaneously, while segregated approach sequentially solves the unknowns in several small segregated steps according to the corresponding physics, e.g. heat transfer, mass transport, solid mechanics, etc., and thus less memory and computing time are required. Which solution procedure to use depends on the robustness of computing for a specific model. A direct solver (PARDISO) is always used for solving the linear system. Newton iterations are performed to resolve non-linearities with a relative tolerance of 10 −4 .
A mesh consisting of 73,160 elements is employed in the model of the base case. Mesh convergence is first performed to guarantee that mesh independence is present for the numerical solutions. The relative errors below 1% between stress components from the current mesh density and a finer mesh with 156,600 elements indicate a further increase of element number has a negligible influence on the results. Same mesh density along with numerical settings for the base case is applied for the stack models with varying number of cells.
All calculations are performed on a workstation with an Intel(R) Xeon(R) W-2135 CPU @ 3.70 GHz 6 Cores, and 128 GB installed RAM. For modelling a new stack configuration, the fully coupled solution procedure may be required for finding the numerical solutions due to the high non-linearities in the turbulence model equations and the convergence issue of iterative coupling using the segregated solution approach, increasing the computational time, i.e. ∼ 27 min to solve a 100-cell stack model. This is reasonable when compared to the ∼ 15 min needed to solve in previous works [16] as we here also model the turbulent flow in the manifolds and have significantly more degrees of freedom because of this as well. For the parametric studies, the previous solutions can be used as initial conditions to define better initial values, which is especially important for the turbulent flows. In this case it is possible to use the segregated solution approach, which significantly reduces the computation time down to 3 min for a 100-cell stack model per run.

Discussions
Loss of contact is a known but less investigated mechanical failure for the SOFC stacks. Some experimental studies reported this failure, such as Pan et al. [43] and Boccaccini et al. [6], who observed cracks formed at the contact spots in a post mortem analysis of an SOFC stack. Blum et al. [3,5] investigated the contact problems based on the many stack testing results obtained from their group. He investigated the cooling from sealing temperature to operation temperature, and showed that the thermal contractions in the sealing and active regions of the stack were different, which could lead to their observed loss of contact. This mismatch is accounted for in the current work,   using a reference temperature, which is different from the operating temperature.
In [44] the stress in the cell-interconnect interface was compared for co-flow and counter-flow by numerical modelling. It was seen that the stress level in co-flow was lower than in counter-flow. This relates to the more linear thermal profile in the co-flow case, which provides smaller stress levels [42]. In this work, we observe a significant influence from the local cooling at the air inlet ( Fig. 6(a)). This cooling induces a local contraction, which results in tensile stresses ( Fig. 7(b)) in all directions and thus also results in a significantintegral ( Fig. 7(d)) at this location. This is not seen in [44], which relates to that there partial internal reforming was used as cooling, and the air flow can thus be minimized.
However, our modelling predictions of high stress level at the inlet (Fig. 7) is consistent with failures in cells observed experimentally at the inlet edge, cf. [45,46]. The local cooling from the air inlet would be problematic for both the cell and contact components, as the local contraction would cause tensile stresses in all directions; in the plane of the cell and out of plane for the loss of contact. So these studies are consistent with the issue observed in this work. In the case studies, we investigate the stresses induced by the thermal gradients due to the uniformly distributed flows at the sides of the stack, which is one of the conditions that lead to the formation of high thermal gradients at the edges of the stack (compare Case C and Case D). The high stress concentrations provide a driving force to generate a crack at the weak points [5]. Generally it is accepted that it is the porous less robust structure of the oxygen electrode-contact layer the weaker component, as the contact layer has a relative lower fracture toughness than other ceramic components [3,6,9,43].
To investigate if determined -integrals will lead to a fracture, they can be compared to the critical energy release rate c , as the -integral is an expression for the energy release rate in the case of linear elastic fracture mechanics. The critical energy release rate of commonly used contact materials are in the order of 1 ∼ 2 J m −2 , and for some of contact materials, e.g. CuMn, LSM/Schott composite, this value can be higher, up to 13 J m −2 , cf. [9,10,20] and [29]. The maximum values of the modelled energy release rate ( = ) in the base case, case A and case B of this work are higher than the upper values measured in experiments, while the maximum values of in case D and case E are in the same magnitude. From a mechanical perspective, the maximum -integral values in both co-flow and counter-flow configurations under the current operating conditions are at some regions high enough to initiate cracking in the contact layer for many of the available contact materials.
Of the different attempts to bring down the high -integral values for the co-flow case -using high flow resistance header for the fuel flow appears to be the more effective one. This is because the fuel is poorly distributed over the stack in the base case, and even with the proposed change, the fuel is still poorly distributed over the width of the stack. This is the main reason for the high -integral values found. For the counter-flow configuration smaller -integral values have been found at the boundaries, which is primarily due to the more evenly distribution of the temperature. The reason is that the two 'cold' inlets are located at either end of the stack. However, a central critical region appeared which relates to that the central region is now relatively warmer, i.e. a more non-linear distribution of temperature along the flow direction. Applying a linear thermal profile, assuming a very good fuel and air distribution, leads to -integral values below 1 J m −2 (results not shown). In the current repeating unit, the stiff interconnect is directly contacted with a stiff cell, which is the main reason for the high stresses and -integrals. Using nickel foam as the contact component on the fuel side will decrease the stresses, as the mesh/foam skeleton would be a soft component, which will allow for thermal expansions without generating high stresses, or in other words it will disperse these stresses, thus reducing the concentration of the elastic strain energy and the risk of cracking [42]. Both show that with a proper design mechanical failures in stacks can be approached.

Conclusions
In this work, we present a novel computational efficient approach for modelling local fractures in full solid oxide cell stacks, by using a multiscale modelling approach, involving homogenization and localization. The homogenization approach brings the possibility of modelling complex structures and multiphysics processes in a full solid oxide cell stack simultaneously without sacrificing the computational time. Primary physical quantities (e.g. temperature, current density, species concentration, pressure drop, and stresses) associated with specific processes can be estimated at a macroscopic scale. However, by this approach the internal geometries of the stack is indirectly included, and it is not possible to directly describe local effects in the sub-structures, such as localized stresses.
In this contribution, we extend the homogenization modelling of a solid oxide cell stack to be capable of examining industrially relevant stack geometries, the effect of the flow in the manifolds and local fractures e.g. at the interface between the cell and interconnects. To do so, localization is used in a novel manner, where a link between the average stresses computed in the homogenized stack model and the so-called local -integral in a sub-model is obtained. This is done by imposing the stresses in the homogenized model as boundary conditions for the sub-model. The -integral is an expression of the energy release rate for linear elasticity. The criteria for fracture is that theintegral exceeds the critical energy release rate for the cell-interconnect assembly, which can be measured using standard methods and solid oxide cell stack components.
Furthermore, it is shown how this link between the average stresses and the -integral can be established once and for all, such that the -integral can be evaluated as an analytical expression associated with the average strain energy. It is found that the strain energy associated with the vertical tensile stresses, which is the primary driving force for crack propagation, as expected.
The model is used to study the fracturing in the contact layer between the oxygen electrode and metallic interconnect in a solid oxide cell stack operating in fuel cell mode. With the chosen base case geometry and a stack of 100 cells, rather high -integrals are initially observed. The impacts of changing the external load on the stack, dimension of manifolds, number of cells, flow resistance of headers and direction of flow on the fracture energy are studied. It is found that the non-uniformity of the fuel flow is the primary reason for the high thermal gradients and local high stresses and consequent highintegrals. This is partially remedied by adding a header region with high flow resistance, reducing the -integrals by a factor of two. It is also shown that using counter flow can reduce the -integrals for the specific stack geometry.
The linking of the -integral to the average stresses by an analytical expression is computationally efficient, and the computational efficiency is thus not influenced. In this work, we also include the turbulent flow in the manifolds. Depending on the solution procedure the stack models could thus be solved in 27 min for each new geometry, and for each consequent simulation in a parametric study in 3 min on a workstation computer. Thus, it is feasible to use the presented homogenization model to assist the design of a solid oxide cell stack, evaluating the compatibility of the metallic interconnect, the electrodes, and the operating conditions including the flow configurations.

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.