Coupled water vapor-liquid phase change heat, oxygen transport and dynamic reaction fuel cell cathode pore scale simulation method, medium and system

Abstract

The present application belongs to the technical field of fuel cell, and particularly relates to a fuel cell cathode pore scale simulation method, medium and system coupling water gas-liquid phase change heat, oxygen transmission and dynamic reaction, comprising: step 1, constructing a calculation domain; step 2, setting initial conditions and boundary conditions; step 3, updating the discrete velocity distribution function of each lattice point; step 4, updating the macroscopic density, macroscopic velocity and macroscopic pressure of water of each lattice point; step 5, updating the viscosity, specific heat capacity and thermal conductivity of all lattice points; step 6, updating the macroscopic temperature of each lattice point; step 7, updating the discrete concentration distribution function of each lattice point; step 8, updating the macroscopic concentration of oxygen of each lattice point; step 9, updating the dynamic reaction rate, and determining the termination of iteration. Compared with the prior art, the present application solves the problem that the prior art is difficult to characterize the multi-component reaction transmission in the real structure of the fuel cell cathode. The present application solves the interaction and micro mechanism.

Description

Technical Field

[0001] This invention belongs to the field of fuel cell technology, specifically relating to a method, medium, and system for simulating the cathode pore size of a fuel cell that couples water-gas-liquid phase change heat transfer, oxygen transport, and dynamic reactions. Background Technology

[0002] Fuel cells, due to their cleanliness and high efficiency, are widely considered a promising new energy device for future transportation applications. However, achieving large-scale commercial application of fuel cells still faces many challenges, one of the core issues being further improving their power density. Key factors limiting fuel cell performance improvement mainly focus on water and thermal management at the cathode. The presence of water is crucial for maintaining good hydration of the proton exchange membrane, but excessive liquid water significantly hinders oxygen transport within the porous electrode, leading to severe water accumulation and reduced battery reaction performance. Simultaneously, the battery operating temperature needs to be controlled within a reasonable range to ensure its stability and efficiency. Since water is both an important material carrier and a key heat carrier in the heat and mass transport process, there is a significant coupling effect between its gas-liquid phase changes, latent heat of phase change, and heat transport behavior. Therefore, studying the microscopic coupling mechanism between water-gas-liquid phase change transport, oxygen transport, heat transport, and dynamic electrochemical reactions within the fuel cell cathode is of great significance for regulating the fuel cell's hydrothermal balance, improving reaction transport processes, and enhancing battery performance.

[0003] However, existing pore-scale simulation methods still have significant limitations in solving the complex multiphysics coupling problem of fuel cell cathodes. On one hand, fuel cells possess complex pore geometry, component distribution, and multiphase interface characteristics, making it difficult for existing methods to simultaneously and accurately characterize the conjugate heat and mass transfer at complex fluid-solid interfaces, as well as the hindering effect of water-gas-liquid phase transitions on oxygen transport. Simultaneously, the dynamic evolution of the water-gas-liquid phase transition state is closely related to local temperature and pressure changes, and its latent heat of phase transition has a significant impact on the heat transfer process; however, existing pore-scale simulation methods often neglect these effects. Furthermore, there are significant coupling relationships between water flow, oxygen diffusion, heat transfer, and dynamic reactions within complex pore structures, and existing methods still fall short in solving multiphysics coupling problems, thus making it difficult to systematically reveal the influence mechanism of microstructural changes on reaction transport characteristics and actual performance. For example, CN119476118A discloses a mesoscale simulation method for the diffusion layer of a fuel cell, establishing a mesoscale simulation framework that considers the coupled transport of water, oxygen, and heat within the diffusion layer in relation to electrochemical reactions. However, this scheme does not consider the influence of the gas-liquid phase change and latent heat effect of water on the reaction transport behavior, and ignores the influence of water distribution on oxygen transport. This makes the results obtained by this method prone to deviating from the local heat accumulation and oxygen transport characteristics of the actual physical process, thereby reducing the characterization accuracy of the actual multi-physics coupled transport process of the fuel cell. Similarly, CN115906700A discloses a simulation method, system, device, and medium for the gas diffusion layer. Although this scheme can obtain relevant data on water vapor within the fuel cell, it does not further consider the conjugate heat transfer process between fluid-solid surfaces in complex porous structures, nor does it introduce the coupling effect of water-liquid phase change behavior on the heat and mass transport process. It focuses more on the statistical analysis of changes in macroscopic parameters of the fuel cell, thus making it difficult to effectively reveal the microscopic coupling mechanism of the multi-component reaction transport process inside the fuel cell. For example, CN121809257A discloses a mesoscale simulation method, medium, and system for dynamic coupling of heat, mass, and reaction within a fuel cell catalyst layer. This method considers the oxygen-heat coupling transport process at the nanoscale within the catalyst layer, but neglects the influence of water generation on oxygen and heat transport. In particular, in actual processes, the liquid form of water can block the effective oxygen transport path, and the gaseous water transport process can also affect the diffusion capacity of oxygen within the pores. Therefore, this method has significant biases in predicting the multi-component reaction transport process within actual fuel cells. Because pore-scale studies need to consider water transport and gas-liquid phase change processes, fluid-solid conjugate heat transfer, complex phase interface evolution and transport differences, and dynamic reaction processes within complex porous structures, these factors interact and jointly affect the overall reaction performance of the fuel cell. The coupling between different processes is difficult, therefore, achieving effective simulation of such problems remains a significant challenge.

[0004] Currently, the following technical challenges remain in the multiphysics pore-scale simulation of coupled water-gas-liquid phase change heat transfer, oxygen transport, and dynamic electrochemical reactions within the fuel cell cathode: (1) The fuel cell cathode has complex pore structure characteristics, component distribution characteristics, and multiphase interface transport differences. Existing models often struggle to simultaneously account for the conjugate heat transfer and reaction transport processes at complex fluid-solid interfaces, as well as the transport barrier effect between water and oxygen, resulting in inaccurate descriptions of the reaction transport processes at complex multi-component boundaries. (2) The gas-liquid phase change behavior of water within the fuel cell cathode is closely related to local temperature and pressure changes, while the latent heat of vaporization has a significant impact on the heat transfer process. Existing simulation methods are still insufficient in characterizing the dynamic evolution of the water-gas-liquid state and the latent heat effect, making it difficult to accurately reveal the true water-thermal coupling effect and its microscopic mechanism within the fuel cell. (3) There are significant coupling relationships and feedback mechanisms among water flow, oxygen diffusion, heat transfer, and dynamic reactions within the complex pore structure of the fuel cell cathode. Existing models still have shortcomings in multiphysics coupling solutions and information interaction, making it difficult to systematically reveal the influence of changes in the microstructure of the fuel cell cathode on the multiphysics reaction transport characteristics and actual performance.

[0005] Therefore, there is an urgent need to develop a pore-scale simulation method that can characterize the true internal structural features of the fuel cell cathode and couple water-gas-liquid phase change heat transfer, oxygen transport, and dynamic reactions, so as to achieve an accurate description of the multi-physics coupling process and its microscopic mechanism within the fuel cell cathode. Summary of the Invention

[0006] The purpose of this invention is to address at least one of the aforementioned problems by providing a method, medium, and system for simulating the pore-scale of a fuel cell cathode that couples water-gas-liquid phase change heat transfer, oxygen transport, and dynamic reactions. This addresses the shortcomings of existing technologies in characterizing the multi-component reaction transport within the actual structure of the fuel cell cathode, while simultaneously considering the conjugate heat transfer and reaction transport processes at complex fluid-solid interfaces, the transport barrier effect between water and oxygen, and the insufficient characterization of the dynamic evolution and latent heat effects of the water-gas-liquid state. This solution can effectively solve for the interactions and microscopic mechanisms among water-gas-liquid phase change heat transfer, oxygen transport, and dynamic reactions within the complex pore structure of the fuel cell cathode.

