A method for evaluating stress characteristics of an operating state of an electric pile based on finite element analysis

Abstract

The application provides a kind of stress characteristic evaluation method of stack operating state based on finite element analysis, comprising: step S1, according to the physical structure modeling assembly of solid oxide stack forms full-size stack model;Step S2, simulates the co-sintering process of porous electrode and thin layer electrolyte in full-size stack model, obtains residual stress and first stress parameter and calculates first elastic stress and first inelastic stress;Step S3, obtains the assembly stress generated in assembly process and second stress parameter and calculates second elastic stress and second inelastic stress;Step S4, heating solid oxide stack by temperature control furnace, obtains thermal stress and third stress parameter and calculates third elastic stress and third inelastic stress;Step S5, integrate elastic stress and inelastic stress as stack operating state stress characteristic evaluation result.The beneficial effect is that the application can break through the limitation under complex stress coupling scene and quantify the stress superposition effect in different stages.

Description

Technical Field

[0001] This invention relates to the technical field of fuel cell stack reliability assessment, and more specifically, to a method for assessing the stress characteristics of fuel cell stack operating conditions based on finite element analysis. Background Technology

[0002] Solid oxide fuel cell stacks are key energy conversion devices that can operate in both fuel cell mode (power generation) and electrolysis mode (hydrogen production). Their performance and lifespan are significantly affected by multi-field coupling stress.

[0003] During long-term operation in fuel cell and electrolysis modes, solid oxide battery stacks are prone to performance degradation or structural failure due to the coupling effect of thermomechanical stress. Existing technologies suffer from the following main shortcomings: Traditional analytical models typically focus only on single stress factors (such as thermal stress or assembly stress), neglecting the synergistic effect of residual stresses (such as sintering stress during manufacturing). Limited by the high-temperature sealed environment, there is a lack of technical means to monitor the internal stress distribution of the stack in real time. Existing research largely relies on extrapolating small-scale short-stack experimental data to industrial-grade long stacks (such as 30-70 stacks), leading to significant errors in practical applications. More importantly, current research is mostly concentrated on a single operating mode (such as power generation mode), lacking a systematic analysis of the complex thermal stress evolution under dual-mode (electrolysis / power generation) switching, making it difficult to comprehensively guide parameter optimization and reliability design under actual operating conditions. These limitations collectively restrict the long-term stable operation and large-scale application of solid oxide battery stacks.

[0004] Therefore, there is an urgent need for a method to evaluate the stress characteristics of fuel cell stack operation that can overcome the limitations of complex stress coupling scenarios and quantify the stress superposition effect at different stages. Summary of the Invention

[0005] The technical problem to be solved by this invention is how to overcome the limitations of complex stress coupling scenarios and quantify the stress superposition effect at different stages. In order to overcome the defects of the above-mentioned existing technologies (or related technologies), this invention provides a method for evaluating the stress characteristics of electric stack operation based on finite element analysis.

[0006] This invention provides a method for evaluating the stress characteristics of fuel cell stacks under operating conditions based on finite element analysis, comprising the following steps:

[0007] Step S1: Based on the physical structure of the solid oxide fuel cell, perform three-dimensional modeling to obtain multiple component parts, and assemble the component parts to form a full-size fuel cell model;

[0008] Step S2: Simulate the co-sintering process of porous electrodes and thin-layer electrolyte in the full-size electric stack model using molecular dynamics simulation, obtain the residual stress and first stress parameter generated during the co-sintering process, and obtain the first elastic stress and first inelastic stress of the full-size electric stack model based on the finite deformation theory and the second Piola-Kirchhoff stress tensor, according to the residual stress and the first stress parameter.

[0009] Step S3: Obtain the assembly stress and second stress parameters generated by each component during the assembly process, and based on the small deformation theory and the linear elastic constitutive relation, obtain the second elastic stress and second inelastic stress of the full-size fuel cell model according to the assembly stress and the second stress parameters.

[0010] Step S4: The solid oxide fuel cell is heated by a temperature-controlled furnace to obtain the thermal stress and third stress parameters generated during the heating process. Based on the small deformation theory and thermomechanical coupling constitutive relationship, the third elastic stress and the third inelastic stress are obtained according to the thermal stress and the third stress parameters.

[0011] Step S5: The first elastic stress, the first inelastic stress, the second elastic stress, the second inelastic stress, the third elastic stress, and the third inelastic stress are used as the stress characteristic evaluation results of the stack operation state.

[0012] The present invention provides a method for evaluating the stress characteristics of fuel cell stacks under operating conditions based on finite element analysis, which has the following advantages compared with existing technologies:

[0013] In this invention, a full-size fuel cell stack model is constructed through step S1, the first elastic stress and the first inelastic stress are calculated under residual stress through step S2, the second elastic stress and the second inelastic stress are calculated under assembly stress through step S3, the third elastic stress and the third inelastic stress are calculated under thermal stress through step S4, and the stress characteristic evaluation results of the fuel cell stack operating state are obtained through step S5. In the above steps, a full-process coupled model of residual stress and assembly stress, i.e., a full-size fuel cell stack model, is established to provide a basic model for the co-sintering process and the assembly process, which overcomes the limitations of complex stress coupling scenarios. Furthermore, a three-stage coupling framework of "residual stress-assembly stress-thermal stress" is proposed to quantify the stress superposition effect of different stages. This is different from the traditional single-stage analysis and can more accurately evaluate the stress characteristic evaluation results of the fuel cell stack operating state.

[0014] In one possible implementation, step S2 involves obtaining the instantaneous temperature during the co-sintering process, the first coefficient of thermal expansion corresponding to different instantaneous temperatures, the reference temperature under no stress, the unit tensor, the intrinsic strain, the first initial strain, the displacement gradient, the first elastic matrix, and the deformation gradient matrix, and incorporating them into the first stress parameter.

