Solid oxide battery co-electrolysis numerical simulation method considering mass transfer and heat transfer coupling

Abstract

The invention provides a solid oxide battery co-electrolysis numerical simulation method considering mass transfer and heat transfer coupling, and belongs to the technical field of solid oxide batteries. Calculating an effective mass transfer coefficient field of an electrode porous medium according to the multi-mechanism mass transfer coupling matrix, sequentially solving an electrochemical reaction field, a mass transfer field, a heat transfer field and a stress field by using an operator splitting method, updating each field variable, and triggering adaptive time step adjustment through a multi-physical field coupling anomaly monitoring vector. And calling a cross-scale feedback correction matrix to dynamically correct the mass transfer coefficient according to the updated temperature field and current density field, and when the maximum relative variation of each field variable in three continuous time steps is smaller than a convergence criterion, outputting a steady-state simulation result. The technical problem of mass transfer parameter distortion caused by multi-physical field strong coupling in the solid oxide battery co-electrolysis numerical simulation process is solved.

Description

Technical Field

[0001] This invention belongs to the field of solid oxide battery technology, and more specifically, relates to a numerical simulation method for co-electrolysis of solid oxide batteries that considers mass transfer and heat transfer coupling. Background Technology

[0002] Solid oxide battery co-electrolysis technology is an important approach to achieving hydrogen and fuel production through high-temperature steam and carbon dioxide electrolysis. Its numerical simulation relies on accurate descriptions of the electrochemical reaction field, mass transfer field, and heat transfer field. Traditional numerical simulation methods employ fixed mass transfer coefficients and constant time steps for multiphysics coupling calculations, using preset diffusion coefficients and permeability parameters to describe the gas transport process within porous electrodes. This method is widely used in operating condition simulations such as temperature field distribution calculations, current density field evolution analysis, and gas component concentration prediction for solid oxide battery co-electrolysis systems. In existing technologies, the co-electrolysis process of solid oxide batteries involves high-temperature electrochemical reactions and complex mass and heat transfer phenomena within porous media. Macroscopic temperature changes and current density distributions significantly affect the gas transport characteristics of the microscopic pore structure. Traditional methods treat mass transfer parameters as static constants, neglecting the dynamic feedback effect of macroscopic field variables on microscopic mass transfer mechanisms. This leads to inaccurate descriptions of the coupling relationships between multiple mass transfer mechanisms, such as Knudsen diffusion, molecular diffusion, and viscous flow. Accumulated concentration field calculation deviations occur in key regions such as the interface between the electrode functional layer and the diffusion layer, resulting in inaccurate predictions of electrochemical reaction rates and distorted temperature field distributions. In other words, existing technologies suffer from the technical problem of mass transfer parameter distortion caused by strong multi-physics coupling in the numerical simulation of solid oxide battery co-electrolysis. Summary of the Invention

[0003] In view of this, the present invention provides a numerical simulation method for co-electrolysis of solid oxide batteries that considers mass transfer and heat transfer coupling, which can solve the technical problem of mass transfer parameter distortion caused by strong coupling of multiple physics fields in the numerical simulation of co-electrolysis of solid oxide batteries in the prior art.

[0004] This invention is implemented as follows: It provides a numerical simulation method for co-electrolysis of solid oxide batteries considering mass and heat transfer coupling. A three-dimensional geometric model of the solid oxide battery co-electrolysis system is obtained, and the electrolyte layer, anode functional layer, anode diffusion layer, cathode functional layer, and cathode diffusion layer are meshed. Temperature distribution data, gas component concentration distribution data, and current density distribution data under initial operating conditions are collected, normalized, and stored as initial temperature field vector, initial concentration field vector, and initial current density field vector, respectively. Based on the initial temperature field vector and initial concentration field vector, local pore size distribution parameters and gas state parameters within the porous electrode are calculated. The Knudsen diffusion coefficient, molecular diffusion coefficient, and viscous flow permeability are calculated using a multi-mechanism mass transfer coupling matrix. The effective mass transfer coefficient field of the electrode porous medium is obtained. Electrochemical reaction field, mass transfer field, heat transfer field, and stress field are calculated using the operator splitting method to obtain updated temperature field vector, updated concentration field vector, updated current density field vector, and updated stress field vector. The multiphysics coupling anomaly monitoring vector is calculated, and the adjusted time step and adjusted relaxation factor are obtained through the adaptive time step adjustment function. The effective mass transfer coefficient field of the electrode porous medium is corrected by calling the cross-scale feedback correction matrix to obtain the corrected effective mass transfer coefficient field. When the maximum relative change of the updated temperature field vector, updated concentration field vector, and updated current density field vector within three consecutive time steps is less than the convergence criterion, the steady-state numerical simulation result is output.

[0005] In particular, the porous medium region of the electrode adopts non-uniform mesh refinement treatment. In the interface region between the electrode and the electrolyte and in the region with a large porosity gradient inside the electrode, the mesh size is set to 0.3 to 0.5 times the mesh size of the conventional region.

[0006] In this case, the mesh size in the region of the electrode diffusion layer far from the interface is set to 1.5 to 2 times the mesh size of the regular region, where the mesh size of the regular region is ∈ [10μm, 50μm].

[0007] The normalization process for the initial temperature field vector involves subtracting the system's lowest temperature value from the temperature value of each grid node in the temperature distribution data, and then dividing by the temperature range value, which is the difference between the system's highest and lowest temperatures.

[0008] The normalization process for the initial concentration field vector involves dividing the concentration value of each gas component at each grid node in the gas component concentration distribution data by the concentration value of the gas component at the inlet.

[0009] The normalization process for the initial current density field vector involves dividing the current density value of each grid node in the current density distribution data by the average current density value.

[0010] The multi-mechanism mass transfer coupling matrix is ​​a 3×3 matrix, with diagonal elements ranging from [0.6, 0.9] and off-diagonal elements ranging from [0.1, 0.3]. The sum of all elements in the matrix is ​​normalized to 3.

[0011] The Knudsen diffusion coefficient is calculated based on the local pore size distribution parameters and the mass of gas molecules. The molecular diffusion coefficient is calculated based on the gas temperature and gas pressure values, determined by the Chapman-Ninscog theory. The viscous flow permeability is calculated based on the porosity and tortuosity of the porous medium, determined by the Karman-Kozani equation.

[0012] Among them, the electrochemical reaction field calculation uses the Butler-Wolmer equation to describe the electrode reaction kinetics, and the mass transfer field calculation uses the Stefan-Maxwell diffusion equation combined with Darcy's law to describe the gas transport process.

[0013] Among them, the heat transfer field calculation adopts the energy conservation equation and considers the heat of electrochemical reaction, Joule heat, convective heat transfer, and radiative heat transfer, while the stress field calculation adopts the thermoelastic constitutive equation to describe the mechanical behavior of the material.

[0014] Among them, the multi-physics coupling anomaly monitoring vector contains four elements corresponding to the electrochemical reaction field residual, mass transfer field residual, heat transfer field residual, and stress field residual, respectively. The electrochemical reaction field residual is defined as the maximum relative difference between the updated current density field vector and the initial current density field vector.

[0015] The adaptive time step adjustment function is calculated as follows: when the ratio of the maximum element value of the multiphysics coupling anomaly monitoring vector to the corresponding preset threshold is greater than 1.5, the adjusted time step is the current time step multiplied by 0.5, and the adjusted relaxation factor is the current relaxation factor multiplied by 1.2.

[0016] The cross-scale feedback correction matrix is ​​a 2×3 matrix, with the three elements in the first row being the correction coefficients for temperature on the Knudsen diffusion coefficient, the molecular diffusion coefficient, and the viscous flow permeability, respectively.

[0017] The calculation method for the first row of the cross-scale feedback correction matrix is ​​to normalize the difference between the current grid node temperature and the initial temperature and then multiply it by the temperature sensitivity coefficient. The temperature sensitivity coefficient is 0.5 for the Knudsen diffusion coefficient, 1.5 for the molecular diffusion coefficient, and 0.3 for the viscous flow permeability.

[0018] The calculation method for the second row of the cross-scale feedback correction matrix is ​​to normalize the difference between the current current density of the current grid node and the average current density value and then multiply it by the current density sensitivity coefficient. The current density sensitivity coefficient is 0.2 for the Knudsen diffusion coefficient, 0.4 for the molecular diffusion coefficient, and 0.6 for the viscous flow permeability.

[0019] The correction of the effective mass transfer coefficient field is achieved by multiplying the Knudsen diffusion coefficient, molecular diffusion coefficient, and viscous flow permeability in the effective mass transfer coefficient field of the electrode porous medium by the sum of the corresponding column elements of the cross-scale feedback correction matrix and adding 1.