[0007] The objective of this invention is achieved through the following technical solution: The first aspect of this invention discloses a method for simulating the cathode pore size of a fuel cell that couples water-gas-liquid phase change heat transfer, oxygen transport, and dynamic reactions, comprising the following steps: Step 1: Construct the diffusion layer structure using a random algorithm and generate a gas flow layer, microporous layer, and catalyst layer that match the diffusion layer to form the overall structure of the fuel cell cathode; binarize the overall structure of the fuel cell cathode and map it proportionally to construct the computational domain; Step 2: Set initial and boundary conditions based on the computational domain constructed in Step 1. The initial conditions include reaction parameters, physical property parameters, macroscopic parameters, and the calculated dynamic reaction rate. The boundary conditions are as follows: the top of the computational domain is set as a constant oxygen concentration boundary, a constant temperature boundary, and a no-slip bounce boundary condition for water flow; the bottom of the computational domain is set as a dynamic reaction boundary; and the left and right boundaries of the computational domain are set as periodic boundaries. The flow-solid boundary within the computational domain is set as a no-slip bounce boundary during the calculation of water flow and oxygen diffusion, and as a conjugate heat transfer boundary during the calculation of heat transfer. Step 3: Based on the initial conditions set in Step 2, construct the discrete velocity distribution function and equilibrium velocity distribution function of water under the target working condition. Then, use the evolution equation to perform collision and migration calculations on all grid points in the computational domain to update the discrete velocity distribution function of each grid point. Step 4: Based on the results of Step 3, update the macroscopic density, macroscopic velocity, and macroscopic pressure of water at each grid point; Step 5: Based on the results of Step 4, update the viscosity, specific heat capacity, and thermal conductivity of all grid points; Step 6: Substitute the results of Step 2, Step 4 and Step 5 into the temperature evolution equation to perform collision and migration calculations on all grid points in the computational domain and update the macroscopic temperature of each grid point. Step 7: Determine the new effective diffusion region based on the macroscopic density of water calculated in Step 4; construct the discrete concentration distribution function and equilibrium concentration distribution function of oxygen under the target operating condition according to the initial conditions set in Step 2; then use the evolution equation to perform collision and migration calculations on all grid points to update the discrete concentration distribution function of each grid point. Step 8: Based on the results of Step 7, calculate and update the macroscopic oxygen concentration at each grid point; Step 9: Substitute the results of Step 6 and Step 8 into Step 2 to update the dynamic reaction rate, and repeat Step 2 to Step 8 for iteration; when the change in macroscopic velocity obtained by iteration does not exceed the set threshold, it is determined that a stable state has been reached, and the iteration is terminated.

[0008] Preferably, in step 2: The dynamic reaction rate includes: water generation rate Oxygen consumption rate and heat generation rate ; Wherein: oxygen consumption rate Based on oxygen reaction flux The calculated water generation rate Based on water reaction flux The calculated heat generation rate Based on oxygen consumption rate Calculated.

[0009] Preferably, in step 3: Determine the initial discrete velocity distribution function : ; Equilibrium velocity distribution function for: ; In the formula: Indicates the first Each discrete velocity direction corresponds to a weight factor, and the sum of all weight factors is 1; This indicates the macroscopic density of water; Indicates macroscopic speed; The spatial vector distribution representing discrete velocities, and The modulus is 1; , This indicates the lattice velocity, which is converted based on the actual speed of sound. The evolution equation is: ; In the formula: This is the collision relaxation term, used to describe the relaxation process of the distribution function from a non-equilibrium state to an equilibrium state under the action of local collisions; Indicates the time step; Represents the discrete force terms for different components; Represents discrete reaction source terms; Indicates discrete additional terms; ; It is the identity matrix. It is a diagonal matrix (which includes the relaxation factor of the model).

[0010] Throughout the computational model, the initial discrete velocity distribution function From the equilibrium velocity distribution function The calculations show that, as time evolves, new equilibrium velocity distribution functions are continuously calculated. Substitute these values ​​into the evolution equation to update the values ​​corresponding to different lattice points. And based on the updated Further calculation of macroscopic density and macro speed Then, the new equilibrium velocity distribution function is calculated using the updated macroscopic velocity and macroscopic density. This iterative calculation continues. Therefore, according to this logic, the equilibrium distribution function is only passed initially. Calculate the corresponding initial discrete velocity distribution function. The subsequent discrete velocity distribution function The calculation is then continuously updated and iterated based on the evolution equation.

[0011] Preferably, in step 4: Macroscopic density of water for: ; Macro speed for: ; In the formula: The discrete velocity distribution function of water; This indicates the macroscopic density of water; The spatial vector distribution representing discrete velocities, and The modulus is 1; Indicates the time step; Represents a force vector; When the lattice point is a fluid region, the macroscopic pressure is... According to macro density and temperature Calculated using state equations; When the lattice point is a solid region, the macroscopic pressure is... .

[0012] Preferably, in step 5: When the lattice points are fluid regions, then: Viscosity for: ; Specific heat capacity for: ; thermal conductivity for: ; In the formula: This represents the density of vapor water under saturated conditions. This represents the density of liquid water in a saturated state. This represents the viscosity of vapor water under saturated conditions. This indicates the viscosity of liquid water under saturated conditions. This represents the specific heat capacity of vapor water under saturated conditions. This is expressed as the specific heat capacity of liquid water under saturated conditions. This represents the thermal conductivity of vapor water under saturated conditions. This represents the thermal conductivity of liquid water under saturated conditions. When the lattice points are solid regions, the density is... Specific heat capacity and thermal conductivity It is a fixed value.

[0013] Preferably, in step 6: The temperature evolution equation is: ; In the formula: Represents the temperature evolution operator; When the lattice point is a fluid region, the temperature evolution operator simultaneously considers the convection term, the heat conduction term, the volume work term controlling the latent heat of vaporization, and the reaction source term. When the lattice point is a solid region, the temperature evolution operator only considers the thermal conductivity term; The macroscopic temperature of each grid point is calculated and updated using the fourth-order Runge-Kutta method.

[0014] Preferably, in step 7: The equilibrium concentration distribution function is: ; In the formula: Indicates the first Each discrete concentration direction has a weighting factor, and the sum of all weighting factors is 1; Indicates the concentration of oxygen; The evolution equation is: ; In the formula: Indicates the time step; Indicates the rate of oxygen consumption; Indicates the relaxation time.

[0015] Preferably, in step 8: The macroscopic concentration of oxygen is: ; In the formula: This represents the discrete concentration distribution function of oxygen.

[0016] A second aspect of the present invention discloses a computer-readable storage medium storing one or more programs that can be executed by one or more processors to implement the steps of the fuel cell cathode pore size simulation method as described in any of the preceding claims.

[0017] A third aspect of the present invention discloses a fuel cell cathode pore size simulation system, including a processor and a memory for executable instructions, wherein the processor executes the instructions to implement the steps of the fuel cell cathode pore size simulation method as described in any of the preceding claims.

[0018] Compared with the prior art, the present invention has the following beneficial effects: 1. This method employs a stochastic algorithm to construct a diffusion layer with realistic structural characteristics, and combines it with a simplified gas flow layer, microporous layer, and catalyst layer to effectively form the overall structure of the fuel cell cathode. Subsequently, by combining characterization methods for fluid-solid conjugate heat transfer boundaries, no-slip rebound transport boundaries, and the effect of water distribution on oxygen transport, the reaction transport behavior of complex boundaries between different components in the fuel cell cathode can be effectively identified, thereby significantly improving the simulation accuracy and applicability of this simulation method for complex boundary problems of multiple components inside the fuel cell cathode.