[0015] In one possible implementation, in step S2, the thermal deformation gradient is obtained based on the first thermal expansion coefficient, the instantaneous temperature, the reference temperature, and the unit tensor; the intrinsic deformation gradient is obtained based on the intrinsic strain and the unit tensor; the initial deformation gradient is obtained based on the first initial strain and the unit tensor; then, the first elastic stress is obtained based on the displacement gradient, the thermal deformation gradient, the intrinsic deformation gradient, the initial deformation gradient, the first elastic matrix, and the deformation gradient matrix; and the first inelastic stress is obtained by subtracting the residual stress from the first elastic stress.

[0016] In one possible implementation, in step S2, the thermal deformation gradient, the intrinsic deformation gradient, and the initial deformation gradient are obtained using the following calculation formulas:

[0017] ;

[0018] ;

[0019] ;

[0020] ;

[0021] in,

[0022] This represents the thermal deformation gradient;

[0023] Represents the unit tensor;

[0024] Indicates thermal strain;

[0025] Indicates the first coefficient of thermal expansion;

[0026] Indicates local temperature;

[0027] This indicates the reference temperature;

[0028] This represents the intrinsic deformation gradient;

[0029] Indicates the intrinsic strain;

[0030] This represents the initial deformation gradient;

[0031] This represents the first initial strain.

[0032] In one possible implementation, in step S2, the first elastic stress and the first inelastic stress are obtained by the following calculation formula:

[0033] ;

[0034] ;

[0035] ;

[0036] ;

[0037] in,

[0038] This refers to the residual stress;

[0039] This represents the first elastic stress;

[0040] This represents the first inelastic stress;

[0041] The determinant of the deformed gradient matrix is ​​represented by .

[0042] This indicates an operation that first finds the inverse of a matrix;

[0043] This indicates that the matrix is ​​first inverted and then transposed.

[0044] This represents the first elasticity matrix;

[0045] Indicates elastic strain;

[0046] Indicates the displacement deformation gradient;

[0047] Represents the unit tensor;

[0048] This represents the displacement gradient;

[0049] Indicates the elastic deformation gradient;

[0050] Indicates the gradient of inelastic deformation;

[0051] This represents the thermal deformation gradient;

[0052] This represents the intrinsic deformation gradient;

[0053] This represents the initial deformation gradient;

[0054] It represents elastic strain.

[0055] In one possible implementation, the second stress parameter includes the first elastic modulus and first Poisson's ratio of the materials of each component during the assembly process at room temperature, the first initial strain during the co-sintering process, and the first stack displacement. In step S3, a second elastic matrix is ​​obtained based on the first elastic modulus and the first Poisson's ratio; the total stack strain is obtained based on the first stack displacement; the inelastic strain of the stack is obtained based on the first initial strain; the elastic strain of the stack is obtained based on the total stack strain and the inelastic strain of the stack; the second elastic stress is obtained based on the elastic strain of the stack and the second elastic matrix; and the assembly stress is subtracted from the second elastic stress to obtain the second inelastic stress.

[0056] In one possible implementation, in step S3, the second elastic stress and the second inelastic stress are obtained by the following calculation formula:

[0057] ;

[0058] ;

[0059] ;

[0060] ;

[0061] ;

[0062] ;

[0063] in,

[0064] This indicates the assembly stress;

[0065] This represents the second elastic stress;

[0066] This represents the second inelastic stress;

[0067] This represents the total strain of the fuel cell stack;

[0068] This indicates the displacement of the first fuel cell stack;

[0069] This indicates that the matrix is ​​first inverted and then transposed.

[0070] This indicates the inelastic strain of the fuel cell stack;

[0071] This represents the first initial strain;

[0072] This represents the elastic strain of the fuel cell stack;

[0073] This represents the second elasticity matrix;

[0074] This represents the first elastic modulus;

[0075] This represents the first Poisson's ratio.

[0076] In one possible implementation, the third stress parameter includes the second elastic modulus and second Poisson's ratio of the materials of each component at room temperature during the heating process, the second initial strain during the assembly process, the second coefficient of thermal expansion of each component, the local temperature, the reference temperature under no stress, and the second stack displacement. In step S4, a third elastic matrix is ​​obtained based on the second elastic modulus and the second Poisson's ratio; the total thermal stress strain is obtained based on the second stack displacement; the inelastic thermal stress strain is obtained based on the second initial strain, the second coefficient of thermal expansion, the local temperature, and the reference temperature; the elastic thermal stress strain is obtained based on the total thermal stress strain and the inelastic thermal stress strain; the third elastic stress is obtained based on the elastic thermal stress strain and the third elastic matrix; and the third inelastic stress is obtained by subtracting the thermal stress from the third elastic stress.

[0077] In one possible implementation, in step S4, the third elastic stress and the third inelastic stress are obtained by the following calculation formula:

[0078] ;

[0079] ;

[0080] ;

[0081] ;

[0082] ;

[0083] ;

[0084] ;

[0085] in,

[0086] This indicates the thermal stress;

[0087] This represents the third elastic stress;

[0088] This represents the third inelastic stress;

[0089] This represents the total strain of the thermal stress;

[0090] This indicates the displacement of the second fuel cell stack;

[0091] This indicates that the matrix is ​​first inverted and then transposed.

[0092] This represents the inelastic strain of the thermal stress;

[0093] This represents the second initial strain;

[0094] This represents the elastic strain of the thermal stress;

[0095] Indicates the thermal strain during operation;

[0096] This represents the second coefficient of thermal expansion;

[0097] This indicates the local temperature;

[0098] This indicates the reference temperature;

[0099] This represents the third elasticity matrix;

[0100] This represents the second elastic modulus;

[0101] This represents the second Poisson's ratio.

[0102] In one possible implementation, step S1 further includes:

[0103] For the pole connector in the solid oxide fuel cell, the channels and ribs contained in the pole connector are homogenized and modeled as a porous channel-rib region to obtain a multiphysics fuel cell CFD model, and the multiphysics fuel cell CFD model replaces the full-size fuel cell model. Attached Figure Description

[0104] Figure 1 This is a flowchart of the steps of the present invention;

[0105] Figure 2 This is a schematic diagram of the assembly process of the solid oxide fuel cell of the present invention;