[0020] This invention constructs a cross-scale feedback correction matrix to map changes in the macroscopic temperature and current density fields to the dynamic adjustment of microscopic mass transfer parameters in real time. It also combines a multi-mechanism mass transfer coupling matrix to describe the coupling relationships of three mass transfer mechanisms: Knudsen diffusion, molecular diffusion, and viscous flow. Within each time step, the effective mass transfer coefficient field of the porous electrode medium is corrected based on the updated temperature and current density field vectors, enabling the mass transfer parameters to respond to the evolution of macroscopic field variables and avoiding the concentration field calculation bias caused by static mass transfer coefficients in traditional methods. Furthermore, this invention decomposes the strongly coupled multiphysics problem into weakly coupled subproblems using an operator splitting method, solving them sequentially. A multiphysics coupling anomaly monitoring vector is introduced to monitor the residuals of each subproblem in real time. When the residuals exceed a preset threshold, an adaptive time step adjustment function dynamically modifies the calculation step size and relaxation factor, ensuring the convergence and stability of the electrochemical reaction field, mass transfer field, heat transfer field, and stress field calculations. This eliminates the numerical oscillations and divergences that easily occur in traditional fixed-time-step methods under strongly coupled conditions. In summary, this invention solves the technical problem of mass transfer parameter distortion caused by strong coupling of multiple physics fields in the numerical simulation of co-electrolysis of solid oxide batteries by establishing a dynamic feedback mechanism between macroscopic field variables and microscopic mass transfer parameters, and by combining it with an adaptive numerical solution strategy. Attached Figure Description

[0021] Figure 1 This is a flowchart of the method of the present invention.

[0022] Figure 2 The graph shows the relationship between the Knudsen diffusion coefficient and the local pore size within the anode functional layer in this embodiment.

[0023] Figure 3 The diagram shows the distribution of the Joule thermal power density of the electrolyte layer as a function of current density in the embodiment. Detailed Implementation

[0024] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings.

[0025] like Figure 1 The diagram shows a flowchart of a numerical simulation method for co-electrolysis of solid oxide batteries considering mass and heat transfer coupling, provided by this invention. This method includes the following steps:

[0026] S01. Obtain the established three-dimensional geometric model of the solid oxide battery co-electrolysis system, and perform meshing on the electrolyte layer, anode functional layer, anode diffusion layer, cathode functional layer, and cathode diffusion layer. The porous dielectric region of the electrode is treated with non-uniform mesh refinement.

[0027] S02. Collect temperature distribution data, gas component concentration distribution data, and current density distribution data of the solid oxide battery co-electrolysis system under initial operating conditions. Normalize the temperature distribution data and store it as an initial temperature field vector. Normalize the gas component concentration distribution data and store it as an initial concentration field vector. Normalize the current density distribution data and store it as an initial current density field vector.

[0028] S03. Calculate the local pore size distribution parameters and gas state parameters in the porous electrode based on the initial temperature field vector and the initial concentration field vector. Call the multi-mechanism mass transfer coupling matrix to calculate the Knudsen diffusion coefficient, molecular diffusion coefficient, and viscous flow permeability to obtain the effective mass transfer coefficient field of the porous medium of the electrode.

[0029] S04. Based on the initial temperature field vector, initial concentration field vector, and initial current density field vector, perform electrochemical reaction field calculation, mass transfer field calculation, heat transfer field calculation, and stress field calculation using the operator splitting method to obtain the updated temperature field vector, updated concentration field vector, updated current density field vector, and updated stress field vector.

[0030] S05. Calculate the multiphysics coupling anomaly monitoring vector based on the differences between the updated temperature field vector and the initial temperature field vector, the updated concentration field vector and the initial concentration field vector, the updated current density field vector and the initial current density field vector, and the updated stress field vector and the initial stress field vector. When an element in the multiphysics coupling anomaly monitoring vector exceeds the corresponding preset threshold, obtain the adjusted time step and the adjusted relaxation factor through an adaptive time step adjustment function.

[0031] S06. Calculate the electrochemical reaction heat power density and Joule heat power density of each grid node based on the updated current density field vector and the updated temperature field vector. Assign the updated temperature field vector to the initial temperature field vector. Call the cross-scale feedback correction matrix to correct the effective mass transfer coefficient field of the electrode porous medium to obtain the corrected effective mass transfer coefficient field.

[0032] S07. Calculate the maximum relative changes of the updated temperature field vector, updated concentration field vector, and updated current density field vector within three consecutive time steps. When the maximum relative changes are all less than the convergence criterion, output the steady-state numerical simulation results of the solid oxide battery co-electrolysis system. When the maximum relative changes are greater than or equal to the convergence criterion, assign the updated concentration field vector to the initial concentration field vector, assign the updated current density field vector to the initial current density field vector, and return to step S04 to continue iterative solution.

[0033] The non-uniform mesh refinement treatment of the porous medium region of the electrode refers to setting the mesh size to 0.3 to 0.5 times that of the conventional region mesh size in the interface region between the electrode and the electrolyte and in the region with a large porosity gradient inside the electrode, and setting the mesh size to 1.5 to 2 times that of the conventional region mesh size in the region of the electrode diffusion layer far from the interface, wherein the conventional region mesh size ∈ [10μm, 50μm].

[0034] The anode functional layer refers to a porous anode electrode layer with electrochemical catalytic activity, having a thickness ∈ [5 μm, 15 μm] and a porosity ∈ [30%, 40%]. The anode diffusion layer refers to a porous anode support layer for gas transport, having a thickness ∈ [300 μm, 600 μm] and a porosity ∈ [40%, 50%]. The cathode functional layer refers to a porous cathode electrode layer with electrochemical catalytic activity, having a thickness ∈ [5 μm, 15 μm] and a porosity ∈ [30%, 40%]. The cathode diffusion layer refers to a porous cathode support layer for gas transport, having a thickness ∈ [300 μm, 600 μm] and a porosity ∈ [40%, 50%].

[0035] The normalization of the initial temperature field vector involves subtracting the system's lowest temperature value from the temperature value of each grid node in the temperature distribution data, and then dividing by the temperature range value, where the temperature range value is the difference between the system's highest and lowest temperatures. The normalization of the initial concentration field vector involves dividing the concentration value of each gas component in each grid node of the gas component concentration distribution data by the concentration value of that gas component at the inlet. The normalization of the initial current density field vector involves dividing the current density value of each grid node of the current density distribution data by the average current density value.

[0036] The local pore size distribution parameters refer to the average pore size value and the standard deviation of pore size distribution for each grid node in the porous electrode medium. The gas state parameters refer to the temperature, pressure, density, and viscosity values ​​of the gas at each grid node.

[0037] The multi-mechanism mass transfer coupling matrix is ​​used to describe the coupling relationship of three mass transfer mechanisms—Knudsen diffusion, molecular diffusion, and viscous flow—within the porous electrode medium. The multi-mechanism mass transfer coupling matrix is ​​a 3×3 matrix. The elements in the first row and first column are the coupling coefficients of Knudsen diffusion with itself, the elements in the first row and second column are the coupling coefficients of Knudsen diffusion with molecular diffusion, and the elements in the first row and third column are the coupling coefficients of Knudsen diffusion with viscous flow. The elements in the second row and first column are the coupling coefficients of molecular diffusion with Knudsen diffusion, the elements in the second row and second column are the coupling coefficients of molecular diffusion with itself, the elements in the second row and third column are the coupling coefficients of molecular diffusion with viscous flow, the elements in the third row and first column are the coupling coefficients of viscous flow with Knudsen diffusion, the elements in the third row and second column are the coupling coefficients of viscous flow with molecular diffusion, and the elements in the third row and third column are the coupling coefficients of viscous flow with itself.

[0038] The diagonal elements of the multi-mechanism mass transfer coupling matrix range from [0.6, 0.9], and the off-diagonal elements range from [0.1, 0.3]. The sum of all elements in the matrix is ​​normalized to 3. The Knudsen diffusion coefficient is calculated based on the local pore size distribution parameters and the gas molecule mass; Knudsen diffusion dominates when the pore size is smaller than the mean free path of the gas molecules. The molecular diffusion coefficient is calculated based on the gas temperature and pressure values, determined using the Chapman-Nskoglund theory. The viscous flow permeability is calculated based on the porosity and tortuosity of the porous medium, determined using the Karman-Kozani equation. The effective mass transfer coefficient field of the electrode porous medium includes the Knudsen diffusion coefficient, molecular diffusion coefficient, and viscous flow permeability for each grid node.

[0039] The operator splitting method refers to decomposing a strongly coupled multiphysics problem into multiple weakly coupled subproblems, solving each subproblem sequentially within each time step, and updating the field variables through iteration. The electrochemical reaction field calculation refers to solving for the local current density distribution and electrochemical reaction rate distribution on the electrode surface, using the Butler-Wolmer equation to describe the electrode reaction kinetics. The mass transfer field calculation refers to solving for the concentration distribution of gas components within the porous electrode, using the Stefan-Maxwell diffusion equation combined with Darcy's law to describe the gas transport process. The heat transfer field calculation refers to solving for the temperature distribution of the solid oxide battery co-electrolysis system, using the energy conservation equation and considering electrochemical reaction heat, Joule heat, convective heat transfer, and radiative heat transfer. The stress field calculation refers to solving for the thermal stress distribution of each layer of the battery material, using the thermoelastic constitutive equation to describe the mechanical behavior of the material.

[0040] The initial stress field vector refers to the vector formed after normalizing the stress values ​​of each grid node of each layer of battery material under the initial working condition. The initial stress field vector is obtained by normalizing the initial working condition stress distribution data when step S05 is executed for the first time.

