Thermomechanical analysis of porous solid oxide fuel cell by using peridynamics

Solid oxide fuel cell (SOFC) is widely used in hybrid marine propulsion systems due to its high power output, excellent emission control and wide fuel suitability. However, the operating temperature in SOFC will rise up to 800–1000 °C due to redox reaction among hydrogen and oxygen ions. This provides a suitable environment for ions transporting through ceramic materials. Under such operation temperatures, degradation may occur in the electrodes and electrolyte. As a result, unstable voltage, low capacity and cell failure may eventually occur. This study presents thermomechanical analysis of a porous SOFC cell plate which contains electrodes, electrolytes and pores. A microscale specimen in the shape of a plate is considered in order to maintain uniform temperature loading and increase the accuracy of estimation. A new computational technique, peridynamics, is utilized to calculate the deformations and stresses of the cell plate. Moreover, the crack formation and propagation are also obtained by using peridynamics. According to the numerical results, damage evolution depends on the electrolyte/electrode interface strength during the charging process. For weak interface strength case, damage emerges at the electrode/electrolyte interface. On the other hand, for stronger interface cases, damage emerges on pore boundaries especially with sharp corner.


Introduction
The emission control in marine industry has become a concerning issue in environment protection and it is under the restriction of IMO regulations.According to Marpol Annes VI requirement, emission of nitrogen oxides should reduce nearly 80% after 2016 as compared with that in 2000.On the other hand, emission of sulfur oxides in 2020 is enforced to remain below 10% of that in 2012 [1].Hence, hybrid marine propulsion system becomes popular in marine industry with its higher power output and lower emissions.Marine batteries as one of the most important components of this propulsion system are under extensive investigation in order to increase the battery performance and reduce the price [2].Fuel cell is one of the most promising power storage systems which has various applications due to high energy conversion efficiency, high fuel flexibility, low emission and low material price [3].
As a common fuel cell type, Solid oxide fuel cell (SOFC) is composed of solid electrodes and electrolyte.Nickel, YSZ and perovskite-based LaMnO 3 or LSCF are the most commonly used materials for anode, electrolyte and cathode, respectively [4,5].Typically, the electrodes and electrolyte are made as layers and stacked on each other in a certain sequence.In order to increase the efficiency of electrochemical reaction, electrolyte material is added inside the electrode layers [3].Electrochemical reactions (oxidation of the hydrogen) that are required for generation of the free electrons occur through the Triple Phase Boundaries (TPBs) which are the intersection lines of the metal/non-metal grains and fuel filled pores.The reaction rate and the generated power in the anode and cathode depend on the active TPBs which can be controlled by the amount of the compositions, grain size, and porosity ratio.The reactions taking place in a hydrogen fueled SOFC are given below; In marine environment, SOFC not only requires hydrogen, but also carbon oxides and some hydrocarbon substance which can be obtained from marine waste gases and marine fuel oil.Fossil fuels like gasoline, natural gas, coal gas and liquid petroleum gas could also be regarded as possible fuel of SOFC [6].The redox reaction in SOFC generates large amount of heat which results in an operating temperature ranging from 800 ℃ to 1000 ℃ [7].High temperature not only provides a suitable environment for ion transportation through ceramic electrolyte, but also catalyses side reactions in SOFC.Impurities in marine fuels like hydrogen sulphide, chlorides, alkali compounds and tars may react with part of the electrode and electrolyte materials, forming stable and inert by-product on the three-phase-boundary.This will lead to performance degradation of SOFC and harmful waste gas emission [6,8,9,10].Hence, gas purifying equipment is mandatorily installed before SOFC sections [6].Under large thermal loading, the electrode and electrolyte layer will deform.Since the electrolyte material is added inside the electrode layer, the mechanical properties in each layer are not homogenous.Besides, due to different material properties of electrodes and electrolyte, differences in mechanical deformation are observed at different parts of SOFC [3].Hence, degradation of SOFC may occur which may cause negative influence on fuel cell performance and even failure.
Fracture and fatigue analysis of SOFC is important in battery design and optimisation.CCM is a widely used methodology in solid mechanics.In CCM, PDEs are utilized to define the mechanical state of each material point [11].However, PDEs may encounter difficulties in modelling fracture and failure since PDEs are not defined over discontinuous domains.In order to predict fracture evolution, fracture initiation instant and propagation direction information may be required which can influence the accuracy of the analysis [12].The numerical studies of fracture in SOFC is very limited in the literature.Laurencin et al. [13] utilized a combined finite element and statistical approach of Weibull to predict cell fracture.Although finite element method (FEM) is a well-established technique for determination of stresses, its capability is relatively limited in failure prediction.There are various techniques proposed for failure prediction in FEM framework including remeshing, virtual crack closure technique (VCCT), cohesive zone elements and extended finite element method (XFEM).Although these techniques are suitable for the analysis of certain problems, they have several limitations such as requirement of pre-existing damage, mesh dependency, etc.Therefore, a new numerical analysis method, called peridynamics [14][15][16][17][18][19][20], can be a suitable alternative to evaluate the formation and evolution of material fracture and failure.Peridynamic theory is based on integro-differential equations which are valid along discontinuous domains without any pre-defined condition and criteria.Moreover, pre-existing damage is not required and the crack path is not effected by the discretization.In this study, peridynamics is utilized for the first time in the literature for the thermomechanical analysis of porous SOFCs to predict fracture initiation and propagation.Furthermore, the effect of the electrolyte/electrode interface on the fracture behavior is analyzed during the charging process.Before presenting fracture analysis of SOFC plate, peridynamic theory is briefly described and a validation test is demonstrated in Section 2 and Section 3, respectively.By comparing the results from FEM and PD, the capability of PD in dealing with a thermomechanical problem is illustrated.The SOFC electrode plate model is provided in Section 4 and detailed investigation and discussion of damage evolution in SOFC electrode plate by using PD are given in Section 5.