[0106] Figure 3 This is a schematic diagram of the geometric structure of the solid oxide fuel cell of the present invention;

[0107] Figure 4 This is a schematic diagram of the "flow channel-rib" homogenization method in the CFD model of the present invention;

[0108] Figure 5 This is a schematic diagram of the main chemical reactions and transport processes of the solid oxide fuel cell of the present invention;

[0109] Figure 6 This is a schematic diagram of the mesh generation for the full-size electric stack stress model of the present invention;

[0110] Figure 7 This is a schematic diagram of the mesh generation for the multiphysics stack CFD model of the present invention;

[0111] Figure 8 This is a comparison of the current-voltage curves between the detailed CFD model and the homogenized CFD model of the one-layer solid oxide fuel cell of the present invention.

[0112] Figure 9 This is a comparison of the H2 mole fraction and temperature distribution between the detailed CFD model and the homogenized CFD model of the one-layer solid oxide fuel cell of the present invention. Detailed Implementation

[0113] First, those skilled in the art should understand that these embodiments are merely illustrative of the technical principles of the present invention and are not intended to limit the scope of protection of the present invention. Those skilled in the art can make adjustments as needed to adapt to specific application scenarios.

[0114] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0115] See Figure 1 and Figure 3This invention discloses a method for evaluating the stress characteristics of a fuel cell stack under operating conditions based on finite element analysis, comprising the following steps:

[0116] Step S1: Based on the physical structure of the solid oxide fuel cell, a three-dimensional model is created to obtain multiple components, and the components are assembled to form a full-size fuel cell model.

[0117] Step S2: The co-sintering process of porous electrodes and thin-layer electrolyte in the full-size electric stack model is simulated by molecular dynamics to obtain the residual stress and first stress parameters generated during the co-sintering process. Based on the finite deformation theory and the second Piola-Kirchhoff stress tensor, the first elastic stress and the first inelastic stress of the full-size electric stack model are obtained according to the residual stress and the first stress parameters.

[0118] Step S3: Obtain the assembly stress and second stress parameters generated by each component during the assembly process, and based on the small deformation theory and linear elastic constitutive relation, obtain the second elastic stress and second inelastic stress of the full-size fuel cell model according to the assembly stress and second stress parameters.

[0119] Step S4: The solid oxide fuel cell is heated by a temperature-controlled furnace to obtain the thermal stress and third stress parameters generated during the heating process. Based on the small deformation theory and thermomechanical coupling constitutive relationship, the third elastic stress and third inelastic stress are obtained according to the thermal stress and third stress parameters.

[0120] Step S5: The first elastic stress, the first inelastic stress, the second elastic stress, the second inelastic stress, the third elastic stress, and the third inelastic stress are used as the stress characteristic evaluation results of the stack operation state.

[0121] In this embodiment of the invention, step S1 further includes:

[0122] For the pole connectors in solid oxide fuel cells, the channels and ribs contained in the pole connectors are homogenized and modeled as porous channel-rib regions to obtain a multiphysics fuel cell CFD model. This multiphysics fuel cell CFD model then replaces the full-size fuel cell model. The multiphysics fuel cell CFD model is shown below. Figure 6 As shown.

[0123] In this embodiment of the invention, considering that the assembly process of the solid oxide fuel cell is one of the important steps, including assembly alignment and bolt tightening, such as... Figure 2As shown, it is worth noting that the battery will bulge and deform in the middle after sintering during the manufacturing process. Therefore, it needs to be flattened with a 2 kg pressure block before the solid oxide battery is assembled. The flattened battery (with the air electrode facing up) is placed on the lower electrode connector of the solid oxide battery and aligned with the reference position. Then, the lower sealant, metal frame, upper sealant and upper electrode connector are placed in sequence. Finally, the assembled solid oxide battery is bolted. The bolt tightening force used in this embodiment is 5844 Newtons.

[0124] In this embodiment of the invention, the core components of the solid oxide fuel cell include a porous fuel electrode Ni-8YSZ and an air electrode La. 0.6 Sr 0.4 Co 0.2 Fe 0.8 O 3–δ –Gd 0.1 Ce 0.9 O 1.95 The solid oxide fuel cell (SOCC) stack, with its dense electrolyte layer 8YSZ and three-phase boundary (TPB) region, is where electronically conductive, ionicly conductive, and gaseous phases coexist, providing the main active sites for electrochemical reactions. The relatively high operating temperature (typically 600 to 1000 degrees Celsius) allows for the direct use of hydrocarbon fuels through in-situ steam reforming without external pretreatment. However, this temperature range also allows for the simultaneous occurrence of multiple chemical reactions within the SOCC stack. When using hydrocarbon fuels, complex reforming reactions (such as so-called co-electrochemical reactions) can occur. Figure 4 As shown, this embodiment mainly considers the reforming reaction of methane in a solid oxide fuel cell, namely the methane steam reforming (MSR) and water-gas shift reforming (WGSR) reactions.

[0125] In this embodiment of the invention, the thermal stress simulation of the solid oxide fuel cell under operating conditions is based on the temperature field predicted by CFD modeling coupled with the electrochemical reaction multiphysics process. However, there is a significant strong nonlinear coupling relationship between the electrochemical reaction and the multiphysics process. In addition, each pole connector contains a complex structure with 66 channels and 68 ribs. When more than 20 repeating elements are assembled in a long fuel cell structure, this will bring significant computational complexity to geometric modeling, mesh generation and interface processing. To solve this problem, while preserving the structural details of the solid oxide fuel cell as much as possible, this embodiment develops and applies a homogenization modeling method, which regards the channels and ribs as porous channel-rib regions for modeling the multiphysics fuel cell CFD model. That is, in the multiphysics fuel cell CFD model, an equivalent homogeneous medium is assumed to represent the non-uniform ribs and channels in the pole connector.

[0126] In this embodiment of the invention, the manufacturing process of the battery cells in the solid oxide stack mainly involves high-temperature sintering and subsequent flattening at room temperature. Specifically, air electrode composite material is deposited on the electrolyte surface of the sintered half cell using screen printing technology (i.e., the first sintering step), and then sintering is performed to form a complete full cell structure (i.e., the second sintering step).