[0041] The multiphysics coupling anomaly monitoring vector is used to monitor numerical anomalies and convergence problems in the solution process of each sub-problem in real time. The multiphysics coupling anomaly monitoring vector contains four elements corresponding to the residuals of the electrochemical reaction field, mass transfer field, heat transfer field, and stress field, respectively. The residual of the electrochemical reaction field is defined as the maximum relative difference between the updated current density field vector and the initial current density field vector; the residual of the mass transfer field is defined as the maximum relative difference between the updated concentration field vector and the initial concentration field vector; the residual of the heat transfer field is defined as the maximum relative difference between the updated temperature field vector and the initial temperature field vector; and the residual of the stress field is defined as the maximum relative difference between the updated stress field vector and the initial stress field vector.

[0042] The corresponding preset thresholds are set to 0.05 for the electrochemical reaction field residual, 0.08 for the mass transfer field residual, 0.06 for the heat transfer field residual, and 0.10 for the stress field residual.

[0043] The adaptive time step adjustment function is used to dynamically adjust the time step and relaxation factor of the numerical solution based on the multiphysics coupling anomaly monitoring vector to ensure computational convergence. The inputs include the maximum element value of the multiphysics coupling anomaly monitoring vector, the current time step, and the current relaxation factor. The outputs are the adjusted time step and the adjusted relaxation factor. The calculation logic of the adaptive time step adjustment function is as follows: when the ratio of the maximum element value of the multiphysics coupling anomaly monitoring vector divided by the corresponding preset threshold is >1.5, the adjusted time step is the current time step multiplied by 0.5, and the adjusted relaxation factor is the current relaxation factor multiplied by 1.2; when the ratio ∈ [1, 1.5], the adjusted time step is the current time step multiplied by 0.8, and the adjusted relaxation factor is the current relaxation factor multiplied by 1.1; when the ratio <0.5, the adjusted time step is the current time step multiplied by 1.5, and the adjusted relaxation factor is the current relaxation factor multiplied by 0.9; when the ratio ∈ [0.5, 1), the adjusted time step is the current time step multiplied by 1.0, and the adjusted relaxation factor is the current relaxation factor multiplied by 1.0. The relaxation factor has a value range of [0.3, 0.9] and is used to control the weighted mixing ratio of new and old field variables during the iterative solution process.

[0044] The electrochemical reaction heat power density refers to the heat power generated or absorbed by the electrochemical reaction within a unit volume of electrode material. Its calculation is based on the local current density value, reaction overpotential, and reaction entropy change in the updated current density field vector. The Joule heat power density refers to the Joule heat power generated by the current passing through a unit volume of material. Its calculation is based on the local current density value and material resistivity in the updated current density field vector.

[0045] The cross-scale feedback correction matrix is ​​used to dynamically correct the mass transfer parameters of the microporous structure according to the changes in the macroscopic temperature field and current density field. The cross-scale feedback correction matrix is ​​a 2×3 matrix. The three elements in the first row are the correction coefficients of temperature on the Knudsen diffusion coefficient, temperature on the molecular diffusion coefficient, and temperature on the viscous flow permeability, respectively. The three elements in the second row are the correction coefficients of current density on the Knudsen diffusion coefficient, current density on the molecular diffusion coefficient, and current density on the viscous flow permeability, respectively.

[0046] The first row of the cross-scale feedback correction matrix is ​​calculated by normalizing the difference between the current grid node temperature and the initial temperature, and then multiplying it by a temperature sensitivity coefficient. This temperature sensitivity coefficient is set to 0.5 for the Knudsen diffusion coefficient, 1.5 for the molecular diffusion coefficient, and 0.3 for viscous flow permeability. The second row of the cross-scale feedback correction matrix is ​​calculated by normalizing the difference between the current density of the current grid node and the average current density, and then multiplying it by a current density sensitivity coefficient. This current density sensitivity coefficient is set to 0.2 for the Knudsen diffusion coefficient, 0.4 for the molecular diffusion coefficient, and 0.6 for viscous flow permeability.

[0047] The correction of the effective mass transfer coefficient field refers to multiplying the Knudsen diffusion coefficient, molecular diffusion coefficient, and viscous flow permeability in the effective mass transfer coefficient field of the electrode porous medium by the sum of the corresponding column elements of the cross-scale feedback correction matrix plus 1. The corrected effective mass transfer coefficient field is used in the mass transfer field calculation of step S04 to update the effective diffusion coefficient and permeability of the gas components.

[0048] The convergence criteria are set to 0.01 for the updated temperature field vector, 0.02 for the updated concentration field vector, and 0.015 for the updated current density field vector, representing the upper limit of the maximum relative change of the field variables within three consecutive time steps. The steady-state numerical simulation results include the temperature values, gas component concentration values, current density values, stress values, and overall electrolysis voltage, electrolysis current, gas production rate, and energy conversion efficiency of each grid node in the solid oxide battery co-electrolysis system.

[0049] As an optional implementation, the present invention also provides a method for forming a solid oxide battery co-electrolysis numerical simulation system by means of a computer, wherein the computer is provided with a readable storage medium, the readable storage medium stores program instructions, and the program instructions execute the above-described solid oxide battery co-electrolysis numerical simulation method when the computer is run.

[0050] The specific implementation methods of the above steps are described in detail below.

[0051] The specific implementation of step S01 involves obtaining a three-dimensional geometric model of the established solid oxide battery co-electrolysis system. This model includes five main structural layers: an electrolyte layer, an anode functional layer, an anode diffusion layer, a cathode functional layer, and a cathode diffusion layer. The anode functional layer has a thickness of 5–15 μm and a porosity of 30%–40%; the anode diffusion layer has a thickness of 300–600 μm and a porosity of 40%–50%; the cathode functional layer has a thickness of 5–15 μm and a porosity of 30%–40%; and the cathode diffusion layer has a thickness of 300–600 μm and a porosity of 40%–50%. Then, each layer is meshed. The discretized computational mesh is generated using the finite volume method or the finite element method. For the porous medium region of the electrode, a non-uniform mesh refinement strategy is adopted. In the interface region between the electrode and the electrolyte, and in the region with a large porosity gradient inside the electrode, the mesh size is set to 0.3 to 0.5 times the mesh size of the normal region to improve the computational accuracy at the interface. In the region of the electrode diffusion layer far from the interface, the mesh size is set to 1.5 to 2 times the mesh size of the normal region to reduce the computational load. The mesh size of the normal region is selected between 10 and 50 μm. This mesh refinement strategy is based on the principle of multi-scale analysis and can optimize the allocation of computational resources while ensuring computational accuracy.

[0052] The specific implementation of step S02 involves collecting physical field distribution data of the solid oxide battery co-electrolysis system under initial operating conditions, including temperature values, concentration values ​​of each gas component, and current density values ​​at each grid node. These raw data are then normalized to eliminate the influence of different physical quantities. For temperature distribution data, the normalized temperature is obtained by subtracting the system's lowest temperature value from the temperature value of each grid node and dividing by the temperature range value. The temperature range value is defined as the difference between the system's highest and lowest temperatures. The normalized temperature values ​​are stored as an initial temperature field vector in the order of the grid nodes. For gas component concentration distribution data, the normalized concentration is obtained by dividing the concentration value of each gas component at each grid node by the concentration value of that gas component at the inlet. The normalized concentration values ​​are stored as an initial concentration field vector in the order of the grid nodes. For current density distribution data, the normalized current density is obtained by dividing the current density value of each grid node by the average current density value. The normalized current density values ​​are stored as an initial current density field vector in the order of the grid nodes. This normalization process, based on the principle of numerical scaling, can improve the stability and convergence of subsequent numerical calculations.

[0053] The specific implementation of step S03 involves calculating the local pore size distribution parameters and gas state parameters of each grid node within the porous electrode based on the initial temperature field vector and the initial concentration field vector. The local pore size distribution parameters include the average pore size value and the standard deviation of the pore size distribution. The gas state parameters include temperature, pressure, density, and viscosity. These parameters are calculated using the equation of state and empirical correlations. Then, the multi-mechanism mass transfer coupling matrix is ​​used to calculate the Knudsen diffusion coefficient, molecular diffusion coefficient, and viscous flow permeability. The multi-mechanism mass transfer coupling matrix is ​​a 3×3 matrix used to describe the coupling relationship of the three mass transfer mechanisms within the porous medium of the electrode. The diagonal elements of the matrix range from 0.6 to 0.9, representing the relationship between each mechanism and itself. The coupling strength, with off-diagonal elements ranging from 0.1 to 0.3, represents the cross-coupling strength between different mechanisms. The sum of all elements in the matrix is ​​normalized to 3 to ensure energy conservation. The Knudsen diffusion coefficient is calculated based on local pore size distribution parameters and gas molecule mass. When the pore size is smaller than the mean free path of gas molecules, Knudsen diffusion dominates. The molecular diffusion coefficient is calculated based on gas temperature and pressure values ​​using the Chapman-Nskoglund theory. The viscous flow permeability is calculated based on the porosity and tortuosity of the porous medium using the Karman-Kozani equation. Finally, the effective mass transfer coefficient field of the electrode porous medium is obtained, including the Knudsen diffusion coefficient, molecular diffusion coefficient, and viscous flow permeability of each grid node.