[0019] 2. This method introduces a state equation to correlate the gas-liquid phase transition process of water with local temperature and pressure changes, realizing the dynamic evolution of the water gas-liquid state under different operating conditions. Furthermore, it considers the influence of the latent heat of water vaporization on the heat transfer process, thereby revealing the water-thermal coupling effect and its microscopic mechanism inside the fuel cell cathode more realistically.

[0020] 3. This method effectively couples the water flow, oxygen diffusion, heat transfer, and dynamic reaction processes within the complex porous structure of the fuel cell cathode, achieving information interaction and feedback between multiple physical fields at the pore scale. By adjusting different parameters, the influence of changes in the fuel cell's microstructure on the multi-physics reaction transport characteristics and its mechanism of action on actual reaction performance can be intuitively characterized, thus providing important support for optimizing the design of the fuel cell cathode structure to improve reaction performance. Attached Figure Description

[0021] Figure 1 This is a flowchart illustrating the method for simulating the cathode pore size of a fuel cell.

[0022] Figure 2 This is an image of the cathode structure of the fuel cell used in Example 2.

[0023] Figure 3 This is an image of the water distribution in the steady state of the fuel cell cathode structure in Example 2.

[0024] Figure 4 This is a temperature distribution image of the fuel cell cathode structure in steady state in Example 2.

[0025] Figure 5 This is an image of the oxygen concentration distribution in steady state of the fuel cell cathode structure in Example 2.

[0026] Figure 6 The curves showing the oxygen consumption rate variation of the fuel cell cathode structure in Example 2 at different time steps are shown. Detailed Implementation

[0027] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments. Unless otherwise specified in the following description, the methods used are known in the art, and any matters not covered herein can be addressed using existing technologies.

[0028] Example 1 This invention proposes a method for simulating the pore-scale structure of a fuel cell cathode by coupling water-gas-liquid phase change heat transfer, oxygen transport, and dynamic electrochemical reactions. This method utilizes a stochastic algorithm to generate structural features consistent with the actual diffusion layer distribution, and combines simplified gas flow layer, micropore layer, and catalyst layer scale features to effectively characterize the overall structural properties of the fuel cell. This method can effectively solve for the interactions and microscopic mechanisms between water-gas-liquid phase change heat transfer, oxygen transport, and dynamic reactions within the complex pore structure of the fuel cell cathode.

[0029] Step 1: A random algorithm is used to construct the carbon fiber distribution characteristics in the diffusion layer, and the parameter ranges of its thickness, porosity, and fiber diameter are controlled. This is combined with a simplified gas flow layer (e.g., a flow channel region of a certain height is set in the top region of the diffusion layer, and a constant oxygen concentration injection condition is set at the top boundary to represent the simplified gas flow layer), a microporous layer (e.g., a micrometer-scale crack structure is considered in the bottom region of the diffusion layer to represent the simplified microporous layer), and a catalyst layer (e.g., an electrochemical reaction boundary is set in the bottom region of the microporous layer to represent the simplified catalyst layer structure) to form the overall structure of the fuel cell cathode. The above simplification methods are all conventional in this field. Then, the generated structure is binarized into a binary matrix corresponding to the pixel points, and proportionally mapped to the lattice points in the computational domain, serving as the basic computational domain for this simulation method.

[0030] Step 2: Based on the basic computational domain determined in Step 1, set the initial and boundary conditions for the simulation process. Specifically, this includes: reaction parameters (reference exchange current density, overpotential, transfer coefficient, reference oxygen concentration, enthalpy of reaction, gas constant, Faraday constant, activation energy), physical property parameters (viscosity, thermal conductivity, specific heat capacity, eccentricity factor, diffusion coefficient, fluid-solid wetting angle), and macroscopic parameters (density of water, molar mass of water, velocity, oxygen concentration, molar mass of oxygen, temperature, pressure), thereby calculating the dynamic reaction rate (water generation rate, oxygen consumption rate, and heat generation rate) under the target operating condition.

[0031] Step 3: Based on the reaction parameters, physical property parameters, macroscopic parameters, and water generation rate determined in Step 2, determine the discrete velocity distribution function, equilibrium velocity distribution function, phase transition related force term, fluid-solid related force term, reaction source term, and additional term for water under the target operating condition. Then, use the evolution equation of the discrete velocity distribution function to perform collision and migration calculations on all lattice points in the computational domain to obtain the new discrete velocity distribution function for each lattice point.

[0032] Step 4: Calculate the macroscopic density, macroscopic velocity, and macroscopic pressure of the water at this point using the discrete velocity distribution function of all grid points calculated in Step 3.

[0033] Step 5: Update the new viscosity, specific heat capacity, and thermal conductivity of all lattice points using the macroscopic density of water calculated in Step 4.

[0034] Step 6: Substitute the heat generation rate calculated in Step 2 into the reaction source term in the temperature evolution equation. Substitute the macroscopic density, macroscopic velocity, and macroscopic pressure of water calculated in Step 4, and the specific heat capacity and thermal conductivity calculated in Step 5 into the heat conduction term, convection term, and volume work term related to latent heat of vaporization in the temperature evolution equation, respectively. Then, traverse all lattice points through the collision migration process to obtain the new macroscopic temperature of each lattice point.

[0035] Step 7: Based on the macroscopic density of water calculated in Step 4, determine the fluid regions not covered by water in all grid points as new effective diffusion regions. Based on the reaction parameters, physical property parameters, macroscopic parameters, and oxygen consumption rate determined in Step 2, determine the discrete concentration distribution function of oxygen under the target operating condition, the equilibrium concentration distribution function, and the reaction source term. Then, using the evolution equation of the discrete concentration distribution function, perform collision and migration calculations on all grid points within the computational domain to update the discrete concentration distribution function of each grid point.

[0036] Step 8: Calculate the macroscopic concentration of oxygen at this time using the new discrete concentration distribution function of each grid point calculated in Step 7.

[0037] Step 9: Substitute the new oxygen concentration and temperature of each grid point in the computational domain obtained in Steps 6 and 8 into Step 2 to update the dynamic reaction rate, and repeat Steps 2, 3, 4, 5, 6, 7, and 8 to start the next round of iteration. When the change in macroscopic rate does not exceed the set threshold, the overall reaction transport process is determined to have reached a steady state, and the loop iteration is terminated.

[0038] Furthermore, the specific steps of step 1 are as follows: (1) A random algorithm is used to construct the diffusion layer structure of the fuel cell. First, a certain proportion of growth cores are randomly distributed in the computational domain. The target diameter, expansion rate and spatial distribution constraints are set so that the growth cores continuously expand and accumulate outward, thereby forming a diffusion layer structure within the target parameter range. Then, a gas flow layer, microporous layer and catalyst layer that match the size of the diffusion layer are generated to form an RGB vector image of the overall structure of the fuel cell cathode.

[0039] (2) The generated RGB vector image of the overall structure of the fuel cell cathode is converted into a binary image, and then converted into a corresponding 0 / 1 binary matrix array according to the spatial arrangement of the pixels. The resulting binary matrix array is then converted into a grid point distribution matrix proportionally, which serves as the basic computational domain of this simulation method, where each grid point is used as a basic computational unit.