[0127] During the manufacturing process of solar cells in solid oxide fuel cells, although the absolute deformation occurs at the millimeter scale, the relative deformation becomes quite large due to the total thickness of the solar cell being 630 micrometers. For example, the deformation displacement caused by sintering can reach 8.04 millimeters, which is about 12.76 times the total thickness of the solar cell. Under such conditions, the Frobenius norm (‖∇u‖) of the displacement gradient tensor may exceed the threshold applicable to the small deformation theory (i.e., ‖∇u‖≪1).

[0128] In this case, geometric nonlinearity needs to be considered when |∇u|>0.1, which makes the Cauchy stress tensor-based approach more complex. σ The linear constitutive model cannot fully describe the nonlinear deformation behavior. Therefore, a finite deformation framework is adopted in this embodiment: the initial undeformed state is defined as the reference configuration, the kinematics is controlled by the total Lagrange formula, and the constitutive equation integrates the second Piola-Kirchhoff stress tensor.

[0129] In this embodiment of the invention, in step S2, the instantaneous temperature during the co-sintering process, the first thermal expansion coefficient corresponding to different instantaneous temperatures, the reference temperature under no stress, the unit tensor, the intrinsic strain, the initial strain, the displacement gradient, the first elastic matrix, and the deformation gradient matrix are obtained and included in the first stress parameter. In step S2, the thermal deformation gradient is obtained based on the first thermal expansion coefficient, the instantaneous temperature, the reference temperature, and the unit tensor; the intrinsic deformation gradient is obtained based on the intrinsic strain and the unit tensor; and the initial deformation gradient is obtained based on the first initial strain and the unit tensor. Subsequently, the first elastic stress is obtained based on the displacement gradient, the thermal deformation gradient, the intrinsic deformation gradient, the initial deformation gradient, the first elastic matrix, and the deformation gradient matrix. The residual stress is subtracted from the first elastic stress to obtain the first inelastic stress.

[0130] In this embodiment of the invention, residual stress can be directly obtained using COMSOL finite element software based on loading conditions (such as temperature change), material constitutive relations (such as elastic modulus), and deformation gradient. The methods for obtaining assembly stress and thermal stress are the same as those for residual stress, and will not be elaborated further. The expression for residual stress is as follows:

[0131] The first elastic stress and the first inelastic stress are obtained by the following calculation formulas:

[0132] ;

[0133] in, This refers to the residual stress; This represents the first elastic stress; This represents the first inelastic stress;

[0134] The first elastic stress is calculated using the following formula:

[0135] ;

[0136] ;

[0137] ;

[0138] in, The determinant of the deformed gradient matrix is ​​represented by . This indicates an operation that first finds the inverse of a matrix; This indicates that the matrix is ​​first inverted and then transposed. This represents the first elasticity matrix; Indicates elastic strain; Indicates the displacement deformation gradient; Represents the unit tensor; This represents the displacement gradient; Indicates the elastic deformation gradient; Indicates the gradient of inelastic deformation; This represents the thermal deformation gradient; This represents the intrinsic deformation gradient; This represents the initial deformation gradient; Indicates elastic strain;

[0139] The intrinsic deformation gradient and the initial deformation gradient are calculated using the following formulas:

[0140] ;

[0141] ;

[0142] in, This represents the intrinsic deformation gradient; Indicates the intrinsic strain; This represents the initial deformation gradient; This represents the first initial strain. When evaluating the entire fuel cell stack sintering process, the initial strain used is the residual strain inside the half-cell after the first sintering. However, when simulating the entire fuel cell stack assembly process, the initial strain is the residual strain inside the entire cell after the second sintering stage. Furthermore, it should be noted that... The strain represents the small deformation strain describing volume change, with its non-zero component appearing only on the diagonal. During the sintering process of the fuel cell stack, particle aggregation and neck formation lead to planar shrinkage of the sintered fuel cell stack. The shrinkage rate along the x-axis (i.e., the horizontal length direction) was emphasized and quantitatively evaluated using molecular dynamics (MD) methods to assess the shrinkage effect of particle clusters, which is further correlated with intrinsic strain.

[0143] The thermal deformation gradient is calculated using the following formula:

[0144] ;

[0145] ;

[0146] in, This represents the thermal deformation gradient; Represents the unit tensor; Indicates thermal strain; Indicates the first coefficient of thermal expansion; Indicates local temperature; The reference temperature is indicated; the sintering process of the battery cell includes multiple stages such as heating, holding and cooling. The inherent temperature gradient formed in each stage will drive the sintered structure to undergo coordinated deformation under the corresponding thermal environment.

[0147] In this embodiment of the invention, during assembly, the flattened battery cells are tightly sealed to other components of the solid oxide fuel cell stack (such as the frame and sealant) via a bolt pre-tightening device. Due to the rigid constraint applied by the bolt pre-tightening force, the deformation of each component is limited to a small deformation range. Assuming that the mechanical behavior of all components within the solid oxide fuel cell stack conforms to the small deformation theory assumption, in this case, the Cauchy stress tensor... σ Effectively utilizes the true stress distribution in current solid oxide fuel cell configurations for description, i.e., a linearized Cauchy strain tensor. ε The linear elastic constitutive relation can reasonably describe the local stress response caused by the pre-tightening of solid oxide fuel cell assembly.

[0148] In this embodiment of the invention, the second stress parameter includes the elastic modulus, Poisson's ratio, initial strain, and stack displacement of the materials of each component during the assembly process at room temperature. In step S3, the second elastic matrix is ​​obtained based on the elastic modulus and Poisson's ratio, the total stack strain is obtained based on the stack displacement, the inelastic strain of the stack is obtained based on the initial strain, the elastic strain of the stack is obtained based on the total stack strain and the inelastic strain of the stack, and the second elastic stress and the second inelastic stress are obtained based on the elastic strain of the stack and the second elastic matrix.

[0149] In this embodiment of the invention, the second elastic stress and the second inelastic stress are calculated using the following formulas:

[0150] ;

[0151] ;

[0152] ;

[0153] ;

[0154] ;

[0155] in, This indicates the assembly stress; This represents the second elastic stress; This represents the second inelastic stress; This represents the total strain of the fuel cell stack; This indicates the displacement of the first fuel cell stack; This indicates that the matrix is ​​first inverted and then transposed. This indicates the inelastic strain of the fuel cell stack; This represents the first initial strain; This represents the elastic strain of the fuel cell stack; This represents the second elastic matrix; under room temperature conditions without a supply of reactive gas, the stack maintains thermal equilibrium due to the lack of electrochemical reaction, therefore the thermal expansion effect can be ignored;

[0156] The second elasticity matrix is ​​calculated using the following formula:

[0157] ;

[0158] in, This represents the second elasticity matrix; This represents the first elastic modulus; This represents the first Poisson's ratio.

[0159] In this embodiment of the invention, during the operation of the fuel cell stack, due to the combined effects of thermophysical / mechanical property mismatch, external structural constraints, and material anisotropy, significant thermomechanical stress will be generated inside the components of the solid oxide fuel cell stack (such as interconnects, sealants, and frames). Furthermore, under the rigid constraint of bolt preload, it is assumed that the components of the solid oxide fuel cell stack follow the linear elastic small deformation theory. Within this theoretical framework, the induced thermal stress distribution of the solid oxide fuel cell stack is described by the Cauchy stress tensor, while the thermomechanical coupling constitutive equation is established using the linearized Cauchy strain tensor.

[0160] In this embodiment of the invention, the third stress parameter includes the second elastic modulus and second Poisson's ratio of the materials of each component at room temperature during the heating process, the second initial strain during the assembly process, the second thermal expansion coefficient of each component, the local temperature, the reference temperature under no stress, and the second electric stack displacement. In step S4, the third elastic matrix is ​​obtained based on the second elastic modulus and the second Poisson's ratio, the total thermal stress strain is obtained based on the second electric stack displacement, the inelastic thermal stress strain is obtained based on the second initial strain, the second thermal expansion coefficient, the local temperature, and the reference temperature, the elastic thermal stress strain is obtained based on the total thermal stress strain and the inelastic thermal stress strain, the third elastic stress is obtained based on the elastic thermal stress strain and the third elastic matrix, and the third inelastic stress is obtained by subtracting the thermal stress from the third elastic stress.

[0161] In this embodiment of the invention, the third elastic stress and the third inelastic stress are calculated using the following formulas:

[0162] ;

[0163] ;

[0164] ;

[0165] ;

[0166] ;

[0167] ;

[0168] ;

[0169] in, This indicates the thermal stress; This represents the third elastic stress; This represents the third inelastic stress; This represents the total strain of the thermal stress; This indicates the displacement of the second fuel cell stack; This indicates that the matrix is ​​first inverted and then transposed. This represents the inelastic strain of the thermal stress; This represents the second initial strain; This represents the elastic strain of the thermal stress; Indicates the thermal strain during operation; This represents the second coefficient of thermal expansion; This indicates the local temperature; This indicates the reference temperature; This represents the third elasticity matrix; This represents the second elastic modulus; The second Poisson's ratio is indicated. The material properties of each component in the solid oxide fuel cell are updated to the evaluation values ​​at the fuel cell operating temperature (e.g., 1023 K). The material properties are shown in Table 1 below:

[0170] Table 1. Material properties of each component in a solid oxide fuel cell.

[0171] ;

[0172] The third elasticity matrix can be obtained by substituting the values ​​corresponding to the second elastic modulus and the second Poisson's ratio in Table 1. The same values ​​can be used to calculate the second elasticity matrix.

[0173] In this embodiment of the invention, the local temperature is predicted by a multiphysics stack CFD model describing the electrochemical reaction coupled with multiple physics processes. Its thermal inhomogeneity stems from the following synergistic mechanisms: First, the electrochemical oxidation reaction at the three-phase boundary (TPB) releases Gibbs free energy, serving as the core heat source driving the temperature gradient formation; second, variations in gas distribution characteristics in the channels and manifolds, as well as changes in stack operating performance, lead to spatial asymmetry in the distribution of reactant gases (fuel and oxidant), indirectly affecting the temperature distribution; third, in the case of using hydrocarbon fuels (e.g., methane) or co-electrochemical reactions, the endothermic effect of the fuel reforming process and the exothermic effect of the electrochemical oxidation reaction form competing heat sources in the fuel electrode, and their spatial mismatch may trigger local temperature oscillations. Furthermore, forced convection-diffusion heat transfer in the porous electrode and the Joule heating effect at the electrode / electrolyte interface (caused by ion / electron conduction losses) together constitute complex factors influencing heat transfer. These thermal behaviors are quantitatively characterized by a heat energy conservation equation that couples heat transfer and integrates the contributions of multiple heat generation / consumption sources, as shown below:

[0174] ;