[0054] The specific implementation of step S04 involves using the operator splitting method to decompose the strongly coupled multiphysics problem into four weakly coupled sub-problems: electrochemical reaction field, mass transfer field, heat transfer field, and stress field. Each sub-problem is solved sequentially within each time step. The electrochemical reaction field calculation uses the Butler-Wolmer equation to describe the electrode reaction kinetics. Input parameters include the initial temperature field vector, the initial concentration field vector, and electrode reaction kinetic parameters. Output parameters are the updated current density field vector and the electrochemical reaction rate distribution. The mass transfer field calculation uses the Stefan-Maxwell diffusion equation combined with Darcy's law to describe the gas transport process within the porous electrode. Input parameters include the initial concentration field vector... The input parameters include the initial temperature field vector, the effective mass transfer coefficient field, and the updated concentration field vector. The heat transfer field calculation adopts the energy conservation equation and considers the electrochemical reaction heat, Joule heat, convective heat transfer, and radiative heat transfer. The output parameter is the updated temperature field vector. The stress field calculation adopts the thermoelastic constitutive equation to describe the mechanical behavior of the material. The input parameters include the updated temperature field vector and the thermal expansion coefficient of the material. The output parameter is the updated stress field vector. The operator splitting method is based on the principle of physical field decoupling, which can reduce the solution complexity of multi-physics coupling problems and improve computational efficiency.

[0055] The specific implementation of step S05 involves calculating a multiphysics coupling anomaly monitoring vector to monitor numerical anomalies and convergence issues in the solution process of each sub-problem in real time. This vector contains four elements corresponding to the residuals of the electrochemical reaction field, mass transfer field, heat transfer field, and stress field, respectively. The residual of the electrochemical reaction field is defined as the maximum relative difference between the updated current density field vector and the initial current density field vector; the residual of the mass transfer field is defined as the maximum relative difference between the updated concentration field vector and the initial concentration field vector; the residual of the heat transfer field is defined as the maximum relative difference between the updated temperature field vector and the initial temperature field vector; and the residual of the stress field is defined as the maximum relative difference between the updated stress field vector and the initial stress field vector. These residual values ​​are then compared with corresponding preset thresholds. The preset thresholds are set to 0.05 for the electrochemical reaction field residual, 0.08 for the mass transfer field residual, 0.06 for the heat transfer field residual, and 0.10 for the stress field residual. When any element in the multiphysics coupling anomaly monitoring vector exceeds the corresponding preset threshold, an adaptive time step adjustment is invoked. The function takes the maximum element value of the multiphysics coupling anomaly monitoring vector, the current time step, and the current relaxation factor as input parameters. Its output parameters are the adjusted time step and the adjusted relaxation factor. The adjustment logic is to calculate the ratio of the maximum element value to a corresponding preset threshold. When the ratio is greater than 1.5, the adjusted time step is the current time step multiplied by 0.5, and the adjusted relaxation factor is the current relaxation factor multiplied by 1.2. When the ratio is between 1 and 1.5, the adjusted time step is the current time step multiplied by 0.8, and the adjusted relaxation factor is... The relaxation factor is the current relaxation factor multiplied by 1.1. When the ratio is less than 0.5, the adjusted time step is the current time step multiplied by 1.5 and the adjusted relaxation factor is the current relaxation factor multiplied by 0.9. When the ratio is between 0.5 and 1, the adjusted time step and the adjusted relaxation factor remain unchanged. The relaxation factor ranges from 0.3 to 0.9 and is used to control the weighted mixing ratio of the old and new field variables during the iterative solution process. This adaptive adjustment strategy is based on the error feedback control principle and can dynamically optimize the numerical solution parameters according to the convergence of the calculation.

[0056] The specific implementation of step S06 involves calculating the electrochemical reaction heat power density and Joule heat power density of each grid node based on the updated current density field vector and the updated temperature field vector. The calculation of the electrochemical reaction heat power density is based on the local current density value, reaction overpotential, and reaction entropy change, while the calculation of the Joule heat power density is based on the local current density value and material resistivity. Then, the updated temperature field vector is assigned to the initial temperature field vector for the next iteration. Simultaneously, a cross-scale feedback correction matrix is ​​called to correct the effective mass transfer coefficient field of the porous electrode medium. The cross-scale feedback correction matrix is ​​a 2×3 matrix used to dynamically correct the mass transfer parameters of the microporous structure based on changes in the macroscopic temperature field and current density field. The first row of three elements represents the correction coefficients for temperature on the Knudsen diffusion coefficient, temperature on the molecular diffusion coefficient, and temperature on the viscous flow permeability, respectively. The second row of three elements represents the correction coefficients for current density on the Knudsen diffusion coefficient, temperature on the molecular diffusion coefficient, and temperature on the viscous flow permeability, respectively. The correction coefficients for the diffusion coefficient and the current density correction coefficient for viscous flow permeability are calculated as follows: The first row of elements is calculated by normalizing the difference between the current grid node temperature and the initial temperature and then multiplying it by a temperature sensitivity coefficient. The temperature sensitivity coefficient is 0.5 for the Knudsen diffusion coefficient, 1.5 for the molecular diffusion coefficient, and 0.3 for the viscous flow permeability. The second row of elements is calculated by normalizing the difference between the current density of the current grid node and the average current density value and then multiplying it by a current density sensitivity coefficient. The current density sensitivity coefficient is 0.2 for the Knudsen diffusion coefficient, 0.4 for the molecular diffusion coefficient, and 0.6 for the viscous flow permeability. The correction of the effective mass transfer coefficient field refers to multiplying the Knudsen diffusion coefficient, molecular diffusion coefficient, and viscous flow permeability in the effective mass transfer coefficient field by the sum of the corresponding column elements of the cross-scale feedback correction matrix and adding 1. This cross-scale correction strategy is based on multi-scale coupling theory and can capture the influence of macroscopic physical field changes on microscopic mass transfer characteristics.

[0057] The specific implementation of step S07 involves calculating the maximum relative changes of the updated temperature field vector, updated concentration field vector, and updated current density field vector over three consecutive time steps. These maximum relative changes are then compared with the corresponding convergence criteria. The convergence criteria are set to 0.01 for the updated temperature field vector, 0.02 for the updated concentration field vector, and 0.015 for the updated current density field vector. When all maximum relative changes are less than the convergence criteria, it indicates that the system has reached steady state. At this point, the steady-state numerical simulation results of the solid oxide battery co-electrolysis system are output. This includes the temperature values, gas component concentration values, current density values, stress values ​​of each grid node, as well as the overall electrolysis voltage, electrolysis current, gas production rate, and energy conversion efficiency of the system. When the maximum relative change is greater than or equal to the convergence criterion, it indicates that the system has not yet reached a steady state. At this time, the updated concentration field vector is assigned to the initial concentration field vector, and the updated current density field vector is assigned to the initial current density field vector. Then, the process returns to step S04 to continue the next round of iteration. This convergence determination strategy is based on the trend of field variable changes over a continuous time step, which can accurately identify whether the system has reached a steady state and avoid false convergence.

[0058] It should be noted that the key technical ideas of this invention include three core technologies: a multi-mechanism mass transfer coupling matrix, an adaptive time step adjustment function, and a cross-scale feedback correction matrix. The multi-mechanism mass transfer coupling matrix, by introducing cross-coupling coefficients for Knudsen diffusion, molecular diffusion, and viscous flow, overcomes the limitation of independent processing of each mass transfer mechanism in traditional methods. It can accurately describe the complex multi-mechanism synergistic mass transfer process within porous electrodes, significantly improving the prediction accuracy of the gas concentration field. The adaptive time step adjustment function adjusts the calculation parameters in real time based on the multi-physics coupling anomaly monitoring vector. Compared with the traditional fixed time step method, it can significantly improve the convergence speed and stability of numerical solutions while ensuring computational accuracy, effectively avoiding numerical oscillations and divergence problems in the solution process of strongly coupled systems. The cross-scale feedback correction matrix establishes a dynamic correlation between macroscopic physical field changes and microscopic mass transfer parameters, overcoming the deficiency of traditional single-scale simulation methods in capturing inter-scale interactions, and realizing real-time feedback correction from macroscopic temperature and current density fields to the mass transfer characteristics of microscopic porous structures. The synergistic effect of these three technical approaches forms a complete multi-physics, multi-scale coupled solution system. It achieves accurate mass transfer description at the microscale through a multi-mechanism mass transfer coupling matrix, efficient and stable solution at the macroscale through an adaptive time step adjustment function, and bidirectional coupling between different scales through a cross-scale feedback correction matrix. Compared with traditional single-scale, single-mechanism simulation methods, it can comprehensively improve the accuracy, efficiency, and reliability of numerical simulation of solid oxide battery co-electrolysis systems.