[0040] Furthermore, the specific steps of step 2 are as follows: (1) Based on the basic computational domain determined in step 1, set the initial and boundary conditions for the simulation process. Specifically, set the initial temperature, water density, and oxygen concentration within the computational domain. Set constant oxygen concentration boundary, constant temperature boundary, and no-slip bounce boundary condition for water flow at the top of the computational domain (gas flow region). Set dynamic reaction boundary at the bottom of the computational domain (catalytic layer region). Set periodic boundaries for the left and right boundaries of the computational domain. Set the flow-solid boundary within the computational domain as a no-slip bounce boundary during the calculation of water flow and oxygen diffusion, and as a conjugate heat transfer boundary during the calculation of heat transfer.

[0041] (2) Set the selected reaction parameters, physical property parameters, and macroscopic parameters within the computational domain. Reaction parameters include: reference exchange current density. overpotential Transmission coefficient Reference oxygen concentration enthalpy of reaction Gas constant Faraday constant ,activation energy Physical properties include: viscosity thermal conductivity Specific heat capacity Eccentricity factor diffusion coefficient Flow-solid wetting angle Macroscopic parameters include: water density. Molar mass of water Quality and speed Oxygen concentration Molar mass of oxygen ,temperature ,pressure .

[0042] (3) Calculate the dynamic reaction rate under the target operating condition based on the set reaction parameters, physical property parameters and macroscopic parameters, specifically including: water generation rate. Oxygen consumption rate and heat generation rate .

[0043] Among them, water generation rate The corresponding oxygen reaction flux needs to be calculated first. Specifically, it can be determined based on the oxygen concentration in the catalytic region. and rate constant The specific calculation formula is as follows: ; rate constant It can be calculated using the following equation: ; ; Water reaction flux Oxygen reaction flux The calculation yielded: .

[0044] Correspondingly, water generation rate and oxygen consumption rate Water reaction flux oxygen reaction flux The calculation yielded: ; ; in The fracture volume of the microporous layer structure is given. This is the effective reaction area; And the rate of heat generation It can be determined by the oxygen consumption rate Further calculations yielded the following: .

[0045] Furthermore, the specific steps of step 3 are as follows: (1) Based on the reaction parameters, physical property parameters, macroscopic parameters, and water generation rate determined in step 2, the initial discrete velocity distribution function of water under the target operating condition is first constructed. Since the computational domain consists of multiple discrete lattice points, each lattice point participates in the evolution calculation as an independent computational unit, and its local particle distribution state is characterized by a discrete distribution function. By updating and statistically analyzing the distribution functions of all lattice points in the computational domain, the overall distribution characteristics of water can be obtained. Furthermore, after discretizing the continuous particle velocity in space into a finite number of velocity directions, the continuous distribution function is expressed as a discrete velocity distribution function corresponding to the discrete velocity directions. The velocity distribution function in equilibrium is defined as the equilibrium velocity distribution function. The initial discrete velocity distribution function of the water is... The specific calculation formula is as follows: ; ; In the formula, The index indicates the direction of the discrete velocity. Indicates the first The specific value of the weight factor corresponding to each discrete velocity direction varies depending on the selection of the discrete velocity set, but the sum of the weight factors of all discrete velocities must equal 1. The spatial vector distribution representing discrete velocities must also satisfy the conditions obtained from all discrete velocities. The modulus is 1; , This represents the lattice velocity, and its value needs to be converted to the actual speed of sound. When the actual speed of sound is 332.532 m / s, .

[0046] (2) Based on the boundary conditions determined in step 2, the bottom of the computational domain is a reactive boundary, where water will be dynamically generated; the left and right sides of the computational domain are periodic boundaries; and the top of the computational domain is set as a no-slip bounce boundary. Solid regions within the computational domain are also set as no-slip bounce boundaries, meaning that when a grid point is identified as a solid region, the conventional collision process is not performed, but a bounce occurs at that grid point. Based on the differences in the above boundary conditions, the discrete velocity distribution function of water in different grid points needs to be assigned corresponding values. (3) By performing migration and collision calculations on all lattice points within the computational domain, the updated discrete velocity distribution function for each lattice point is obtained. Since the time and spatial discretization in this model are not independent but interconnected through the discrete velocities of the particles, particle motion can be divided into two stages: migration and collision. Between two adjacent time steps, a particle migrates from one lattice node to its neighboring lattice point and collides with other lattice points upon reaching that point. These migration and collision processes together determine the evolution of the distribution function, and its evolution equation can be expressed as: ; In the formula, by adjusting the position ,time By solving for the discrete velocity distribution function and its equilibrium velocity distribution function, the velocity distribution at adjacent positions can be obtained. ,time The discrete velocity distribution function is given by, where, For time step; item The discrete force terms representing different components specifically include the gas-liquid phase transition force term and the fluid-solid wettability force term for water; It represents the discrete reaction source term, specifically including the rate of water formation in the dynamic reaction; This indicates discrete additional terms, specifically including mechanical consistency and surface tension.

[0047] in A matrix represented by a specific numerical value: ; It is the identity matrix. It is a diagonal (relaxed) matrix, represented as: ; in , , .

[0048] in This is a preset constant (an empirical parameter). Adjusting this value can reduce the spurious velocity (false velocity) of the model. Viscosity of water at different lattice points Related, expressed as: .

[0049] Can be determined by density velocity vector ,as well as Velocity components in different directions , Specifically, it means: ; in These are preset constants (belonging to free parameters), while It can be derived from the force term vector ,as well as Force components in different directions , Specifically, it means: ; in Specifically, it includes force terms related to water vapor-liquid phase change and fluid-solid wettability, expressed as: ; in This indicates the strength of the gas-liquid phase interaction of water. This represents the potential function when the lattice points are fluid regions with water distribution. The specific calculation is as follows: ; in Indicates the spatial step size. Representing macroscopic pressure, the above method can effectively represent the gas-liquid phase transition state of water. When the lattice point is a solid region and is not in contact with any fluid region, ; This refers to the situation where a lattice point is a solid region but in contact with the surrounding fluid region. It is necessary to calculate the average potential function of the surrounding fluid region. The wettability of the fluid-solid surface can be characterized using the above method, specifically calculated as follows: ; ; in and It is a constant when the fluid-solid wetting angle is 90°. or ;when ,and It can be increased when further reduction is needed. ,when ,and It can be improved when further enlargement is needed. . For indicator functions, Represents solids. Represents fluid.

[0050] and Specifically, it includes the rate of water formation during dynamic reactions. Specifically, it is expressed as: ; and Specifically, this includes mechanical consistency and surface tension, which are expressed as follows: ; ; in It is a mechanical stability adjustment parameter that can be used to modify mechanical stability conditions. This is a surface tension adjustment parameter, which can be used to adjust surface tension. To adjust the auxiliary parameters for mechanical stability, by setting... and The value of can be further obtained. .

[0051] Furthermore, the specific steps of step 4 are as follows: The discrete velocity distribution function of each grid point is evolved iteratively to obtain a new discrete velocity distribution function. Combined with the spatial vector distribution of discrete velocities The macroscopic density of water at different grid points within the computational domain is obtained. and macro speed The specific calculation is as follows: ; .

[0052] Then, based on the macroscopic density of each fluid region grid point... and temperature The macroscopic pressure at each grid point is further calculated using the equation of state. The specific calculation is as follows: ; ; ; ; in As the eccentricity factor, The critical temperature. The critical pressure. This is the gas constant. When the lattice points are solid regions, .