Peridynamic Theory
Peridynamic theory was originally introduced by Silling in 2000 [11].Peridynamics is a new continuum mechanics formulation.Hence, it is based on continuum assumption where each domain is composed of infinite number of infinitely small material points.According to peridynamics, each material point can interact with other material points inside a certain distance , horizon (see Figure 1).According to the Newton's Second Law, the acceleration of each material point can be obtained by taking into account all the interactions associated with that material point, x where H x represents horizon and b refers to body force density.The pairwise force density, f, for each interaction depends on the initial position and displacement of associated material points.The initial length of the interaction (bond) can be defined as and the relative displacement of the bond after deformation can be expressed as Then, the stretch of a bond can be calculated as For a linear elastic isotropic material, the pairwise force function can be defined as where α is the coefficient of thermal expansion, T avg represents the average temperature of material points x and x', and μ is the failure parameter of the peridynamic bond.c is named as bond constant and for two-dimensional cases, it can be expressed as where E refers to elastic modulus and h is the thickness of the two-dimensional domain.
Figure 1.Material point x interacts with neighbouring points within H x [12].
The concept of LEFM was introduced by Griffith in 1921 which has become one of the solutions when dealing with the fracture problem by using CCM [21].According to LEFM, the propagation of cracks increases the crack surface energy, but decreases the elastic energy of particles.When the strain energy release rate is less than twice of the crack surface energy, then the crack can be arrested [22].Hence, the strain energy release rate has become an additional criterion for the governing equation of CCM.In LEFM, the numerical prediction of crack behaviour depends on the loading, geometry and discretisation, and accurate calculation of the strain energy release rate remains a challenge by using classical continuum mechanics [12].On the other hand, in PD, the nucleation and evolution of the fracture are represented based on bond breakage.Displacement of material points will lead to deformation of bonds.If the material is brittle and the stretch of each bond exceeds a critical value, s c , i.e. critical stretch, the bond breaks and there is no longer interaction between the material points associated with this bond.Hence, the failure parameter can be defined to describe the state of a bond as Then, the local damage of the material point x can be expressed as In other words, local damage is the percentage of the remaining (unbroken) bonds of material point x.φ = 0 means particle x is undamaged and φ = 1 refers to particle x loosing all of its interactions and becoming fully damaged.As more material points being damaged, the crack nucleation and propagation can be observed.
The value of critical stretch depends on critical energy release rate, G c [11,12].G c can be expressed according to peridynamics as On the other hand, G c can also be calculated from fracture toughness K [23] which can be expressed for plane stress condition as By substituting Eqs.(2-10) in Eq. ( 1), the displacement of each material point can be calculated.Moreover, the pressure state and stresses can also be evaluated as part of post-processing operation.Hydrostatic pressure can be derived from Cauchy stress and first Piola-Kirchoff stress.The first Piola-Kirchoff stress can be obtained in peridynamics as Then, Cauchy stress can be calculated as Where   and Hence, the two-dimensional hydrostatic stress of the material point x can be expressed as Moreover, the pressure of material point x can be defined as Since any stress tensor is the summation of hydrostatic stress and deviatoric stress, the deviatoric stress can be obtained as Hence, the two-dimensional von-Mises stress can be calculated as [23] Analytical solution of Eq. ( 1) is generally not possible.Therefore, numerical solution schemes are widely used.In this study, meshless discretization approach is utilised by discretizing the solution domain into finite number of volumes by defining a point at the centre of each volume.The time integration can be performed by using either explicit or implicit time integration schemes.If a static condition is analysed as in this study, the solution can be obtained by using an adaptive dynamic relaxation scheme.