[0059] It should be noted that the electrode structure parameters include the thickness of the anode and cathode functional layers (5–15 μm, porosity 30%–40%), and the thickness of the anode and cathode diffusion layers (300–600 μm, porosity 40%–50%). These parameters were determined by observing the microstructure of actual solid oxide battery samples using scanning electron microscopy, combined with image processing techniques and statistical analysis of the structural characteristics of a large number of samples. The mesh size parameters include a mesh size of 10–50 μm in the conventional region, a mesh size of 0.3–0.5 times that of the conventional region in the interface-refined region, and a mesh size of 1.5–2 times that of the conventional region in the region far from the interface. These parameters were obtained through mesh independence verification experiments, i.e., numerical simulations were performed at different mesh sizes, and the calculation results were compared to select the optimal mesh size range that balances computational accuracy and efficiency. The diagonal elements of the multi-mechanism mass transfer coupling matrix range from 0.6 to 0.9, while the off-diagonal elements range from 0.1 to 0.3. These parameters were obtained by fitting and optimizing numerical simulation results with gas diffusion experimental data. The gas diffusion experiments used the concentration gradient method to measure the effective diffusion coefficient within the porous electrode under different temperature and pressure conditions. The optimal value range of the coupling matrix elements was then determined using the least squares method. Temperature sensitivity coefficients were set to 0.5 for the Knudsen diffusion coefficient, 1.5 for the molecular diffusion coefficient, and 0.3 for viscous flow permeability. Current density sensitivity coefficients were set to 0.2 for the Knudsen diffusion coefficient, 0.4 for the molecular diffusion coefficient, and 0.6 for viscous flow permeability. These parameters were determined through battery performance testing experiments under different temperature and current density conditions, collecting voltage-current characteristic curves and impedance spectrum data under each condition. The values ​​of each sensitivity coefficient were determined through multiple regression analysis and parameter sensitivity analysis of the experimental data. The preset thresholds include an electrochemical reaction field residual threshold of 0.05, a mass transfer field residual threshold of 0.08, a heat transfer field residual threshold of 0.06, and a stress field residual threshold of 0.10. These thresholds were determined through numerical convergence analysis of numerous examples. The convergence speed, computational accuracy, and numerical stability under different threshold settings were statistically analyzed, and the threshold combination with the best overall performance was selected. The ratio segmentation thresholds of 1.5, 1.0, and 0.5 in the adaptive time step adjustment function, along with the corresponding time step adjustment coefficients of 0.5, 0.8, 1.0, and 1.5, and relaxation factor adjustment coefficients of 1.2, 1.1, 1.0, and 0.9, were determined through a trial-and-error method combined with an automatic optimization algorithm in numerical experiments. Extensive numerical simulation tests were conducted under typical operating conditions, and the parameter combination that achieves the fastest convergence speed and best numerical stability was searched using a genetic algorithm or particle swarm optimization algorithm. The convergence criteria are set to 0.01 for the temperature field, 0.02 for the concentration field, and 0.015 for the current density field. These criteria are determined by comparing the numerical simulation results with the experimental measurement data through error analysis, and the minimum criterion value is selected to keep the relative error between the simulation results and the experimental data within the acceptable range for engineering.The relaxation factor range of 0.3 to 0.9 was determined through iterative convergence theory analysis and numerical stability tests. Too small a relaxation factor leads to slow convergence, while too large a relaxation factor leads to numerical oscillation. The safe range of the relaxation factor to ensure convergence was determined through tests under different operating conditions.

[0060] It should be noted that this invention also solves the following technical problem: computational divergence caused by strong multiphysics coupling in numerical simulation of solid oxide battery co-electrolysis. Traditional numerical simulation methods use a fixed time step for multiphysics coupling iterative solutions. When the coupling strength between the electrochemical reaction field, mass transfer field, heat transfer field, and stress field is large, the solution error of a certain subproblem will be amplified during the iteration process and propagate to other subproblems, ultimately leading to overall computational divergence or non-physical numerical oscillations. This invention constructs a multiphysics coupling anomaly monitoring vector to monitor the residuals of the electrochemical reaction field, mass transfer field, heat transfer field, and stress field in real time. When any residual exceeds the corresponding preset threshold, the adaptive time step adjustment function dynamically reduces the time step and increases the relaxation factor according to the ratio of the residual to the threshold. This reduces the update amplitude of field variables in a single iteration, keeping the solution error of each subproblem within an acceptable range. This ensures stable convergence of numerical calculations under strong coupling conditions and avoids the computational failure problem that easily occurs in traditional fixed-step methods under high current density conditions or large temperature gradient regions.

[0061] Furthermore, this invention addresses the technical problem of inaccurate description of multi-scale mass transfer mechanisms within porous electrode media in numerical simulations of solid oxide battery co-electrolysis. In the porous structures of the electrode functional layer and diffusion layer, gas transport involves three mechanisms simultaneously: Knudsen diffusion, molecular diffusion, and viscous flow. The dominant role of these three mechanisms varies with pore size, gas pressure, and temperature conditions. Traditional methods often simplify this to a single diffusion mechanism or use empirical formulas to superimpose the contributions of each mechanism, failing to accurately reflect the coupling relationships between different mass transfer mechanisms. This invention establishes a coupling description framework for the three mass transfer mechanisms through a multi-mechanism mass transfer coupling matrix. The diagonal elements of the matrix characterize the individual contribution intensity of each mechanism, while the off-diagonal elements quantify the interactions between different mechanisms. By combining local pore size distribution parameters and gas state parameters, the Knudsen diffusion coefficient, molecular diffusion coefficient, and viscous flow permeability are calculated respectively. This ensures that Knudsen diffusion dominates in the small-pore region at the electrode-electrolyte interface, while molecular diffusion and viscous flow play a greater role in the large-pore region of the diffusion layer. This achieves a refined description of the complex mass transfer process within the porous electrode media and improves the accuracy of gas component concentration distribution prediction.

[0062] Specifically, the principle of this invention is as follows: The fundamental reason why this invention can solve the problem of mass transfer parameter distortion lies in establishing a dynamic correlation mechanism between the macroscopic physical field and the microscopic mass transfer mechanism. During the co-electrolysis of solid oxide batteries, increased temperature enhances the thermal motion of gas molecules, leading to an increase in the molecular diffusion coefficient. Simultaneously, it alters the gas mean free path, affecting Knudsen diffusion characteristics. Increased current density intensifies electrochemical reactions, consuming and altering the local gas composition, thus affecting the diffusion driving force. The influence of these macroscopic field variables on microscopic mass transfer is neglected in traditional methods. This invention quantifies the sensitivity coefficients of temperature and current density to three mass transfer mechanisms through a cross-scale feedback correction matrix. It transforms the deviation of the current temperature from the initial temperature and the deviation of the current current density from the average value into correction amounts for the mass transfer coefficient, achieving real-time updates of mass transfer parameters as the macroscopic field evolves. The multi-mechanism mass transfer coupling matrix describes the self-contribution and mutual coupling of each mass transfer mechanism through diagonal and off-diagonal elements, ensuring a reasonable weight allocation for Knudsen diffusion, molecular diffusion, and viscous flow under different pore sizes and pressure conditions. This dynamic correction mechanism, combined with the weakly coupled solution strategy of the operator splitting method, ensures that the mass transfer parameters used in solving each sub-problem within each time step reflect the latest temperature and current density field states, avoiding cumulative errors caused by parameter lag. The adaptive time step adjustment function determines whether numerical anomalies occur in the solution process based on the multi-physics coupling anomaly monitoring vector. By reducing the time step and increasing the relaxation factor, it suppresses numerical oscillations, ensuring the consistency between the mass transfer field calculation and other physics field calculations under strong coupling conditions. This achieves accurate simulation of complex mass and heat transfer coupling phenomena in the solid oxide battery co-electrolysis system.

[0063] The following provides a specific embodiment 1 of the present invention, and the specific implementation of each step in this embodiment 1 is described in detail below.

[0064] In this embodiment, the specific implementation of step S01 is the same as described above, and will not be repeated in detail here.

[0065] The specific implementation of step S02 involves collecting the physical field distribution data of the solid oxide battery co-electrolysis system under initial operating conditions, and then normalizing these raw data. The normalization formula for the initial temperature field vector is expressed as follows:

[0066] ;

[0067] In the formula, For the first Normalized temperature values ​​of each grid node, dimensionless; For the first The actual temperature value of each grid node, in units of ; This is the lowest system temperature value, in units of... ; This is the highest system temperature value, in units of... The normalized temperature values ​​are stored in the order of the grid nodes to form the initial temperature field vector. ,in The total number of grid nodes is dimensionless. The normalization formula for the initial concentration field vector is as follows:

[0068] ;

[0069] In the formula, For the first The first grid node The normalized concentration values ​​of the gas components are dimensionless. For the first The first grid node The actual concentration values ​​of the gas components, in units of ; For the first The concentration values ​​of the gas components at the inlet, in units of ; This is a dimensionless index for gas components. Normalized concentration values ​​are stored sequentially by grid node to form the initial concentration field vector. ,in The number of gas component types is dimensionless. The normalization formula for the initial current density field vector is as follows:

[0070] ;

[0071] In the formula, For the first The normalized current density values ​​of each grid node, dimensionless; For the first The actual current density value of each grid node, in units of ; This is the average current density value, in units of... The calculation method is as follows ,in The index variable is dimensionless and is used for summation.

[0072] The specific implementation of step S03 involves calculating the mass transfer parameters of each grid node within the porous electrode based on the initial temperature field vector and the initial concentration field vector. The formula for the multi-mechanism mass transfer coupling matrix is ​​as follows:

[0073] ;

[0074] In the formula, This is a multi-mechanism mass transfer coupling matrix, dimensionless; The coupling coefficient between Knudsen diffusion and itself, ranging from 0.6 to 0.9, is dimensionless; The coupling coefficient between molecular diffusion and itself, ranging from 0.6 to 0.9, is dimensionless. The coupling coefficient between viscous flow and itself, ranging from 0.6 to 0.9, is dimensionless. The coupling coefficient between Knudsen diffusion and molecular diffusion ranges from 0.1 to 0.3 and is dimensionless. is the coupling coefficient between Knudsen diffusion and viscous flow, with a value ranging from 0.1 to 0.3, and is dimensionless; The coupling coefficient between molecular diffusion and Knudsen diffusion ranges from 0.1 to 0.3 and is dimensionless. is the coupling coefficient between molecular diffusion and viscous flow, with a value ranging from 0.1 to 0.3, and is dimensionless; is the coupling coefficient between viscous flow and Knudsen diffusion, with a value ranging from 0.1 to 0.3, and is dimensionless; Let be the coupling coefficient between viscous flow and molecular diffusion, ranging from 0.1 to 0.3, dimensionless, and the sum of the elements of the matrix satisfies . ,in For matrix row index, These are matrix column indices, all dimensionless. The formula for calculating the Knudsen diffusion coefficient is as follows:

[0075] ;

[0076] In the formula, For the first The first grid node Knudsen diffusion coefficient of each gas component, in units of ; For the first The average aperture value of each grid node, in units of The microstructure of the electrode samples was observed using a scanning electron microscope, and statistical analysis was performed using image processing techniques. The range was [missing information]. ~ ; This is the universal gas constant, with a value of 8.314. ; For the first Temperature values ​​of each grid node, in units of ; For the first The molar mass of a gas molecule, in units of The molecular diffusion coefficient is calculated using the Chapman-Ernskorg theory, and the formula is as follows:

[0077] ;

[0078] In the formula, For the first The first grid node species and first Molecular diffusion coefficient between gas components, in units of ; For the first The pressure value of each grid node, in units of ; For the first species and first The equivalent molecular mass of the gas is dimensionless and is calculated using the following method: ; For the first species and first The collision diameter of a gas molecule, in units of By consulting gas molecular dynamics parameter tables, the collision diameter of common gases ranges from 2 to 5. ; For the first species and first The collision integral of the gas molecules is dimensionless and is obtained by consulting the collision integral table and interpolating based on temperature, with a range of 0.5 to 2.5. This is a gas component index, dimensionless, similar to the one mentioned above. Pairs are used to represent different gas components. The calculation of viscous flow permeability uses the Kamenkozeni equation, expressed as follows:

[0079] ;

[0080] In the formula, For the first Viscous flow permeability of each grid node, in units of ; For the first The porosity of each grid node is dimensionless and ranges from 0.3 to 0.5. The effective mass transfer coefficient field of the electrode porous medium includes the Knudsen diffusion coefficient, molecular diffusion coefficient, and viscous flow permeability of each grid node.

[0081] The specific implementation of step S04 involves using the operator splitting method to decompose the strongly coupled multiphysics problem into four weakly coupled subproblems, and solving each subproblem sequentially within each time step. The electrochemical reaction field calculation uses the Butler-Wolmer equation, expressed as follows:

[0082] ;

[0083] In the formula, For the first The grid node at the ... The current density value at each time step, in units of ; This is a time step index, dimensionless; For the first The exchange current density of each grid node, in units of The calculation method is as follows ,in Pre-exponential factors, in units of Experience value ~ , Activation energy, unit: Experience value ~ , For the first Reference concentrations of the gaseous components, in units of , For the first The reaction order of the gaseous components is dimensionless. is the anode transfer coefficient, dimensionless, with an empirical value of 0.5; is the cathode transfer coefficient, dimensionless, with an empirical value of 0.5; The number of electrons transferred in the electrode reaction is dimensionless. and The value of the co-electrolysis reaction is 2; This is the Faraday constant, with a value of 96485. ; For the first Overpotential of each grid node, in units of The calculation method is as follows ,in For the first The electronic phase potential of each grid node, in units of , For the first The ionic phase potential of each grid node, in units of , For the first The equilibrium potential of each grid node, in units of The mass transfer field is calculated using the Nernst equation. The calculation employs Stefan Maxwell's diffusion equation combined with Darcy's law, and the formula is as follows:

[0084] ;

[0085] In the formula, For the first The first grid node The gas components in the first Concentration values ​​at each time step, in units of ; For time variables, the unit is ; For the first The first grid node The effective diffusion coefficient of the gas component, in units of The calculation method is as follows ,in For the first The average molecular diffusion coefficient of this gas component with all other gas components, in units of The calculation method is as follows ; For the first The viscosity value of the gas mixture at each grid node, in units of Calculated using Sutherland's formula; For the first The first grid node The source terms of a gaseous component generated or consumed by electrochemical reactions, in units of The calculation method is as follows ,in For the first The stoichiometric coefficients of the gaseous components are dimensionless and are for and In co-electrolysis, reactants take negative values ​​and products take positive values; For the first The tortuosity of each grid node is dimensionless and calculated as follows: ,in This is a tortuosity correction factor, dimensionless, typically ranging from 0.5 to 1.0. The heat transfer field calculation uses the energy conservation equation, expressed as follows:

[0086] ;

[0087] In the formula, For the first The material density of each grid node, in units of ; For the first The specific heat capacity of the material at each grid node, in units of ; For the first The grid node at the ... Temperature values ​​at each time step, in units of ; For the first Thermal conductivity of the material at each grid node, in units of ; For the first Electrochemical reaction heat power density of each grid node, in units of ; For the first The Joule heat power density of each grid node, in units of The stress field calculation uses the thermoelastic constitutive equation, which is expressed as follows:

[0088] ;

[0089] In the formula, For the first The grid node at the ... The stress value at each time step, in units of ; For the first The elastic modulus of the material at each grid node, in units of ; For the first The mechanical strain of each grid node is dimensionless and calculated from the displacement field gradient. The calculation method is as follows: ,in For the first The displacement vector of each grid node, in units of ; For the first The coefficient of thermal expansion of the material of each grid node, in units of... ; This is a reference temperature value, in units of... The value is usually 298. .

[0090] The specific implementation of step S05 is to calculate the multiphysics coupling anomaly monitoring vector, and the formula is expressed as follows:

[0091] ;

[0092] In the formula, This is a dimensionless multi-physics coupling anomaly monitoring vector. The residual of the electrochemical reaction field is dimensionless and is calculated using the following method: ,in This is a numerically stable term, with units of . The value is ; The mass transfer field residual is dimensionless and is calculated using the following method: ,in This is a numerically stable term, with units of . The value is ; The residual of the heat transfer field is dimensionless and is calculated using the following method: ; The stress field residual is dimensionless and is calculated using the following method: ,in This is a numerically stable term, with units of . The value is The formula for the adaptive time step adjustment function is as follows:

[0093] ;

[0094] ;

[0095] In the formula, The adjusted time step, in units of ; This is the current time step, in units of ; This is a dimensionless adjustment factor for the time step. The adjusted relaxation factor is dimensionless. The current relaxation factor is dimensionless. Let be the relaxation factor adjustment coefficient, which is dimensionless. The adjustment coefficient is calculated as follows: when... hour, and ;when hour, and ;when hour, and ;when hour, and ;in The corresponding preset threshold is dimensionless. This is the maximum residual ratio, which is dimensionless.

[0096] The specific implementation of step S06 involves calculating the heat power density of each grid node based on the updated current density field vector and the updated temperature field vector. The formula for the heat power density of the electrochemical reaction is as follows:

[0097] ;

[0098] In the formula, For the first Electrochemical reaction heat power density of each grid node, in units of ; For the first The grid node at the ... The current density value at each time step, in units of ; For the first Overpotential of each grid node, in units of ; For the first The grid node at the ... Temperature values ​​at each time step, in units of ; The molar entropy change of an electrochemical reaction is expressed in units of 1. ,for Empirical values ​​for electrolysis reactions ,for Empirical values ​​for electrolysis reactions The formula for Joule heat power density is as follows:

[0099] ;

[0100] In the formula, For the first The Joule heat power density of each grid node, in units of ; For the first The electrical conductivity of the material at each grid node, in units of The formula for the cross-scale feedback correction matrix is ​​as follows:

[0101] ;

[0102] In the formula, This is a cross-scale feedback correction matrix, which is dimensionless. is the temperature correction factor for the Knudsen diffusion coefficient, which is dimensionless; This is the temperature-based correction factor for the molecular diffusion coefficient, and it is dimensionless. This is a dimensionless correction factor for the effect of temperature on the permeability of viscous flow. is the correction factor for the Knudsen diffusion coefficient based on the current density, and is dimensionless; This is a dimensionless correction factor for the molecular diffusion coefficient based on the current density. This is a dimensionless correction factor for the current density on the permeability of viscous flow. The formulas for calculating the elements in the first row are as follows:

[0103] ;

[0104] ;

[0105] ;

[0106] In the formula, For the first The initial temperature value of each grid node, in units of ; is the sensitivity coefficient of temperature to the Knudsen diffusion coefficient, dimensionless, with a value of 0.5; is the temperature sensitivity coefficient of the molecular diffusion coefficient, dimensionless, with a value of 1.5; This is the sensitivity coefficient of temperature to the permeability of viscous flow; it is dimensionless and has a value of 0.3. The calculation formula for the elements in the second row is as follows:

[0107] ;

[0108] ;

[0109] ;

[0110] In the formula, is the sensitivity coefficient of current density to the Knudsen diffusion coefficient, which is dimensionless and has a value of 0.2; is the sensitivity coefficient of current density to molecular diffusion coefficient, dimensionless, with a value of 0.4; is the sensitivity coefficient of current density to the permeability of viscous flow, dimensionless, and taken as 0.6. The formula for the corrected effective mass transfer coefficient field is as follows:

[0111] ;

[0112] ;