[0053] Furthermore, the specific steps of step 5 are as follows: Based on the calculated macroscopic density of water Further calculations yielded the new viscosity at each fluid region's grid point. Specific heat capacity and thermal conductivity The specific calculation is as follows: ; ; ; in , , , Let represent the density, viscosity, specific heat capacity, and thermal conductivity of water vapor under saturated conditions, respectively. , , , Let these represent the density, viscosity, specific heat capacity, and thermal conductivity of liquid water in its saturated state, respectively. For the lattice points in the solid region, the density, specific heat capacity, and thermal conductivity are fixed at... , , .

[0054] Furthermore, the specific steps of step 6 are as follows: (1) Based on the boundary conditions determined in step 2, the bottom of the computational domain has a reactive boundary, which will have a dynamic heat generation process; the left and right sides of the computational domain are periodic boundaries; and the top of the computational domain is set as a constant temperature boundary. The heat transfer process between the fluid-solid interface in the computational domain adopts a conjugate heat transfer boundary. (2) Based on the heat generation rate determined in step 2, the macroscopic density, macroscopic velocity, macroscopic pressure of water determined in step 4, and the specific heat capacity and thermal conductivity determined in step 5, respectively, substitute the heat conduction term, convection term, and volume work term and reaction source term controlling the latent heat of vaporization in the temperature evolution equation. In the specific solution process, first express the temperature evolution equation in a first-order form with respect to time: ; in, This represents the temperature evolution operator within the fluid region. Represented as: ; Among them, items For convection terms, term For thermal conductivity, the term To control the volume work term for latent heat of vaporization, the term... This is the reaction source term. In the solid region, only the thermal conductivity term is considered; therefore, in the solid region… Represented as: .

[0055] In the temperature evolution calculation, the fourth-order Runge-Kutta method is used to perform collision migration traversal of all grid points to obtain the new macroscopic temperature at each grid point, specifically expressed as follows: ; ; ; ; ; in The new macroscopic temperature at each grid point is calculated for the next time step.

[0056] Furthermore, the specific steps of step 7 are as follows: (1) Based on the boundary conditions determined in step 2, the bottom of the computational domain has a reaction boundary, which will have a dynamic oxygen consumption process. The left and right sides of the computational domain are periodic boundaries, the top of the computational domain is set as a constant oxygen concentration boundary, and the solid region in the computational domain is set as a no-slip rebound boundary condition. Based on the water distribution calculated in step 4, the fluid regions covered by water in all grid points are also set as no-slip rebound boundary conditions, while the fluid regions not covered by water in the grid points are used as new effective oxygen diffusion regions.

[0057] (2) Based on the reaction parameters, physical property parameters, macroscopic parameters and oxygen consumption rate determined in step 2, determine the discrete concentration distribution function of oxygen under the target operating condition. The equilibrium concentration distribution function of each grid point is also discretized into the same number of discrete equilibrium concentration distribution functions. When it is the initial state, the corresponding discrete concentration distribution function and equilibrium discrete concentration distribution function Represented as: ; ; In the formula, This represents the index of the discrete concentration direction. For the first The weighting factors corresponding to each discrete concentration direction, and the sum of the weighting factors of all discrete quantities must equal 1; Subsequently, using the evolution equation of the discrete concentration distribution function, collision and migration calculations are performed on all grid points in the computational domain to obtain the updated discrete concentration distribution function for each grid point. The specific evolution equation is expressed as follows: ; In the formula, by adjusting the position Location, Time By solving for the discrete concentration distribution function and its equilibrium concentration distribution function, we can obtain its distribution at adjacent locations. Location, Time Discrete concentration distribution function under; terms Represents the reaction source term in the oxygen diffusion process, where The rate of oxygen consumption; The relaxation time is related to the oxygen diffusion coefficient. Directly related, can be represented as: ; in Represents the lattice correlation coefficient; This is the weighting factor corresponding to the 0th discrete concentration direction.

[0058] Furthermore, the specific steps of step 8 are as follows: The new discrete concentration distribution function for each grid point is calculated using step 7. Substituting these values ​​into the following formula, the new macroscopic oxygen concentration can be calculated. Specifically, it is expressed as: .

[0059] Furthermore, the specific steps of step 9 are as follows: When the process reaches the [number]th In the next iteration, all grid points are traversed until the equation is satisfied: ; In the formula, This is the steady-state threshold for macroscopic velocity changes; when this condition is met, the entire process is considered to have reached a steady state, and the iterative cycle terminates. The output at this point includes: steady-state water distribution image, steady-state temperature distribution image, steady-state oxygen concentration distribution image, and average reaction rate change curves at different time steps.

[0060] Example 2 like Figure 1 The diagram shows a flowchart of a fuel cell cathode pore-scale simulation method that couples water-gas-liquid phase change heat transfer, oxygen transport, and dynamic reactions.

[0061] Step 1: A stochastic algorithm is used to construct the carbon fiber distribution characteristics in the diffusion layer, controlling the parameter ranges of its thickness, porosity, and fiber diameter. Combined with a simplified gas flow layer, microporous layer, and catalyst layer, the overall structure of the fuel cell cathode is formed. Then, the generated structure is binarized into a binary matrix corresponding to the pixels, and proportionally mapped to the lattice points within the computational domain. This matrix serves as the basic computational domain for simulation in this method, mainly including: (1) A random algorithm is used to construct the diffusion layer structure of the fuel cell cathode. First, a certain proportion of growth cores are randomly distributed in the computational domain. The target diameter, expansion rate and spatial distribution constraints are set so that the growth cores continuously expand and accumulate outward, thereby forming a diffusion layer structure within the target parameter range. Then, a gas flow layer, microporous layer and catalyst layer that match the size of the diffusion layer are generated to form an RGB vector image of the overall structure of the fuel cell cathode. The specific parameters are shown in Table 1.

[0062] Table 1. Main parameters of the computational domain (2) Use image processing methods such as MATLAB and OpenCV to convert the generated RGB vector image of the overall structure of the fuel cell cathode into a binary image. For example Figure 2 As shown, black represents the solid region of the fuel cell cathode, and green represents the flow region of the fuel cell cathode. The pixels are then arranged spatially to form a corresponding 0 / 1 binary matrix array. This binary matrix array is then proportionally converted into a grid distribution matrix, which serves as the basic computational domain for this simulation method, with each grid point acting as a basic computational unit.

[0063] Step 2: Based on the basic computational domain determined in Step 1, set the initial and boundary conditions for the simulation process. Specifically, this includes: reaction parameters (reference exchange current density, overpotential, transfer coefficient, reference oxygen concentration, enthalpy of reaction, gas constant, Faraday constant, activation energy), physical property parameters (viscosity, thermal conductivity, specific heat capacity, eccentricity factor, diffusion coefficient, fluid-solid wetting angle), and macroscopic parameters (water density, water molar mass, velocity, oxygen concentration, oxygen molar mass, temperature, pressure). This allows for the calculation of the dynamic reaction rate (water generation rate, oxygen consumption rate, and heat generation rate) under the target operating condition, mainly including: (1) Based on the basic computational domain determined in step 1, set the initial and boundary conditions for the simulation process. Specifically, set the initial temperature within the computational domain. water density and oxygen concentration As shown below: ; ; .

[0064] Subsequently, constant oxygen concentration boundary, constant temperature boundary, and no-slip bounce boundary conditions for water flow are set at the top of the computational domain for the gas flow region, and dynamic reaction boundary is set at the bottom of the computational domain for the catalytic layer region. The left and right boundaries of the computational domain are set as periodic boundaries, wherein the constant oxygen concentration and constant temperature boundaries are specifically set as follows: ; .