Validation Case
Before we carry out the thermomechanical analysis of porous electrode plate, it is essential to validate the current formulation.Therefore, a simple thermomechanical example is considered as shown in Figure 2.  A square plate is selected as the specimen of our validation case.Before we start the simulation, the plate is free of thermal loading and the initial temperature of the plate is T 0 (0 ℃).The plate is then subjected to thermal loading on its boundaries with temperature T 1 (0.15 ℃).Hence, it is expected that the specimen will expand as a result of this loading condition.An imaginary material, A, is used for simplicity and material properties are given in Table 1.Two different methods are utilized to simulate the thermomechanical process.Hence, peridynamic results are compared against CCM based approach FEM.FEM analysis is performed by utilizing the commercially available software, ANSYS, due to its sophisticated thermal-structural analysis capability.Since the validation test is a thermal expansion problem, 8-node coupled thermal-mechanical element PLANE 223 is selected as element type using plane stress assumption.PD analysis is performed by using an in-house code based on the information in Section 2 and implemented in MATLAB.
In order to reduce the computational time, only a quarter of the geometry is modelled (shown as red contour in Figure 2), since the plate geometry and loading distribution are symmetric.As shown in Figure 3, once the plate is subjected to temperature loading, the boundary regions of the plate expand while the centre of the plate remains undeformed.Results have shown a good agreement between FEM and PD according to the distributions of temperature and displacements.Hence, it can be concluded that PD has the capability for performing accurate thermomechanical analysis.

Thermomechanical Analysis of Porous Electrode Plate
Generally, the SOFC unit is a multilayer structure which is composed of ceramic and metallic materials [24].Electrode layers and electrolyte layers are stacked in a certain sequence.Besides, electrode materials are designed to have excellent catalytic activity and some amount of electrolyte material is inserted inside the electrode layers.In order to obtain a detailed visualization of particle distribution in electrode layers, various advanced technologies such as scanning electron microscopy and X-ray computed tomography have been applied [3,25].In this study, an anode layer (Nickel-YSZ) is considered and the coloured image is shown in Figure 4.In Figure 4, the yellow colour represents the electrolyte (YSZ) material and the green colour represents the electrode (Nickel) material.Black colour in this figure represents pore.Based on this information, the anode layer can be discretised into finite number of points with a certain volume.Detailed discretisation is shown in Figure 5.The dimensions of the anode layer are given as 12 μm × 15.9 μm × 0.12 μm.Since the thickness is relatively small as compared to the length and width scales, the plane stress assumption is utilised.The operating temperature of SOFC varies from 800 ℃ to 1000 ℃ and multiple cases with different thermal loading conditions are considered in order to evaluate the stress evolution and fracture occurrence.

Numerical Results and Discussion
In this case study, we first performed thermomechanical analysis of SOFC layer without considering failure to obtain displacement and von Mises stress distributions to compare the peridynamic results against FEM results.Then, we concentrated on crack evolution.Three different bond strength configurations for the electrode-electrolyte interface (weak interface, uniform interface and strong interface) were considered.Dimensions of electrode plate specimen are specified as 12 μm × 15.9 μm.The top edge of the plate is free whereas the left and right edges of the plate are constrained in horizontal direction.Bottom edge is fully constrained.Since the dimensions of the plate are small, it is assumed that the temperature is uniform throughout the plate.Material properties of electrode and electrolyte are given in Table 2.

Thermomechanical simulation without failure
Normal operating temperature of SOFC varies from 800 ℃ to 1000 ℃.High temperature may cause degradation of the electrode plate.In this case, thermomechanical analysis was performed by disregarding failure.Under uniform thermal loading at 800 ℃, operating condition of the electrode and origin of the damage can be investigated.Displacement and von-Mises stress distribution based on PD are shown in Figure 6a, c.Besides, results are also obtained by using FEM for comparison.According to Table 3 and Figure 6, the results from PD and FEM have reached an agreement for displacements in both horizontal and vertical directions and von-Mises stress.Since some of the boundaries of the plate specimen are constrained, the pores tend to shrink.Hence, high stresses concentrate around the pore regions especially those with sharp tips. Figure 7 show the distribution of von-Mises stress and they indicate the possible positions that cracks may form.The results have an overall agreement between PD and FEM.Slight differences between the two approaches may be due to different visualization tools used for PD and FEM analyses, and utilizing two different ways to calculate stresses in PD and FEM analyses.