[0175] The left-hand side of the equation is due to macroscopic fluid motion (velocity field) The convective term of heat transport caused by enthalpy describes enthalpy ( ) The advection of ); the first term on the right is due to the temperature gradient ( The diffusion term of heat conduction driven by thermal conductivity follows Fourier's law, which is modified by local thermal conductivity; the second term on the right is the heat source term, quantified by the volumetric heat energy generation or dissipation rate. It is the density of the mixed gas. It is the overall speed. It is the specific heat capacity of the fluid. It is the effective thermal conductivity. It is a general heat source item.

[0176] In this embodiment of the invention, the heat source term in the porous electrode can be calculated as follows:

[0177] ;

[0178] in, and These are the effective electronic / ionic conductivity, respectively. , and These are the electronic, ionic, and reversible electrochemical potentials, respectively. The reversible electrochemical potential of the electrode can be expressed as the galvanic potential modulated by the operating temperature, pressure, and gas component concentration. It represents the electrochemically active specific surface area; Current density; To activate the over-points.

[0179] In this embodiment of the invention, if a hydrocarbon fuel (such as methane) is used, the heat source term of the fuel electrode will include a contribution from the reforming reaction, expressed as follows: ,in and The reaction rates are methane steam reforming (MSR) and water-gas shift reforming (WGSR), respectively. and These represent the corresponding enthalpy changes, reaction rates of MSR and WGSR, respectively. and It can be obtained based on the following calculation formula:

[0180] ;

[0181] ;

[0182] in, and These are the chemical reaction rate constants for the forward reactions of methane steam reforming (MSR) and water-gas shift reforming (WGSR), respectively. and These are the equilibrium constants for MSR and WGSR, respectively.

[0183] In this embodiment of the invention, the heat source term in the electrolyte layer can be calculated using the following formula:

[0184] ;

[0185] The heat source term in interconnects and current collectors can be calculated using the following formula:

[0186] .

[0187] In this embodiment of the invention, considering that elementary reactions in electrochemistry (including charge transfer processes) can be modeled through detailed electrochemical kinetics or rate-limiting steps, the modified Butler-Wolmer equation derived from elementary heterogeneous electrochemical reactions (assuming each air electrode and fuel electrode reaction has a single rate-limiting step) is applied to effectively represent the overall reaction. The current density described by the modified Butler-Wolmer equation... With activation overpotential The relationship between these two conditions applies to both hydrogen fuel cells and hydrogen production via water electrolysis, and the specific expression is as follows:

[0188] ;

[0189] ;

[0190] in, and These are the current densities of the fuel electrode and the air electrode. It corresponds to the exchange current density for each electrochemical reaction. It is the charge transfer coefficient.

[0191] In this embodiment of the invention, for fuel cell mode (e.g., hydrogen fuel), the exchange current density is... and The fuel electrode can be represented as:

[0192] ;

[0193] ;

[0194] in, and These are empirical constants used to fit numerical results to experimental data. and It can be calculated using the following formula:

[0195] ;

[0196] ;

[0197] in, and It refers to the pre-factor ( , ), and It is the activation energy ( , ), It is the surface activity density ( ), It is the adhesion probability ( ), It is the molar mass of the gas;

[0198] The activation overpotential can be expressed as:

[0199] .

[0200] In this embodiment of the invention, a homogenization method is used for the channel-rib structure to simplify the multi-physics transmission process. The homogenized region of the channel-rib is conceptualized as an equivalent porous medium with directional permeability, wherein the porosity ( The ratio of the gas phase volume to the total control volume is defined as the homogenization model in the main flow direction (i.e., the x-axis, characteristic length). ) and rib height direction (i.e., z-axis, characteristic height) Maintaining geometric consistency with the original heterogeneous structure in the horizontal direction (i.e., the y-axis, total width) On the surface, the gas channels and ribbed solids are periodically distributed, with the gas channels occupying a proportion of the width. The ribs occupy the remaining Porosity of homogenized regions It can be calculated using the following formula:

[0201] ;

[0202] in, Let be the volume of the gas flow channel, and This is the total volume of the gas flow channel and ribs (i.e., the total volume of the area under consideration).

[0203] By introducing the orthogonal anisotropic permeability tensor Adjusting this tensor allows for the capture of the directional characteristics of gas flow. Adjusting the effective permeability tensor can simulate the flow of gas under different pressure conditions through physically repeating channel geometry. For rectangular channels, along... x Axial permeability component It can be represented as follows:

[0204] ;

[0205] in, The hydraulic diameter of the gas flow channel. This is the shape correction factor for the flow channel cross-section;

[0206] ;

[0207] ;

[0208] The pressure drop of gas with average velocity flowing in laminar flow in a straight channel can be described by applying the Darcy-Weisbach equation:

[0209] ;

[0210] in, The Darcy friction factor is calculated using the following formula:

[0211] ;

[0212] in, It is the Reynolds number.

[0213] In this embodiment of the invention, within the assumed porous channel-rib region, the gas flow path is simplified to a straight channel, ignoring geometric curvature (the tortuosity factor is explicitly set to 1). This simplification indicates that the flow resistance is controlled by the fluid viscosity effect, and structural interactions within the porous region do not contribute additionally, resulting in effective thermal conductivity. Following the mixing rule: multiply the vapor phase thermal conductivity by the porosity, and add the rib material thermal conductivity multiplied by the rib volume fraction, as shown below:

[0214] ;

[0215] in, The thermal conductivity of the ribs (i.e., interconnects) The thermal conductivity of the corresponding gas flowing through the channel;

[0216] Effective conductivity Calculated as the product of the fin volume fraction and its intrinsic conductivity:

[0217] ;

[0218] in, It is the electrical conductivity of the rib.

[0219] In this embodiment of the invention, the entire modeling process is divided into three stages: The first stage is full-size cell co-sintering and subsequent flattening simulation. In this stage, the sintered half-cell and LSCF-GDC air electrode are further co-sintered to form a complete point pair, and the entire stack flattening process is implemented. After sintering, the stress and strain of the functional layer are set as the initial stress and strain of the corresponding electrode during the flattening process, fully considering the thermodynamic changes caused by sintering. The second stage is the stack assembly process. At this time, the stress and strain of each functional layer of the flattened solid oxide stack are used as the initial conditions for solid oxide stack assembly. The third stage is the stack operation process. The stress and strain of the flattened solid oxide stack and the assembly of each component of the solid oxide stack are set as the initial stress and strain for solid oxide stack operation. These initial conditions are thermodynamically coupled with the temperature distribution obtained from the multiphysics stack CFD model.

[0220] In this embodiment of the invention, the boundary condition for evaluating the residual stress of the solid oxide fuel cell is as follows: during full cell sintering, the cell surface is set as a free surface to suppress the rigid body motion of the entire cell, and the air electrode... The temperatures were set at 293.15 K (heating stage), 1348.15 K (holding stage), and 1348.15 K (cooling stage), respectively. The thermophysical properties of the electrode materials were derived from the molecular dynamics (MD) simulation predictions for each stage. During flattening, a boundary force was applied to the top of the air electrode, while the remaining surfaces were treated as free surfaces to suppress rigid body motion of the entire battery. The battery's thermophysical properties were consistent with the MD predictions at the end of sintering. The boundary conditions for evaluating the assembly stress of the solid oxide fuel cell were as follows: mechanical constraints were set by fixing the bottom end plate. The solid oxide fuel cell was tightly fixed by the upper and lower end plates and eight bolts with a preload of 5844 N. The thermophysical properties of each component were the corresponding values ​​at room temperature, and the battery's thermophysical properties remained consistent with the MD predictions at the end of sintering. The boundary conditions for evaluating the thermal stress of the solid oxide fuel cell were similar to those during the assembly process, but the thermophysical properties of each component were updated values ​​at the predicted local temperature. Based on the operating temperature of the solid oxide fuel cell stack, the model assumptions include: the battery is completely flat after being flattened; all components in the solid oxide fuel cell stack are in ideal contact; there are no cracks in the solid oxide fuel cell stack; and the materials of all components in the solid oxide fuel cell stack are linearly elastic and isotropic.

[0221] In this embodiment of the invention, after assembly, the solid oxide fuel cell stack is placed in a temperature-controlled furnace, with the upper and lower surfaces covered with ceramic fiber insulation material, and the peripheral surfaces exposed. At this time, the upper surface of the top plate and the lower surface of the bottom plate are set to thermal insulation conditions, and the sealant and the peripheral surfaces of the frame are allowed to undergo convective and radiative heat exchange with the environment. The model assumes that: the solid oxide fuel cell stack is in steady-state operation; the reactive sites in the porous electrode are uniformly distributed, and the electronic and ion conduction phases are continuous and homogeneous; the NiO-YSZ of the fuel electrode is completely reduced to Ni-YSZ after being heated to the operating temperature, and the volume change caused by the Ni-NiO redox cycle is ignored.

[0222] In embodiments of the present invention, such as Figure 5 As shown, based on the geometric characteristics of the components of the solid oxide fuel cell stack, a differentiated discretization scheme is adopted. In the full-size fuel cell stack model, a parametric scanning algorithm is used to generate a structured hexahedral mesh in orthogonal geometric regions (such as air electrodes and interconnects); an unstructured triangular mesh method is used in the four rounded transition regions of the fuel electrode and electrolyte layer; and an adaptive tetrahedral mesh refinement technique is used for complex geometries such as the metal frame and fasteners.

[0223] In this embodiment of the invention, in the CFD model of the fuel cell stack used to evaluate temperature and other multi-physics parameters, a boundary layer refinement method is implemented for the channels / ribs and manifolds to analyze various transport phenomena. However, numerical experiments show that if a detailed modeling method is used to solve all multi-physics transport processes for solid oxide fuel cell stacks with more than 10 sections, the computer memory required by the nonlinear solver will limit the computational scalability. By applying the channel-rib homogenized fuel cell stack CFD modeling method, the computing power of the same workstation can be expanded to 70 fuel cell stacks. The specific number of meshes and the corresponding computation time for each model are shown in Table 2 below.

[0224] Table 2. Comparison of total grid count and corresponding computation time for detailed and simplified CFD models of solid oxide fuel cells.

[0225] ;

[0226] As shown in Table 2 above, the advantages of using a simplified CFD model of the fuel cell stack to solve for the temperature distribution are very obvious.

[0227] In this embodiment of the invention, a comparative analysis was conducted on the homogenized CFD model and the multiphysics stack CFD model developed in this embodiment for a single-layer solid oxide fuel cell stack. The focus was on evaluating the consistency of the two models in terms of polarization characteristics and key physical parameter distribution patterns to verify the feasibility of the channel-rib homogenization modeling method. The model verification conditions were set as fuel cell operating conditions: fuel gas corresponding to Case A (0.5 SLM) was supplied to the fuel manifold, and air of 1.5 SLM was introduced into the air manifold, maintaining the operating temperature at 1073 K.

[0228] like Figure 7 and Figure 8 As shown, numerical results indicate that within the battery operating voltage range of 0.6–1.05 V, the polarization curves obtained from the detailed model and the simplified model agree well, with a maximum relative error of 5.4% for the current (22.19 A vs. 20.99 A) and a root mean square error (RMSE) of 4.3% within the test range. Both models show a decrease in hydrogen concentration near the air manifold inlet region. The minimum mole fraction is 46.9% for the homogenized model and 47.8% for the detailed model, with a relative deviation of 1.20%. Both models also exhibit the same thermal distribution pattern. The homogenized CFD model predicts the highest temperature in the central stack region to be 1093.1 K, while the detailed CFD model predicts 1093.7 K, with a relative error of 0.15%. The comparison results show that the simplified stack model based on the channel-rib homogenization method maintains good prediction accuracy for multi-physics processes while reducing the computation time from 11 hours to 0.7 hours, effectively overcoming the scalability limitations of traditional finite element CFD methods in long stack simulations.

[0229] In the description of this invention, the references to "one embodiment," "some embodiments," "in this embodiment," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples. Moreover, without contradiction, those skilled in the art can combine and integrate the different embodiments or examples described in this specification, as well as the features of different embodiments or examples.

[0230] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A stress characteristic evaluation method of an operating state of an electric pile based on finite element analysis, characterized by, The method comprises the following steps: Step S1, obtaining a plurality of component parts by three-dimensional modeling according to the physical structure of the solid oxide stack, and assembling the component parts to form a full-size stack model; Step S2, obtaining residual stress and first stress parameters generated in the co-sintering process of the porous electrode and the thin-layer electrolyte in the full-size stack model by molecular dynamics simulation, and obtaining first elastic stress and first inelastic stress of the full-size stack model according to the residual stress and the first stress parameters based on the finite deformation theory and the second Piola-Kirchhoff stress tensor; Step S3, obtaining assembly stress and second stress parameters generated in the assembly process of the component parts, and obtaining second elastic stress and second inelastic stress of the full-size stack model according to the assembly stress and the second stress parameters based on the small deformation theory and the linear elastic constitutive relation; Step S4, heating the solid oxide stack by a temperature control furnace, obtaining thermal stress and third stress parameters generated in the heating process, and obtaining third elastic stress and third inelastic stress according to the thermal stress and the third stress parameters based on the small deformation theory and the thermal-mechanical coupling constitutive relation; Step S5, taking the first elastic stress, the first inelastic stress, the second elastic stress, the second inelastic stress, the third elastic stress and the third inelastic stress as the stack operating state stress characteristic evaluation result.

2. The stress signature evaluation method of claim 1, wherein, In the step S2, the instantaneous temperature in the co-sintering process, the first thermal expansion coefficient corresponding to different instantaneous temperatures, the reference temperature under no stress, the unit tensor, the intrinsic strain, the first initial strain, the displacement gradient, the first elastic matrix and the deformation gradient matrix are obtained and included in the first stress parameters.

3. The stress signature evaluation method of claim 2, wherein, In the step S2, the thermal deformation gradient is obtained according to the first thermal expansion coefficient, the instantaneous temperature, the reference temperature and the unit tensor, the intrinsic deformation gradient is obtained according to the intrinsic strain and the unit tensor, the initial deformation gradient is obtained according to the first initial strain and the unit tensor, then the first elastic stress is obtained according to the displacement gradient, the thermal deformation gradient, the intrinsic deformation gradient, the initial deformation gradient, the first elastic matrix and the deformation gradient matrix, and the first inelastic stress is obtained by subtracting the residual stress from the first elastic stress.

4. The stress signature evaluation method of claim 3, wherein, In the step S2, the thermal deformation gradient, the intrinsic deformation gradient and the initial deformation gradient are obtained by the following calculation formula: ； ； ； ； Wherein, represents the thermal deformation gradient; denotes the unit tensor; represents thermal strain; denotes the first coefficient of thermal expansion; represents the local temperature; represents the reference temperature; denotes the intrinsic deformation gradient; denotes the intrinsic strain; represents the initial deformation gradient; denotes the first initial strain.

5. The stress signature evaluation method of claim 3, wherein, In the step S2, the first elastic stress and the first inelastic stress are obtained by the following calculation formula: ； ； ； ； Wherein, represents the residual stress; represents the first elastic stress; represents the first non-elastic stress; denotes the determinant of the deformation gradient matrix; denotes a pre-inverse matrix operation on the matrix; denotes a transpose matrix operation on the matrix after pre-inversion of the matrix; denotes the first elastic matrix; represents the elastic strain; Indicates the displacement deformation gradient; denotes the unit tensor; representing the displacement gradient; represents the gradient of elastic deformation; denotes the gradient of inelastic deformation; represents the thermal deformation gradient; denotes the intrinsic deformation gradient; represents the initial deformation gradient; represents the elastic strain.

6. The finite element analysis-based stress profile evaluation method of a state of operation of an electrical stack according to claim 1, characterized in that, The second stress parameter includes a first elastic modulus, a first Poisson's ratio, a first initial strain and a first stack displacement of each component material at room temperature during the assembly process, in the step S3, a second elastic matrix is obtained according to the first elastic modulus and the first Poisson's ratio, a total stack strain is obtained according to the first stack displacement, an inelastic stack strain is obtained according to the first initial strain, an elastic stack strain is obtained according to the total stack strain and the inelastic stack strain, the second elastic stress is obtained according to the elastic stack strain and the second elastic matrix, and the second inelastic stress is obtained by subtracting the assembly stress from the second elastic stress.

7. The finite element analysis-based stress profile evaluation method of a state of operation of an electrical stack according to claim 6, characterized in that, The second elastic stress and the second inelastic stress are obtained by the following calculation formula in the step S3: ； ； ； ； ； ； Wherein, represents the assembly stress; represents the second elastic stress; represents the second inelastic stress; represents the total strain of the stack; represents the first stack displacement; denotes a transpose matrix operation on the matrix after pre-inversion; represents the inelastic strain of the stack; represents the first initial strain; represents the elastic strain of the stack; denotes the second elastic matrix; represents the first elastic modulus; represents the first Poisson's ratio.

8. The finite element analysis-based stress profile evaluation method of a state of operation of an electrical stack according to claim 1, characterized in that, The third stress parameter includes a second elastic modulus, a second Poisson's ratio, a second initial strain during the assembly process, a second thermal expansion coefficient of each component, a local temperature, a reference temperature without stress and a second stack displacement, in the step S4, a third elastic matrix is obtained according to the second elastic modulus and the second Poisson's ratio, a total thermal stress strain is obtained according to the second stack displacement, an inelastic thermal stress strain is obtained according to the second initial strain, the second thermal expansion coefficient, the local temperature and the reference temperature, an elastic thermal stress strain is obtained according to the total thermal stress strain and the inelastic thermal stress strain, a third elastic stress is obtained according to the elastic thermal stress strain and the third elastic matrix, and the third inelastic stress is obtained by subtracting the thermal stress from the third elastic stress.

9. The finite element analysis-based stress profile evaluation method of a state of operation of an electrical stack according to claim 8, characterized in that, The third elastic stress and the third inelastic stress are obtained by the following calculation formula in the step S4: ； ； ； ； ； ； ； Wherein, representing said thermal stress; represents the third elastic stress; represents the third inelastic stress; denotes the total strain of the thermal stress; represents the second stack displacement; denotes a transpose matrix operation on the matrix after pre-inversion; denotes the thermal stress inelastic strain; represents the second initial strain; represents the thermal stress elastic strain; represents the running thermal strain; represents the second coefficient of thermal expansion; represents the local temperature; represents the reference temperature; denotes the third elastic matrix; G' represents the second elastic modulus; represents the second Poisson's ratio.

10. The finite element analysis-based stress profile evaluation method of a state of operation of an electrical stack according to claim 1, characterized in that, The step S1 further includes: For the pole connector in the solid oxide stack, the channels and ribs contained in the pole connector are modeled as a porous channel-rib area to obtain a multi-physical field stack CFD model, and the multi-physical field stack CFD model is replaced by the full-size stack model.

Images