[0065] The fluid-solid boundary within the computational domain is set as a no-slip bounce boundary during the calculation of water flow and oxygen diffusion, and as a conjugate heat transfer boundary during the calculation of heat transfer. Specifically, the conjugate heat transfer boundary is defined as follows: ; ; in It is represented as a solid region but in contact with the surrounding fluid region. It is represented as a fluid region but in contact with the surrounding solid region. This is represented as a solid region that is not in contact with the surrounding fluid region. It is represented as a fluid region that does not contact the surrounding solid region.

[0066] (2) Set the selected reaction parameters, physical property parameters, and macroscopic parameters within the computational domain. The reaction parameters include: reference exchange current density. overpotential Transmission coefficient Reference oxygen concentration enthalpy of reaction Gas constant Faraday constant ,activation energy The physical properties include: viscosity. thermal conductivity Specific heat capacity Eccentricity factor diffusion coefficient Flow-solid wetting angle The macroscopic parameters include: the density of water. Molar mass of water ,speed Oxygen concentration Molar mass of oxygen ,temperature ,pressure The macroscopic parameters need to be continuously updated through evolutionary calculations, and the specific values ​​of the remaining parameters are shown in Table 2. Table 2 Basic Parameters (3) Calculate the dynamic reaction rate under the target operating condition based on the set reaction parameters, physical property parameters and macroscopic parameters, specifically including: water generation rate. Oxygen consumption rate and heat generation rate Among them, the oxygen consumption rate The corresponding oxygen reaction flux needs to be calculated first. Specifically, it can be determined based on the oxygen concentration in the catalytic region. and rate constant The specific calculation formula is as follows: ; rate constant It can be calculated using the following equation: ; .

[0067] Water reaction flux Oxygen reaction flux The calculation yielded: .

[0068] Correspondingly, water generation rate and oxygen consumption rate Water reaction flux oxygen reaction flux The calculation yielded: ; .

[0069] And the rate of heat generation It can be determined by the oxygen consumption rate Further calculations yielded the following: .

[0070] Step 3: Based on the reaction parameters, physical property parameters, macroscopic parameters, and water generation rate determined in Step 2, determine the discrete velocity distribution function, equilibrium velocity distribution function, phase transition related force term, fluid-solid related force term, reaction source term, and additional term for water under the target operating condition. Then, using the evolution equation of the discrete velocity distribution function, collision and migration calculations are performed on all lattice points in the computational domain to obtain the new discrete velocity distribution function for each lattice point, mainly including: (1) Based on the reaction parameters, physical property parameters, macroscopic parameters and water generation rate determined in step 2, construct the discrete velocity distribution function of water under the target operating condition. and equilibrium velocity distribution function The initial discrete velocity distribution function of the water is The density of gaseous water in step (2) Substituting the initial equilibrium distribution function, the initial velocity... If set to 0, then and It can be calculated using the following formula: .

[0071] Subsequently, during iterative processes at different times, the latest calculated macroscopic density of water was... Macro speed Substitute the equilibrium velocity distribution function Solve using the calculation formula: ; In the formula, The values ​​are set from 0 to 8, containing a total of 9 discrete quantities; These are weighting factors corresponding to different discrete velocity directions, with specific values ​​as follows: ; In the equilibrium velocity distribution function , This represents the lattice velocity, in this example. ; This represents the spatial vector distribution of discrete velocities, in this example. It can be represented as: .

[0072] (2) Based on the boundary conditions determined in step 2, the bottom of the computational domain is the reaction boundary, where a dynamic water generation process will occur, and the left and right sides of the computational domain are periodic boundaries. The specific calculation is as follows: ; .

[0073] The upper edge of the computational domain is set as a no-slip bounce boundary, and the solid regions within the computational domain are also set as no-slip bounce boundaries. This means that if a cell at a certain location is determined to be a boundary entity, the normal collision process is omitted, and the density will be bounced. The specific calculation is as follows: .

[0074] (3) By performing migration and collision calculations on all lattice points within the computational domain, the updated discrete velocity distribution function for each lattice point is obtained. Since the time and spatial discretization in this model are not independent but interconnected through the discrete velocities of the particles, particle motion can be divided into two stages: migration and collision. Between two adjacent time steps, a particle migrates from one lattice node to its neighboring lattice point and collides with other lattice points upon reaching that point. These migration and collision processes together determine the evolution of the distribution function, and its evolution equation can be expressed as: ; In the formula, by adjusting the position ,time By solving for the discrete velocity distribution function and its equilibrium velocity distribution function, the velocity distribution at adjacent positions can be obtained. ,time The discrete velocity distribution function is given below, where in this example... ;item The discrete force terms representing different components specifically include the gas-liquid phase transition force term and the fluid-solid wettability force term for water; It represents the discrete reaction source term, specifically including the rate of water formation in the dynamic reaction; This represents discrete additional terms, specifically including mechanical consistency and surface tension. A matrix represented by a specific numerical value: ; It is the identity matrix. It is a diagonal (relaxed) matrix, represented as: ; in , , In this example, .

[0075] in Viscosity of water at different lattice points Related, expressed as: .

[0076] Can be determined by density velocity vector ,as well as Velocity components in different directions , Specifically, it means: ; In this example, .

[0077] and It can be derived from the force term vector ,as well as Force components in different directions , Specifically, it means: ; in Specifically, it includes force terms related to water vapor-liquid phase change and fluid-solid wettability, expressed as: ; in This indicates the strength of the gas-liquid phase interaction of water. This represents the potential function when the lattice points are fluid regions with water distribution. The specific calculation is as follows: ; in Representing macroscopic pressure, the above method can effectively represent the gas-liquid phase transition state of water. When the lattice point is a solid region and is not in contact with any fluid region, When a lattice point is a solid region but in contact with the surrounding fluid region, then... It is necessary to calculate the average potential function of the surrounding fluid region. The wettability of the fluid-solid surface can be characterized using the above method, specifically calculated as follows: ; ; in and It is a constant when the fluid-solid wetting angle is 90°. or ;when ,and It can be increased when further reduction is needed. ,when ,and It can be improved when further enlargement is needed. In this example, since the fluid-solid wetting angle is set to 120°, therefore... .

[0078] and Specifically, it includes the rate of water formation during dynamic reactions. Specifically, it is expressed as: ; In this example, .

[0079] and Specifically, this includes mechanical consistency and surface tension, which are expressed as follows: ; ; in It can be used to adjust mechanical stability. It can be used to adjust surface tension, in this example, , .

[0080] Step 4: Calculate the macroscopic density, macroscopic velocity, and macroscopic pressure of the water at this point using the discrete velocity distribution function of all grid points calculated in Step 3. This mainly includes: The discrete velocity distribution function of each grid point is evolved iteratively to obtain a new discrete velocity distribution function. Combined with the spatial vector distribution of discrete velocities The macroscopic density of water at different grid points within the computational domain is obtained. and macro speed The specific calculation is as follows: ; .

[0081] Then, based on the macroscopic density of each fluid region grid point... and temperature The macroscopic pressure at each grid point is further calculated using the equation of state. The specific calculation is as follows: ; ; ; ; In this example , , , Thus, the critical pressure can be determined. When the lattice points are solid regions, .

[0082] Step 5: Update the new viscosity, specific heat capacity, and thermal conductivity of all lattice points using the macroscopic density of water calculated in Step 4. This mainly includes: Based on the calculated macroscopic density of water Further calculations yielded the new viscosity at each fluid region's grid point. Specific heat capacity and thermal conductivity The specific calculation is as follows: ; ; ; in , , , Let represent the density, viscosity, specific heat capacity, and thermal conductivity of water vapor under saturated conditions, respectively. , , , Let these represent the density, viscosity, specific heat capacity, and thermal conductivity of liquid water in its saturated state, respectively. The density, specific heat capacity, and thermal conductivity of the lattice points in the solid region are fixed at... , , The specific values ​​of the above parameters in this example are shown in Table 2.