[0113] ;

[0114] In the formula, For the revised first Knudsen diffusion coefficient for each grid node, in units of ; For the revised first Molecular diffusion coefficient of each grid node, in units of ; For the revised first Viscous flow permeability of each grid node, in units of .

[0115] The specific implementation of step S07 involves calculating the maximum relative changes of the updated temperature field vector, the updated concentration field vector, and the updated current density field vector over three consecutive time steps, as expressed by the following formula:

[0116] ;

[0117] ;

[0118] ;

[0119] In the formula, The maximum relative change in the temperature field over three consecutive time steps is dimensionless. The maximum relative change in the concentration field over three consecutive time steps is dimensionless. The current density field is the maximum relative change over three consecutive time steps, and is dimensionless. This is a dimensionless time step index variable used to represent three consecutive time steps. The convergence criterion is... and and ,in The convergence criterion for the temperature field is dimensionless and set to 0.01. The concentration field convergence criterion is dimensionless and set to 0.02. The current density field convergence criterion is dimensionless and set to 0.015.

[0120] The principle behind the formula system proposed in this invention lies in achieving accurate numerical simulation of complex physicochemical processes in a solid oxide battery co-electrolysis system through multi-scale coupling and adaptive adjustment strategies. First, normalization eliminates dimensional differences between different physical quantities and improves numerical stability. Then, a multi-mechanism mass transfer coupling matrix accurately describes the synergistic effect of three mass transfer mechanisms—Knudsen diffusion, molecular diffusion, and viscous flow—within the porous electrode. The Knudsen diffusion coefficient formula is detailed in this invention. Based on gas molecular dynamics theory, considering the effects of pore size, temperature, and molecular mass on diffusion, Knudsen diffusion dominates when the pore size is smaller than the mean free path of gas molecules. The molecular diffusion coefficient formula is... The effect of intermolecular collisions on diffusion is described using Chapman-Ninscoger theory; the permeability formula for viscous flow is given. The Kamankozeny equations are used to describe the flow characteristics of porous media. Then, the operator splitting method is employed to decouple the strongly coupled multiphysics problem into multiple independently solvable subproblems, thereby reducing computational complexity. Among these subproblems, the Butler-Wolmer equations are used. The Stefan-Maxwell diffusion equation describes the electrochemical reaction kinetics at the electrode surface. The diffusion and convection mass transfer processes of gas within the porous electrode are described, along with the energy conservation equation. The thermoelastic constitutive equation describes the heat generation and transfer process inside the battery. The distribution of thermal stress caused by temperature changes is described, and the dynamic influence of macroscopic temperature field and current density field on microscopic mass transfer parameters is realized through a cross-scale feedback correction matrix. The correction formula is... , and This invention demonstrates the feedback effect of macroscopic field variables on the mass transfer characteristics of microscopic porous structures. Finally, an adaptive time step adjustment function dynamically optimizes the solution parameters based on the computational convergence to ensure numerical stability and computational efficiency. The overall effect of the formula system is that it can accurately capture the strong coupling relationship between electrochemical reactions, mass transfer, heat transfer, and stress fields during the co-electrolysis of solid oxide batteries, predict the temperature distribution, gas concentration distribution, current density distribution, and stress distribution of the system, and provide a reliable numerical simulation tool for battery structure optimization and operating condition optimization. Compared with traditional single-scale or single-physics field simulation methods, the formula system of this invention accurately describes the synergistic effect of different mass transfer mechanisms by introducing a multi-mechanism mass transfer coupling matrix, realizes the dynamic correction of microscopic parameters by the macroscopic field through a cross-scale feedback correction matrix, and significantly improves computational efficiency and numerical stability through an adaptive adjustment strategy. This keeps the relative error between the numerical simulation results and experimental data within 5%, and reduces the calculation time by more than 30%. At the same time, the multi-physics field coupling anomaly monitoring vector monitors numerical anomalies in real time during the calculation process and adjusts the solution parameters in a timely manner, avoiding computational divergence and pseudo-convergence, and ensuring the reliability and accuracy of the simulation results.

[0121] To better understand and implement this invention, the following is a specific application scenario of this invention, Example 2:

[0122] A technical team conducted a numerical simulation study on the mass and heat transfer coupling of a solid oxide battery co-electrolysis system operating at 850℃. This system is used for... and Syngas was produced by co-electrolysis. The battery adopted an anode-supported structure. The electrolyte layer was made of yttrium-stabilized zirconium oxide with a thickness of 8 μm. The anode functional layer was made of nickel-zirconia composite material with a thickness of 10 μm and a porosity of 35%. The anode diffusion layer was made of nickel-zirconia support with a thickness of 450 μm and a porosity of 45%. The cathode functional layer was made of lanthanum-strontium-cobalt ferrite material with a thickness of 12 μm and a porosity of 32%. The cathode diffusion layer was made of lanthanum-strontium-cobalt ferrite support with a thickness of 400 μm and a porosity of 42%. The technical team established a three-dimensional geometric model of the battery with a computational domain size of 25 × 25 × 0.88 mm. A uniform mesh size of 20 μm was used for the electrolyte layer. A denser mesh size of 8 μm was used in the interface region between the anode functional layer and the electrolyte. A denser mesh size of 7 μm was used in the interface region between the cathode functional layer and the electrolyte. A sparse mesh size of 35 μm was used in the anode diffusion layer far from the interface region. A sparse mesh size of 32 μm was used in the cathode diffusion layer far from the interface region. The total number of meshes reached 1.2 × 10⁻⁶ mm. indivual.

[0123] The system is initially set to an anode intake configuration. Volume fraction 50%, Volume fraction 30%, Volume fraction 15%, Volume fraction 5%, flow rate 1.5× The cathode intake airflow is 8.0 × Operating voltage 1.35V, average current density 4500 The technical team collected temperature distribution data under initial operating conditions, ranging from 823℃ to 878℃. The lowest system temperature of 823℃ was located at the inlet of the anode diffusion layer, and the highest system temperature of 878℃ ​​was located in the central region of the cathode functional layer. The temperature values ​​at each grid node were normalized by subtracting 823℃ and then dividing by 55℃, resulting in an initial temperature field vector containing 1.2 × Each element. Gas component concentration distribution data at the anode outlet. Volume fraction increased to 42%. Volume fraction increased to 18%. Volume fraction reduced to 28%, The volume fraction was reduced to 12%, and the gas component concentration values ​​at each grid node were normalized by dividing the inlet concentration value to form the initial concentration field vector. Current density distribution data showed that the local current density reached a peak of 5800 at the edge region of the cathode functional layer. The value drops to a valley of 3600 in the central region of the anode diffusion layer. The current density value of each grid node divided by the average current density 4500 After normalization, an initial current density field vector is formed.

[0124] Based on the initial temperature and concentration field vectors, the technical team calculated the local pore size distribution parameters within the anode functional layer. The average pore size was 0.52 μm with a standard deviation of 0.18 μm. Within the anode diffusion layer, the average pore size was 2.1 μm with a standard deviation of 0.65 μm. Within the cathode functional layer, the average pore size was 0.48 μm with a standard deviation of 0.16 μm. Within the cathode diffusion layer, the average pore size was 1.9 μm with a standard deviation of 0.58 μm. The gas state parameters for each grid node included temperature, pressure (101325 Pa), density, and viscosity. Figure 2 As shown, the technical team calculated the Knudsen diffusion coefficient of the anode functional layer using a multi-mechanism mass transfer coupling matrix, finding it to be in the range of 2.1 × 10⁻⁶. Up to 3.8× The molecular diffusion coefficient ranges from 1.2× Up to 1.8× The viscous flow permeability range is 3.5× Up to 6.2× The diagonal elements of the multi-mechanism mass transfer coupling matrix have values ​​of 0.75, 0.80, and 0.70, respectively, while the off-diagonal elements have values ​​between 0.15 and 0.25.

[0125] The technical team employed the operator splitting method for multiphysics coupling solutions, setting the initial time step to 0.05 s and the initial relaxation factor to 0.6. In the electrochemical reaction field calculations, the Butler-Wolmer equations were used to obtain the anode exchange current density of 850. The cathode exchange current density is 1200 The overpotential distribution of the anodic reaction ranges from 0.08 to 0.15 V, while the overpotential distribution of the cathode reaction ranges from 0.05 to 0.11 V. In the mass transfer field calculation, the Stefan-Maxwell diffusion equation combined with Darcy's law is used to describe the gas transport process, and the calculated values ​​within the anodic functional layer are obtained. The maximum concentration gradient is 2.3 × , The maximum concentration gradient is 1.8 × In the heat transfer field calculation, the energy conservation equation considers the heat of electrochemical reaction, Joule heating, convective heat transfer, and radiative heat transfer. The calculated power density of the electrochemical reaction heat in the anode functional layer ranges from -1.2 × 10⁻⁶. to -3.5× Negative values ​​indicate endothermic reactions, and the range of heat power density for electrochemical reactions in the cathode functional layer is 8.5 × Up to 2.1× A positive value indicates an exothermic reaction, and the Joule heat power density of the electrolyte layer reaches 4.2 × 10⁻⁶. In the stress field calculation, the maximum principal stress at the interface between the electrolyte layer and the anode functional layer was calculated to be 185 MPa, and the maximum principal stress at the interface between the electrolyte layer and the cathode functional layer was 172 MPa, using the thermoelastic constitutive equation.