Thermomechanical simulation with crack propagation
According to Figure 7, the maximum von-Mises stress on the SOFC plate is around 150 MPa under 800 ℃ thermal loading which is larger than the yield strength of electrode material.Hence, damage may occur at this temperature.In this section, three cases will be presented to demonstrate damage formation and evolution based on the strength of electrode-electrolyte interface, i.e. weak interface, uniform interface and strong interface.Critical stretch values for the bonds crossing electrode-electrolyte interface are assumed to be 0.5s c as weak interface, s c as uniform interface and 2s c as strong interface.Thermal loading increases gradually from 0 ℃ to 1000 ℃ in 10000 steps.
Damage states for the three different interface strength cases in one single charging process are shown in Figures 8-10.The formation and evolution of cracks differ in all three cases.Damage will initiate on high hydrostatic stress region of the plate; such as pores, especially with sharp boundary.As the temperature increases, the amount of damage increases.Maximum damage at 1000 ℃ has reached 16%, 3.1% and 2.5% for weak, uniform and strong interface configurations, respectively.For the weak interface case, due to low critical stretch value, damage prefers to evolve and propagate along the electrode-electrolyte interface.This will lead to the separation between electrode particles and electrolyte particles.However, for uniform and strong interface cases, the separation of electrode and electrolyte is not obvious.Most of the damages form at electrolyte regions and pore boundaries where high hydrostatic stresses are located.Since the fracture toughness of the YSZ is much weaker as compared to Nickel, electrolyte particles damage much easier than electrode particles.Moreover, the damage starts at lower temperatures for the weak interface case.
According to the simulations described above, damage evolution depends on electrode and electrolyte interface configuration and SOFC electrode plate geometry.By applying PD, the damage situation, probable position of fracture and probable direction of fracture evolution of the plate can be predicted in a single charging process.During frequent cycling operations, cracks may form and propagate due to increase in damage.Electron flows from high electric potential region to low electric potential region through electrolyte particles.However, once the electrolyte particles are damaged, electron diffusion path may be blocked by these cracks and damages.Therefore, they may have to find alternative path to transfer to low electric potential region.This may cause an increase of inner resistance and reduction of electrical capacity and thermal stability, which will lead to performance degradation and even failure of the fuel cell.

Conclusions
In order to increase efficiency of electrochemical reaction in fuel cells, electrolyte materials are mixed with electrode material.Since the operating temperature varies from 800 ℃ to 1000 ℃, fuel cell may be damaged due to different material properties of electrode and electrolyte.Under high temperature loading, electrode plate deforms and stresses concentrate at some of the pore boundaries, especially those with sharp tips.Hence, damage may occur at these regions.Damage may also evolve along electrode-electrolyte interfaces especially if the interface strength is weaker.This may degrade the cell performance and even cause failure of the fuel cell.Therefore, fracture analysis of SOFC is essential and peridynamics as a new computational technique is used for the first time in the literature.First, the capability of peridynamics is validated by considering a simple square plate specimen under uniform loading condition.Peridynamic results agree very well with that of finite element method.Then, damage evolution in a porous electrode plate subjected to uniform temperature loading is investigated for different interface strengths of the interface region between electrode and electrolyte.For the weaker interface case, damage starts at lower temperature than the stronger interface cases.Moreover, since the fracture toughness of electrolyte material (YSZ) is weaker than electrode material (Nickel), electrolyte is damaged earlier while electrode remains undamaged.As a future work, this study will be extended to analyze a 3-Dimensional structure to obtain a more realistic model.This work is currently under progress and will be a subject of a future publication.Moreover, only a single charging process is considered.The proposed study can easily be improved by considering multiple charging processes to determine the number of cycles causing total failure of the specimen.

Table 1 .Figure 3 .
Figure 3. Results of thermal expansion: (a) Temperature distribution by FEM (b) Temperature distribution by PD (c) Displacement in x direction by FEM (d) Displacement in x direction by PD (e) Displacement in y direction by FEM (f) Displacement in y direction by PD.

Figure 6 .Figure 7 .
Figure 6.Comparison between PD and CCM: (a) Displacement in horizontal direction by PD, (b) Displacement in horizontal direction by FEM, (c) Displacement in vertical direction by PD, and (d) Displacement in vertical direction by FEM.

Table 2 .
Material properties of anode layer.

Table 3 .
Comparison of maximum deformation and von-Mises stress.