[0083] Step 6: Substitute the heat generation rate calculated in Step 2 into the reaction source term of the temperature evolution equation. Substitute the macroscopic density, macroscopic velocity, and macroscopic pressure of water calculated in Step 4, and the specific heat capacity and thermal conductivity calculated in Step 5 into the heat conduction term, convection term, and volume work term related to latent heat of vaporization in the temperature evolution equation, respectively. Then, traverse all lattice points through a collision migration process to obtain the new macroscopic temperature at each lattice point, mainly including: (1) Based on the boundary conditions determined in step 2, the bottom of the computational domain has a reactive boundary, which will have a dynamic heat generation process; the left and right sides of the computational domain are periodic boundaries; and the top of the computational domain is set as a constant temperature boundary, specifically expressed as follows: ; ; .

[0084] For the heat transfer process between the fluid-solid interface in the computational domain, a conjugate heat transfer boundary is used, specifically represented as follows: ; ; in It is represented as a solid region but in contact with the surrounding fluid region. It is represented as a fluid region but in contact with the surrounding solid region. This is represented as a solid region that is not in contact with the surrounding fluid region. It is represented as a fluid region that does not contact the surrounding solid region.

[0085] (2) Based on the heat generation rate determined in step 2, the macroscopic density, macroscopic velocity, macroscopic pressure of water determined in step 4, and the specific heat capacity and thermal conductivity determined in step 5, respectively, substitute the heat conduction term, convection term, and volume work term and reaction source term controlling the latent heat of vaporization in the temperature evolution equation. In the specific solution process, first express the temperature evolution equation in a first-order form with respect to time: ; in, This represents the temperature evolution operator within the fluid region. Represented as: ; Among them, items For convection terms, term For thermal conductivity, the term To control the volume work term for latent heat of vaporization, the term... This is the reaction source term. In the solid region, only the thermal conductivity term is considered; therefore, in the solid region… Represented as: .

[0086] In the temperature evolution calculation, the fourth-order Runge-Kutta method is used to perform collision migration traversal of all grid points to obtain the new macroscopic temperature at each grid point, specifically expressed as follows: ; ; ; ; ; in The new macroscopic temperature at each grid point is calculated for the next time step.

[0087] Step 7: Based on the macroscopic density of water calculated in Step 4, determine the fluid regions not covered by water in all grid points as new effective diffusion regions. Based on the reaction parameters, physical property parameters, macroscopic parameters, and oxygen consumption rate determined in Step 2, determine the discrete concentration distribution function, equilibrium distribution function, and reaction source term of oxygen under the target operating condition. Then, using the evolution equation of the discrete concentration distribution function, perform collision and migration calculations on all grid points in the computational domain to update the discrete concentration distribution function of each grid point, mainly including: (1) Based on the boundary conditions determined in step 2, the bottom of the computational domain has a reaction boundary, which involves a dynamic oxygen consumption process. The upper side of the computational domain is set as a constant oxygen concentration boundary, specifically expressed as follows: ; .

[0088] The left and right sides of the computational domain are periodic boundaries, specifically represented as follows: ; .

[0089] Both the solid region and the fluid region covered by water in the computational domain are set to no-slip bounce boundary conditions, specifically as follows: .

[0090] (2) Based on the reaction parameters, physical property parameters, macroscopic parameters and oxygen consumption rate determined in step 2, determine the discrete concentration distribution function of oxygen under the target operating condition. The equilibrium concentration distribution function of each grid point is also discretized into the same number of discrete equilibrium concentration distribution functions. When it is the initial state, the corresponding discrete concentration distribution function and equilibrium discrete concentration distribution function Represented as: ; ; In the formula, The values ​​are set from 0 to 4, containing a total of 5 discrete quantities. For the first The weighting factors corresponding to each discrete concentration direction, where The specific value is The other specific values ​​are: .

[0091] Subsequently, using the evolution equation of the discrete concentration distribution function, collision and migration calculations are performed on all grid points in the computational domain to obtain the updated discrete concentration distribution function for each grid point. The specific evolution equation is expressed as follows: ; In the formula, by adjusting the position Location, Time By solving for the discrete concentration distribution function and its equilibrium concentration distribution function, we can obtain its distribution at adjacent locations. Location, Time Discrete concentration distribution function; This represents the spatial vector distribution of discrete concentrations, in this example. It can be represented as: .

[0092] item Represents the reaction source term in the oxygen diffusion process, where The rate of oxygen consumption; The relaxation time is related to the oxygen diffusion coefficient. Directly related, the specific values ​​in this example are shown in Table 2, and the calculation is as follows: ; In this example , .

[0093] Step 8: Calculate the macroscopic oxygen concentration at this point using the new discrete concentration distribution function calculated in Step 7 for each grid point. This mainly includes: The new discrete concentration distribution function for each grid point is calculated using step 7. Substituting these values ​​into the following formula, the new macroscopic oxygen concentration can be calculated. Specifically, it is expressed as: .

[0094] Step 9: Substitute the new oxygen concentration and temperature of each grid point within the computational domain obtained in Steps 6 and 8 into Step 2 to update the dynamic reaction rate, and repeat Steps 2, 3, 4, 5, 6, 7, and 8 to begin the next iteration. When the macroscopic rate change does not exceed a set threshold, the overall reaction transport process is considered to have reached a steady state, and the iteration cycle terminates. This mainly includes: When the process reaches the [number]th In the next iteration, all grid points are traversed until the equation is satisfied: ; When this condition is met, the entire process is considered to have reached a steady state, and the iterative loop terminates. The output at this point includes: a steady-state water distribution image, a steady-state temperature distribution image, a steady-state oxygen concentration distribution image, and oxygen consumption rate change curves at different time steps. Figure 3 , Figure 4 , Figure 5 , Figure 6 ).

[0095] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0096] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0097] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0098] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0099] In summary, this invention discloses a pore-scale simulation method for fuel cell cathodes that couples water-gas-liquid phase change heat transfer, oxygen transport, and dynamic reactions. This method utilizes a stochastic algorithm to generate actual diffusion layer structural features and combines simplified gas flow layer, micropore layer, and catalyst layer scale features to effectively characterize the overall structural properties of the fuel cell cathode. Based on determined initial and boundary conditions, discrete velocity evolution equations considering water-gas-liquid phase change flow and dynamic reaction water production processes are established; temperature evolution equations considering latent heat of water vaporization, fluid-solid conjugate heat transfer, and dynamic reaction heat production processes are established; and discrete concentration evolution equations considering the water's hindering effect on oxygen transport and dynamic oxygen consumption processes are established. By effectively coupling water-gas-liquid phase change heat transfer, oxygen transport, and dynamic reaction processes within the complex pore structure of the fuel cell cathode, information interaction and dynamic feedback between multiple physics fields at the pore scale are realized. The dynamic evolution processes of the velocity field, temperature field, concentration field, and reaction rate can be obtained in real time, thereby revealing the hydrothermal coupling effect and its microscopic mechanism within the fuel cell cathode. Once a steady state is reached, the method can output curves showing the water distribution, temperature distribution, oxygen concentration distribution, and reaction rate changes inside the fuel cell cathode. This method is designed for the pore-scale structural characteristics of real fuel cells and is suitable for studying the coupling relationships and information interaction mechanisms between water flow, oxygen diffusion, heat transfer, and dynamic reaction processes within complex pore structures. It can intuitively characterize the influence of microstructural changes in the fuel cell cathode on multiphysics reaction transport characteristics and its mechanism of action on actual reaction performance, thus providing important support for fuel cell structural optimization design and enhanced reaction performance.