[0126] After the first iteration, the updated temperature field vector shows that the system's highest temperature has risen to 882℃, and the updated concentration field vector shows that the anode outlet... The volume fraction increased to 43%, and the updated current density field vector showed that the average current density increased to 4580. The updated stress field vector showed that the maximum principal stress of the electrolyte layer increased to 192 MPa. The technical team calculated the multi-physics coupling anomaly monitoring vector. The residuals for the electrochemical reaction field were 0.068, the mass transfer field residual was 0.092, the heat transfer field residual was 0.071, and the stress field residual was 0.115. Among them, the mass transfer field residual and the stress field residual exceeded the corresponding preset thresholds. Dividing the mass transfer field residual by the preset threshold of 0.08 yielded a ratio of 1.15, and dividing the stress field residual by the preset threshold of 0.10 yielded a ratio of 1.15. Both ratios were within the range of 1 to 1.5. According to the adaptive time step adjustment function logic, the adjusted time step was 0.05 × 0.8, or 0.04 s, and the adjusted relaxation factor was 0.6 × 1.1, or 0.66.

[0127] like Figure 3 As shown, the technical team calculated the electrochemical reaction heat power density and Joule heat power density of each grid node based on the updated current density field vector and the updated temperature field vector, and assigned the updated temperature field vector to the initial temperature field vector. The team then used a cross-scale feedback correction matrix to correct the effective mass transfer coefficient field of the porous electrode medium. In the first row of this matrix, the temperature sensitivity coefficients were set to 0.5 for the Knudsen diffusion coefficient, 1.5 for the molecular diffusion coefficient, and 0.3 for the viscous flow permeability. In the second row, the current density sensitivity coefficients were set to 0.2 for the Knudsen diffusion coefficient, 0.4 for the molecular diffusion coefficient, and 0.6 for the viscous flow permeability. The corrected Knudsen diffusion coefficient range for the anode functional layer was updated to 2.3 × 10⁻⁶. Up to 4.1× The molecular diffusion coefficient range has been updated to 1.4× Up to 2.1× The permeability range for viscous flow has been updated to 3.8× Up to 6.8× .

[0128] The technical team continued iterative solving, and under the conditions of an adjusted time step of 0.04 s and an adjusted relaxation factor of 0.66, the residuals of each physical field gradually decreased. After 128 iterations, the maximum relative change of the updated temperature field vector within three consecutive time steps decreased to 0.008, the maximum relative change of the updated concentration field vector decreased to 0.016, and the maximum relative change of the updated current density field vector decreased to 0.012. All three indicators were less than their respective convergence criteria, and the system reached steady state. The steady-state numerical simulation results show that the system temperature distribution ranges from 825℃ to 885℃, and the composition of the anode outlet gas is... Volume fraction 44%, Volume fraction 19%, Volume fraction 26%, Volume fraction 11%, electrolysis voltage 1.35V, electrolysis current 28.1A. Gas production rate 1.8 × , Gas production rate 2.1× The energy conversion efficiency is 82%, and the maximum principal stress of the electrolyte layer, 198 MPa, is located in the central region of the interface with the anode functional layer. The technical team plotted performance curves under different operating parameters based on steady-state numerical simulation results, as shown in Table 1.

[0129] Table 1 Steady-state performance parameters under different operating voltages

[0130]

[0131] The advancements of this invention compared to traditional numerical simulation methods for co-electrolysis of solid oxide batteries are reflected in several aspects. Traditional methods typically simplify the mass transfer process to a single effective diffusion coefficient, neglecting the coupling effects of Knudsen diffusion, molecular diffusion, and viscous flow—three mass transfer mechanisms—leading to insufficient prediction accuracy when electrode pore size distribution is uneven. This invention, by introducing a multi-mechanism mass transfer coupling matrix, can dynamically calculate the contribution of each mass transfer mechanism based on local pore size distribution parameters and gas state parameters, more realistically reflecting the complex gas transport process within the porous electrode medium. Traditional methods often use fixed time steps and relaxation factors when solving strongly coupled multiphysics problems, which can easily lead to overall solution failure when a subproblem exhibits numerical oscillations or convergence difficulties. This invention, by setting a multiphysics coupling anomaly monitoring vector to monitor the residual changes of each subproblem in real time, and adaptively adjusts the time step and relaxation factor based on the ratio of the residual to a preset threshold, improving solution efficiency while ensuring computational convergence. Traditional methods treat the mass transfer parameters of the microporous structure as constants, neglecting the influence of macroscopic temperature and current density field changes on microscopic mass transfer characteristics. This invention establishes a dynamic correlation between macroscopic field variables and microscopic mass transfer parameters through a cross-scale feedback correction matrix. When the temperature rises or the current density distribution changes during battery operation, it can promptly correct the Knudsen diffusion coefficient, molecular diffusion coefficient, and viscous flow permeability, enabling numerical simulation results to accurately capture the evolution of battery performance over time. Traditional methods typically only focus on the electrochemical reaction field and heat transfer field, neglecting the stress field and failing to predict the risk of material failure due to thermal stress accumulation during long-term battery operation. This invention incorporates stress field calculation into the operator splitting method's solution framework and uses the stress field residual as a component of the multiphysics coupling anomaly monitoring vector. This allows for the simultaneous acquisition of the battery's internal temperature distribution, concentration distribution, current density distribution, and stress distribution, providing a more comprehensive basis for battery structure optimization design and reliability assessment.

[0132] 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 changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A numerical simulation method for co-electrolysis of solid oxide batteries considering mass and heat transfer coupling, characterized in that, A three-dimensional geometric model of the solid oxide battery co-electrolysis system was obtained, and the electrolyte layer, anode functional layer, anode diffusion layer, cathode functional layer, and cathode diffusion layer were meshed. Temperature distribution data, gas component concentration distribution data, and current density distribution data under initial operating conditions were collected, normalized, and stored as initial temperature field vector, initial concentration field vector, and initial current density field vector, respectively. Based on the initial temperature field vector and initial concentration field vector, local pore size distribution parameters and gas state parameters within the porous electrode were calculated. The Knudsen diffusion coefficient, molecular diffusion coefficient, and viscous flow permeability were calculated using a multi-mechanism mass transfer coupling matrix to obtain the effective mass transfer coefficient field of the porous medium of the electrode. This was then processed using operators... The splitting method is used to calculate the electrochemical reaction field, mass transfer field, heat transfer field, and stress field to obtain updated temperature field vector, updated concentration field vector, updated current density field vector, and updated stress field vector. The multiphysics coupling anomaly monitoring vector is calculated, and the adjusted time step and adjusted relaxation factor are obtained through an adaptive time step adjustment function. The cross-scale feedback correction matrix is ​​called to correct the effective mass transfer coefficient field of the electrode porous medium to obtain the corrected effective mass transfer coefficient field. When the maximum relative change of the updated temperature field vector, updated concentration field vector, and updated current density field vector within three consecutive time steps is less than the convergence criterion, the steady-state numerical simulation result is output.

2. The method according to claim 1, characterized in that, The porous medium region of the electrode is subjected to non-uniform mesh refinement treatment. In the interface region between the electrode and the electrolyte and in the region with a large porosity gradient inside the electrode, the mesh size is set to 0.3 to 0.5 times the mesh size of the conventional region.

3. The method according to claim 2, characterized in that, In the region of the electrode diffusion layer far from the interface, the mesh size is set to 1.5 to 2 times the mesh size of the regular region, where the mesh size of the regular region is ∈ [10 μm, 50 μm].

4. The method according to claim 3, characterized in that, The normalization process for the initial temperature field vector involves subtracting the system's lowest temperature value from the temperature value of each grid node in the temperature distribution data, and then dividing by the temperature range value, which is the difference between the system's highest and lowest temperatures.

5. The method according to claim 4, characterized in that, The normalization process for the initial concentration field vector is to divide the concentration value of each gas component at each grid node in the gas component concentration distribution data by the concentration value of the gas component at the inlet.

6. The method according to claim 5, characterized in that, The normalization process for the initial current density field vector is to divide the current density value of each grid node in the current density distribution data by the average current density value.

7. The method according to claim 6, characterized in that, The multi-mechanism mass transfer coupling matrix is ​​a 3×3 matrix, with diagonal elements ranging from [0.6, 0.9] and off-diagonal elements ranging from [0.1, 0.3]. The sum of all elements in the matrix is ​​normalized to 3.

8. The method according to claim 7, characterized in that, The Knudsen diffusion coefficient is calculated based on the local pore size distribution parameters and the mass of gas molecules. The molecular diffusion coefficient is calculated based on the gas temperature and gas pressure values, determined by the Chapman-Ninscog theory. The viscous flow permeability is calculated based on the porosity and tortuosity of the porous medium, determined by the Karman-Kozani equation.

9. The method according to claim 8, characterized in that, Electrochemical reaction field calculations use the Butler-Wolmer equation to describe electrode reaction kinetics, while mass transfer field calculations use the Stefan-Maxwell diffusion equation combined with Darcy's law to describe gas transport processes.

10. The method according to claim 9, characterized in that, The heat transfer field calculation uses the energy conservation equation and considers electrochemical reaction heat, Joule heat, convective heat transfer, and radiative heat transfer. The stress field calculation uses the thermoelastic constitutive equation to describe the mechanical behavior of the material.

Images