[0100] The above description of the embodiments is provided to enable those skilled in the art to understand and use the invention. It will be apparent to those skilled in the art that various modifications can be made to these embodiments, and the general principles described herein can be applied to other embodiments without inventive effort. Therefore, the present invention is not limited to the above embodiments, and any improvements and modifications made by those skilled in the art based on the disclosure of the present invention without departing from the scope of the invention should be within the protection scope of the present invention.

Claims

1. A method for simulating the cathode pore size of a fuel cell that couples water-gas-liquid phase change heat transfer, oxygen transport, and dynamic reactions, characterized in that, Includes the following steps: Step 1: Construct the diffusion layer structure using a random algorithm and generate a gas flow layer, microporous layer, and catalyst layer that match the diffusion layer to form the overall structure of the fuel cell cathode; binarize the overall structure of the fuel cell cathode and map it proportionally to construct the computational domain; Step 2: Set initial and boundary conditions based on the computational domain constructed in Step 1. The initial conditions include reaction parameters, physical property parameters, macroscopic parameters, and the calculated dynamic reaction rate. The boundary conditions are as follows: the top of the computational domain is set as a constant oxygen concentration boundary, a constant temperature boundary, and a no-slip bounce boundary condition for water flow; the bottom of the computational domain is set as a dynamic reaction boundary; and the left and right boundaries of the computational domain are set as periodic boundaries. The flow-solid boundary within the computational domain is set as a no-slip bounce boundary during the calculation of water flow and oxygen diffusion, and as a conjugate heat transfer boundary during the calculation of heat transfer. Step 3: Based on the initial conditions set in Step 2, construct the discrete velocity distribution function and equilibrium velocity distribution function of water under the target working condition. Then, use the evolution equation to perform collision and migration calculations on all grid points in the computational domain to update the discrete velocity distribution function of each grid point. Step 4: Based on the results of Step 3, update the macroscopic density, macroscopic velocity, and macroscopic pressure of water at each grid point; Step 5: Based on the results of Step 4, update the viscosity, specific heat capacity, and thermal conductivity of all grid points; Step 6: Substitute the results of Step 2, Step 4 and Step 5 into the temperature evolution equation to perform collision and migration calculations on all grid points in the computational domain and update the macroscopic temperature of each grid point. Step 7: Determine the new effective diffusion region based on the macroscopic density of water calculated in Step 4; construct the discrete concentration distribution function and equilibrium concentration distribution function of oxygen under the target operating condition according to the initial conditions set in Step 2; then use the evolution equation to perform collision and migration calculations on all grid points to update the discrete concentration distribution function of each grid point. Step 8: Based on the results of Step 7, calculate and update the macroscopic oxygen concentration at each grid point; Step 9: Substitute the results of Step 6 and Step 8 into Step 2 to update the dynamic reaction rate, and repeat Step 2 to Step 8 for iteration; when the change in macroscopic velocity obtained by iteration does not exceed the set threshold, it is determined that a stable state has been reached, and the iteration is terminated.

2. The method for simulating the cathode pore size of a fuel cell coupled with water-gas-liquid phase change heat transfer, oxygen transport, and dynamic reaction as described in claim 1, is characterized in that... In step 2: The dynamic reaction rate includes: water generation rate Oxygen consumption rate and heat generation rate ; Wherein: oxygen consumption rate Based on oxygen reaction flux The calculated water generation rate Based on water reaction flux The calculated heat generation rate Based on oxygen consumption rate Calculated.

3. The method for simulating the cathode pore size of a fuel cell coupled with water-gas-liquid phase change heat transfer, oxygen transport, and dynamic reactions according to claim 1, characterized in that, In step 3: Equilibrium velocity distribution function for: ； In the formula: Indicates the first Each discrete velocity direction corresponds to a weight factor, and the sum of all weight factors is 1; This indicates the macroscopic density of water; Indicates macroscopic speed; The spatial vector distribution representing discrete velocities, and The modulus is 1; , This indicates the lattice velocity, which is converted based on the actual speed of sound. The evolution equation is: ； In the formula: Indicates the time step; Represents the discrete force terms for different components; Represents discrete reaction source terms; Indicates discrete additional terms; ； It is the identity matrix. It is a diagonal matrix.

4. The method for simulating the cathode pore size of a fuel cell coupled with water-gas-liquid phase change heat transfer, oxygen transport, and dynamic reaction as described in claim 1, is characterized in that... In step 4: Macroscopic density of water for: ； Macro speed for: ； In the formula: The discrete velocity distribution function of water; This indicates the macroscopic density of water; The spatial vector distribution representing discrete velocities, and The modulus is 1; Indicates the time step; Represents a force vector; When the lattice point is a fluid region, the macroscopic pressure is... According to macro density and temperature Calculated using state equations; When the lattice point is a solid region, the macroscopic pressure is... .

5. The method for simulating the cathode pore size of a fuel cell coupled with water-gas-liquid phase change heat transfer, oxygen transport, and dynamic reaction as described in claim 1, characterized in that, In step 5: When the lattice points are fluid regions, then: Viscosity for: ； Specific heat capacity for: ； thermal conductivity for: ； In the formula: This represents the density of vapor water under saturated conditions. This represents the density of liquid water in a saturated state. This represents the viscosity of vapor water under saturated conditions. This indicates the viscosity of liquid water under saturated conditions. This represents the specific heat capacity of vapor water under saturated conditions. This is expressed as the specific heat capacity of liquid water under saturated conditions. This represents the thermal conductivity of vapor water under saturated conditions. This represents the thermal conductivity of liquid water under saturated conditions. When the lattice points are solid regions, the density is... Specific heat capacity and thermal conductivity It is a fixed value.

6. The method for simulating the cathode pore size of a fuel cell coupled with water-gas-liquid phase change heat transfer, oxygen transport, and dynamic reaction as described in claim 1, is characterized in that... In step 6: The temperature evolution equation is: ； In the formula: Represents the temperature evolution operator; When the lattice point is a fluid region, the temperature evolution operator simultaneously considers the convection term, the heat conduction term, the volume work term controlling the latent heat of vaporization, and the reaction source term. When the lattice point is a solid region, the temperature evolution operator only considers the thermal conductivity term; The macroscopic temperature of each grid point is calculated and updated using the fourth-order Runge-Kutta method.

7. The method for simulating the cathode pore size of a fuel cell coupled with water-gas-liquid phase change heat transfer, oxygen transport, and dynamic reaction as described in claim 1, is characterized in that... In step 7: The equilibrium concentration distribution function is: ； In the formula: Indicates the first Each discrete concentration direction has a weighting factor, and the sum of all weighting factors is 1; Indicates the concentration of oxygen; The evolution equation is: ； In the formula: Indicates the time step; Indicates the rate of oxygen consumption; Indicates the relaxation time.

8. The method for simulating the cathode pore size of a fuel cell coupled with water-gas-liquid phase change heat transfer, oxygen transport, and dynamic reaction as described in claim 1, is characterized in that... In step 8: The macroscopic concentration of oxygen is: ； In the formula: This represents the discrete concentration distribution function of oxygen.

9. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores one or more programs, which can be executed by one or more processors to implement the steps of the fuel cell cathode pore size simulation method as described in any one of claims 1-8.

10. A fuel cell cathode pore-scale simulation system, characterized in that, The device includes a processor and a memory for executable instructions, wherein the processor executes the instructions to implement the steps of the fuel cell cathode pore size simulation method as described in any one of claims 1-8.

Images