A mesoscale simulation method for fuel cell diffusion layers

Abstract

The present invention belongs to the field of fuel cell technology, and specifically relates to a mesoscale simulation method for a fuel cell diffusion layer, comprising the following steps: constructing a calculation domain for the diffusion layer structure; calculating the electrochemical reaction rate under target operating conditions within the calculation domain; determining the discrete velocity distribution function of each grid point within the calculation domain; calculating the saturation, macroscopic pressure, and macroscopic velocity of water and oxygen; determining the mass transfer process control equation for each grid point within the calculation domain and calculating the oxygen concentration; determining the fluid-solid coupling heat transfer process control equation for each grid point within the calculation domain and calculating the macroscopic temperature; and judging whether the water saturation and macroscopic velocity meet a threshold. Compared with the prior art, the present invention solves the problem of the difficulty of simulating the precise coupling and mutual feedback between the mesoscale multi-physical fields of the prior art battery diffusion layer. This solution constructs a mesoscale fine simulation method for multiphase flow mass transfer, dynamic electrochemical reaction, and heat transfer within the flow layer, diffusion layer, and catalytic layer of a fuel cell.

Description

Technical Field

[0001] The present invention belongs to the technical field of fuel cells, and in particular relates to a mesoscale simulation method for a fuel cell diffusion layer. Background Art

[0002] Global climate change has shifted the focus from fossil fuels to clean and renewable energy, and new energy technologies are becoming a global solution. Fuel cells, in particular, as electrochemical energy conversion devices, offer advantages such as high efficiency, fast response, high power density, and zero environmental pollution. However, fuel cell performance depends critically on the diffusion and utilization of fuel and oxidant fluids at the anode and cathode, as well as the effective management of water generated at the electrodes. This is because sufficient water content maintains high proton conductivity, but under high loads, water may condense into droplets in the porous electrodes. If not adequately removed, liquid water will eventually accumulate in the diffusion and catalytic layers, hindering the diffusion of reactants and causing flooding, ultimately reducing overall fuel cell performance. Research suggests that optimizing the electrode diffusion layer can effectively mitigate flooding. The functions of the fuel cell diffusion layer include removing accumulated water from the fuel cell, providing a pathway for reactants to diffuse to the catalyst layer, providing structural support for the reaction zone, and facilitating electron and heat transfer. Further mitigation of flooding can be achieved by designing different structural features of the diffusion layer or modifying its surface properties.

[0003] However, existing simulation methods primarily focus on macroscopic or systemic research. For example, CN116467913A discloses a numerical simulation method for solid oxide fuel cells, establishing a macroscopic simulation method for coupled heat transfer, mass transfer, and electrochemical reactions within the stack, enabling rapid simulation of solid oxide fuel cells. Another example is the computer simulation method for mass, momentum, energy, and charge transfer in proton exchange membrane fuel cells disclosed in CN118472326A, which establishes a physical transfer method between different layers within the electrode, enabling simulation of the transport and migration characteristics of various fluids within the electrode. These approaches can all achieve fuel cell performance analysis and parameter control, but they are primarily based on continuously adjusting parameters to find patterns and then determine a more reasonable electrode design. This approach lacks a deep understanding of the interactions between multiphase flow and complex electrode structures, the coupling processes between heat transfer, mass transfer, and electrochemical reactions, and the microscopic mechanisms. Therefore, the design of the solution is often influenced by experience. Mesoscale research, on the other hand, requires simultaneous consideration of the interactions under multiple dynamic conditions, such as fluid dynamics, mass transfer, heat transfer, and electrochemical reactions, allowing for a deeper understanding of the mechanisms and influencing mechanisms that affect the electrochemical performance of fuel cells. However, due to the mutual coupling of gases, liquids, solids and electrochemical processes, their behavior is particularly complex at the mesoscale. Therefore, it is very difficult to develop a simulation method for the precise coupling and mutual feedback between the mesoscale multi-physical fields in the fuel cell diffusion layer.

[0004] At present, there are the following technical difficulties in the fine-scale simulation methods for the diffusion layer of fuel cells: (1) It is difficult to accurately characterize the details of the electrode microstructure, material composition and interface phenomena, and it is impossible to effectively establish the correlation between physical behaviors from micro to macro scales. (2) It is necessary to simultaneously consider the interaction between multiple physical fields such as multiphase fluid-solid interaction, heat transfer, mass transfer and electrochemical reactions, but the current simulation methods are difficult to accurately couple these complex processes. (3) At present, the electrochemical modeling of fuel cells is mostly based on a constant reaction rate, but actual fuel cells operate under dynamic load and temperature fluctuation conditions. How to effectively express the dynamic water management and oxygen transfer behavior in the diffusion layer is very important.

[0005] Therefore, at this stage, there is an urgent need for a sophisticated simulation method that is applicable to the research on multi-physics field coupled water management in the diffusion layer of fuel cells. Summary of the Invention

[0006] The purpose of the present invention is to provide a mesoscale simulation method for a fuel cell diffusion layer in order to solve at least one of the above problems, so as to solve the problem in the prior art that due to the mutual coupling of gases, liquids, solids and electrochemical processes, their behaviors are particularly complex at the mesoscale, and the simulation method of accurate coupling and mutual feedback between the mesoscale multi-physical fields of the battery diffusion layer is very difficult. This solution constructs a mesoscale fine simulation method for multiphase flow mass transfer, dynamic electrochemical reaction and heat transfer inside the flow layer, diffusion layer and catalytic layer of the fuel cell based on the structural characteristics of the fuel cell diffusion layer, combined with the simplified scale characteristics of other functional layers (flow layer and catalytic layer). This method is oriented towards the heterogeneous electrode structure characteristics of the real diffusion layer material, and is suitable for studying the overall dynamic electrochemical reaction performance of the fuel cell based on the diffusion layer and in combination with other functional layers, as well as the interaction and transmission mechanism of the multi-component flow, heat transfer and mass transfer processes in the porous electrode contained therein.

[0007] The purpose of the present invention is achieved through the following technical solutions:

[0008] A mesoscale simulation method for a fuel cell diffusion layer comprises the following steps:

[0009] Step 1: Construct a calculation domain of a diffusion layer structure, wherein the top of the calculation domain is the flow layer and the bottom of the calculation domain is the catalytic layer;

[0010] Step 2: Calculate the electrochemical reaction rate under the target working conditions within the calculation domain;

[0011] Step 3: Determine the discrete velocity distribution function of each grid point in the computational domain;

[0012] Step 4: Calculate the water and oxygen saturation, macroscopic pressure, and macroscopic velocity;

[0013] Step 5: Determine the mass transfer process control equation for each grid point in the computational domain and calculate the oxygen concentration;

[0014] Step 6: Determine the governing equations of the fluid-solid coupled heat transfer process at each grid point in the computational domain and calculate the macroscopic temperature;

[0015] Step 7: Determine whether the water saturation and macro velocity meet the set thresholds:

[0016] If the conditions are met, it is determined that a stable state has been reached, and the water saturation change curve, the migration distribution image of oxygen and water in the electrode at steady state, the velocity gradient image and the temperature distribution image are output;

[0017] If not, it is determined that the steady state has not been reached, and the macroscopic temperature and oxygen concentration calculated in steps 5 and 6 are substituted into step 2 for iteration.

[0018] Preferably, the mesoscale simulation method is carried out under the following settings: in the initial state, the flow layer, diffusion layer and catalytic layer are in an oxygen-saturated state; a constant oxygen concentration injection and a constant temperature distribution are maintained at the top of the calculation domain, and an electrochemical reaction with dynamic changes in oxygen concentration and temperature occurs at the bottom of the calculation domain; the left and right sides of the calculation domain are set as periodic boundary conditions, and the solid wall surface within the calculation domain adopts a no-slip rebound boundary.

[0019] Preferably, step 1 comprises the following steps:

[0020] (1) Slices of the porous media structure of the diffusion layer of the electrode material in the fuel cell are collected and image-scanned to generate a flow layer and a catalytic layer that match the size of the diffusion layer, forming an RGB vector image of the overall electrode structure;

[0021] (2) Distinguish the structural features of different levels in the electrode structure and convert the RGB vector image into a binary image;

[0022] (3) Convert the binary image into a binary matrix array according to the arrangement of pixel points as the calculation domain.

[0023] Preferably, step 2 comprises the following steps:

[0024] (1) Determine the initial distribution of oxygen and water in the computational domain;

[0025] (2) Setting the macroscopic parameters and basic physical property parameters of the computational domain;

[0026] (3) Calculate the oxygen consumption rate and water generation rate in the electrochemical reaction under the target operating conditions:

[0027] Oxygen consumption rate R 氧 (x,t):

[0028]

[0029] Reaction rate constant k elec (x,t) is calculated according to the following equation:

[0030]

[0031] Where i0 is the exchange current density, F is the Faraday constant, C 氧,ref is the reference oxygen concentration, α is the transfer coefficient, R is the universal gas constant, T(x,t) is the macroscopic temperature, and η is the overpotential

[0032] Water production rate R 水 (x,t):

[0033] R 水 (x,t)=-2.0*R 氧 (x,t).

[0034] Preferably, step 3 comprises the following steps:

[0035] (1) Determine the initial discrete velocity distribution functions of oxygen and water

[0036]

[0037] Where k represents the serial number of the discrete quantity, a represents the serial number of different fluids in the computational domain; w k are weight factors for different discrete quantities; c represents the grid velocity; represents the spatial distribution of discrete quantities; ρ a Indicates the macroscopic density of different phase fluids; Indicates the macroscopic velocity of different phase fluids;

[0038] (2) Constructing the interaction forces between fluids and between fluids and walls

[0039]

[0040] in,

[0041]

[0042] In the subscripts, a and β represent different fluids; for example, if a represents water, β represents oxygen; if a represents oxygen, β represents water;

[0043] Among them, ψ a (x,t)=1-exp(-ρ a (x,t)), G 氧-水 Indicates the strength of oxygen-water interaction;

[0044]

[0045] in, represents the interaction strength between the fluid phase a and the solid phase s, s(x) is an indicator function, s(x) = 1 represents solid, s(x) = 0 represents fluid;

[0046]

[0047] (3) Establish the evolution equation of the discrete distribution function considering the dynamic electrochemical reaction rate:

[0048]

[0049] Among them, [(M a ) -1 Λ a M a ] represents the relaxation time, represents the force terms of different phase fluids, represents the source term of the electrochemical reaction of different phase fluids, and Δt is the time step;

[0050]

[0051] Λ a is the diagonal relaxation matrix;

[0052]

[0053] R a (x, t) represents the electrochemical reaction rate, including the oxygen consumption rate and water production rate.

[0054] Preferably, Λ a At least one

[0055] Among them, relaxation time and fluid viscosity υ a satisfy:

[0056]

[0057] Λ a middle According to the specific simulation process, it is taken as a constant value through iteration, or as a value related to the fluid viscosity υ a Related values ​​(corresponding to different working conditions), where at least one Need to be consistent with the fluid viscosity a Related.

[0058] Preferably, in step 4:

[0059] Macro speed

[0060]

[0061] in, represents the initial relaxation time, and a represents the sequence number of different fluids in the computational domain;

[0062] represents the discrete velocity distribution function, k represents the sequence number of the discrete quantity;

[0063] Represents the spatial distribution of discrete quantities;

[0064]

[0065] Water saturation Sat 水 (x,t):

[0066]

[0067] Among them, m 水 represents the molar mass of water;

[0068] Oxygen saturation: Sat 氧 (x,t):

[0069]

[0070] Among them, m 氧 represents the molar mass of oxygen;

[0071] Macroscopic pressure P(x,t):

[0072]

[0073] in, c represents the grid velocity; ρ 氧 (x,t) represents the macroscopic density of oxygen; ρ 水 (x, t) represents the macroscopic density of water; G 氧-水 represents the interaction strength between oxygen and water; ψ 氧 (x, t) represents the potential function of oxygen; ψ 水 (x,t) represents the potential function of water.

[0074] Preferably, in step 5:

[0075] The mass transfer process control equation is:

[0076]

[0077] Among them, R 氧(x, t) is the oxygen consumption rate in the electrochemical reaction rate obtained in step 2, is the macroscopic velocity obtained in step 4, D 氧 is the diffusion coefficient of oxygen, C 氧 (x, t) is the oxygen concentration.

[0078] Preferably, in step 6:

[0079] The governing equation for the fluid-solid coupled heat transfer process is:

[0080]

[0081] Where T(x,t) represents the macroscopic temperature;

[0082] S T (x, t) = ΔH·R 水 (x,t) / (ρ all (x,t)c p (x, t)), ΔH represents the heat of reaction, R 水 (x, t) is the water generation rate in the electrochemical reaction rate obtained in step 2, ρ all (x, t) represents the density of different grid points in the computational domain, c p (x, t) represents the constant-pressure specific heat capacity of the multi-component fluid and solid at different lattice points;

[0083] In the flow area:

[0084] ρ all (x,t)=ρ 水 (x,t)+ρ 氧 (x,t),ρ 水 (x,t) represents the macroscopic density of water, ρ 氧 (x, t) represents the macroscopic density of oxygen;

[0085] represents the specific heat capacity of water at constant pressure, represents the specific heat capacity of oxygen at constant pressure; Sat 氧 (x, t) represents oxygen saturation, Sat 水 (x, t) represents water saturation;

[0086] λ(x,t)=n v (x,t)*λ 氧 +(1.0-n v (x,t))*λ 水 ,λ 水 represents the thermal conductivity of water, λ 氧 represents the oxygen thermal conductivity;

[0087] In solid regions:

[0088] ρ all (x,t)=ρ 固体 (x,t),ρ 固体 (x, t) represents the density value of solid particles;

[0089] represents the specific heat capacity of solid at constant pressure;

[0090] λ(x,t)=λ 固体 ,λ 固体 represents the thermal conductivity of solid.

[0091] Preferably, in step 7:

[0092] judge:

[0093] |Sat 水 (x,n)-Sat 水 (x,n-1)|<ε1;

[0094]

[0095] Among them, Sat 水 (x,n) represents water saturation, represents the macro velocity, ε1 is the steady-state threshold of water saturation, and ε2 is the steady-state threshold of macro velocity.

[0096] Preferably, in step 7, the output result is obtained by calculating the water saturation and oxygen saturation values ​​of all grid points in the domain, the velocity gradient values ​​in the domain in different spatial directions, the temperature distribution of the domain in different spatial directions, and the pressure distribution of the domain in different spatial directions.

[0097] Compared with the prior art, the present invention has the following beneficial effects:

[0098] 1. This method uses the structural characteristics of the diffusion layer of real electrode materials and combines them with the flow layer and catalytic layer of the same scale characteristics to carry out simulation research. It can more realistically reflect the complex structural characteristics inside the fuel cell electrode, thereby improving the accuracy and authenticity of the simulation method.

[0099] 2 This method combines multiple physical parameters such as the migration and distribution characteristics of water and oxygen, flow rate, pressure, concentration, temperature, and electrochemical reaction rate to achieve high-precision simulation of the multi-physical field coupled evolution process of oxygen and water in the diffusion layer. It not only studies the microscopic reaction mechanism and promotion mechanism inside the fuel cell, but also represents the overall macroscopic parameter changes of the fuel cell.

[0100] 3. This method establishes a fuel cell dynamic water management simulation method to simulate the generation and migration of water and its actual impact on the dynamic electrochemical reaction of the fuel cell. It can help optimize the design of the diffusion layer, alleviate water flooding, and improve the overall performance of the fuel cell. BRIEF DESCRIPTION OF THE DRAWINGS

[0101] Figure 1 Schematic diagram of the process of the fuel cell diffusion layer mesoscale fine simulation method of the present invention;

[0102] Figure 2 This is an image of the electrode structure of the fuel cell in Example 1;

[0103] Figure 3 is the discrete velocity model of each grid point in Example 1;

[0104] Figure 4 The water saturation variation curve of the fuel cell diffusion layer and the fluidized layer in Example 1;

[0105] Figure 5 This is an image of the migration and distribution of oxygen and water in the electrode of the fuel cell in steady state in Example 1;

[0106] Figure 6 This is the velocity gradient image in the electrode of the fuel cell in steady state in Example 1;

[0107] Figure 7 This is an image of the temperature distribution inside the electrode of the fuel cell in Example 1 when it is in steady state. DETAILED DESCRIPTION

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

[0109] To overcome the shortcomings of the existing technology, the present invention constructs a mesoscale fine simulation method for multiphase flow and mass transfer, dynamic electrochemical reactions, and heat transfer within the flow, diffusion, and catalytic layers of a fuel cell based on the structural characteristics of the fuel cell diffusion layer and combined with simplified scale characteristics of other functional layers (flow layer and catalyst layer). This method is targeted at the heterogeneous electrode structure characteristics of the real diffusion layer material and is suitable for studying the overall dynamic electrochemical reaction performance of the fuel cell based on the diffusion layer and combined with other functional layers within the porous electrode, as well as the interaction and transmission mechanism of the multi-component flow, heat, and mass transfer processes contained therein. Specifically, it includes the following steps:

[0110] Step 1: Using the structural characteristics of the diffusion layer of the real electrode material, combined with the simplified flow layer and catalyst layer, an RGB vector image of the fuel cell structure is formed. This image is then binarized into a binary matrix of corresponding pixels and an iterative step size is set to convert the pixels into corresponding grid points of equal proportion, which serves as the basic calculation domain.

[0111] Step 2: In the basic calculation domain, the electrochemical reaction rate of oxygen and water under the target working condition is calculated based on the initial distribution of oxygen and water in the calculation domain, macroscopic parameters (density and viscosity of water and oxygen, oxygen concentration, velocity, temperature, pressure), and basic physical properties (interaction strength coefficient between oxygen and water, interaction strength coefficient between different phase fluids and the wall, Faraday constant, exchange current density, universal gas constant, thermal conductivity, diffusion coefficient, constant-pressure specific heat capacity, and overvoltage).

[0112] Step 3: In the basic computational domain, based on the initial distribution, macroscopic parameters, basic physical properties, and electrochemical reaction rates of oxygen and water determined in Step 2, determine the initial discrete velocity distribution functions of oxygen and water under the target operating conditions. Then, by performing migration and collision on the evolution equation of the discrete velocity distribution function, calculate and traverse all grid points to obtain a new discrete velocity distribution function for each grid point.

[0113] Step 4: Calculate the saturation of water and oxygen, macroscopic pressure and macroscopic velocity at each grid point using the new discrete velocity distribution function calculated in step 3;

[0114] Step 5: Substitute the basic physical properties and the calculated electrochemical reaction rate of water from step 2 into the source term of the mass transfer process control equation. Substitute the macroscopic velocity calculated from step 4 into the mass transfer process control equation. Then, traverse all grid points through the migration collision process to obtain the new oxygen concentration at each grid point.

[0115] Step 6: Substitute the basic physical properties and the calculated electrochemical reaction rate of water from step 2 into the source term of the heat transfer process control equation. Substitute the macroscopic velocity calculated from step 4 into the fluid-solid coupling heat transfer process control equation and traverse all grid points through the migration collision process to obtain the new macroscopic temperature of each grid point.

[0116] Step 7: Substitute the new macroscopic temperature and oxygen concentration calculated in steps 5 and 6 into step 2 to obtain a new electrochemical reaction rate, and repeat steps 2, 3, 4, 5, and 6 to perform a new round of iterative calculations. When the changes in the water saturation and velocity gradient values ​​do not exceed the set thresholds, it is determined that the overall electrochemical reaction of the fuel cell has reached a stable state, and the iteration cycle ends.

[0117] Furthermore, the specific steps of step 1 are as follows:

[0118] (1) Collect the porous medium structure slices of the actual electrode material diffusion layer, use image representation and processing methods to scan the electrode slices, generate the flow layer and catalytic layer that match the diffusion layer size, and form an RGB vector image of the overall electrode structure;

[0119] (2) Distinguish the structural features of different levels of electrodes and convert the electrode structure RGB vector image into a binary image;

[0120] (3) The binary image is converted into a corresponding 0, 1 binary matrix array according to the arrangement position of the pixel points, and the iterative step size is set to convert the pixel points into a matrix array with the corresponding grid point distribution in equal proportion, which serves as the basic calculation domain of this simulation method, in which each grid point is regarded as a basic calculation unit.

[0121] Furthermore, the specific steps of step 2 are as follows:

[0122] (1) Determine the initial distribution state of oxygen and water in the basic calculation domain of the electrode. Specifically, the electrode flow layer, diffusion layer and catalytic layer are initially set to an oxygen saturation state. A constant oxygen concentration is injected through the flow layer at the top inlet of the flow layer, and diffuses into the diffusion layer. Finally, it reaches the catalytic layer to undergo redox reaction, consume oxygen and generate water. The temperature of the top inlet of the flow layer is also kept constant.

[0123] (2) Set the macroscopic parameters and basic physical parameters selected in the calculation domain, including the density of water and oxygen ρ 水 , ρ 氧 and viscosity υ 水 、υ 氧 , oxygen concentration C 氧 , temperature T; and basic physical parameters: interaction strength coefficient G between oxygen and water 氧-水 , the interaction strength coefficient between oxygen, water and the wall Faraday constant F, exchange current density i0, universal gas constant R, transfer coefficient α, reference oxygen concentration C 氧,ref , thermal conductivity λ 水 ,λ 氧 ,λ 固体 , specific heat capacity at constant pressure Overpotential η and diffusion coefficient:

[0124] (3) Calculate the oxygen consumption rate R in the electrochemical reaction based on the set macroscopic parameters and basic physical parameters 氧 (x, t) and the water generation rate R 水 (x,t). R 氧 (x, t) can be calculated based on the local oxygen concentration C at the catalytic site. 氧 and the reaction rate constant k elec (x, t) is determined, and the specific calculation formula is as follows:

[0125]

[0126] Reaction rate constant k elec (x,t) can be calculated using the following equation:

[0127]

[0128] Where i0 is the exchange current density, F is the Faraday constant, C 氧,ref is the reference oxygen concentration, α is the transfer coefficient, R is the universal gas constant, T(x,t) is the macroscopic temperature distribution, and η is the overpotential. Based on the above parameter values, the oxygen reaction consumption rate R can be determined. 氧 (x, t). The reaction rate of water production R 水 (x, t) can be calculated based on the oxygen consumption rate R 氧 (x, t) is obtained, and the specific calculation formula is as follows:

[0129] R 水 (x,t)=-2.0*R 氧 (x,t)

[0130] Furthermore, the specific steps of step 3 are as follows:

[0131] (1) According to the initial distribution, macroscopic parameters, basic physical parameters and electrochemical reaction rate of oxygen and water determined in step 2, the initial discrete distribution function of oxygen and water under the target working condition is determined; since each grid point is used as a basic calculation unit, the total fluid distribution in the calculation domain is obtained by traversing and calculating the particle distribution of all grid points in the entire basic calculation domain, and the particle distribution of different grid points is represented by the distribution function. When the particle distribution of a certain grid point is in equilibrium, the distribution function is called the equilibrium distribution function; since microscopic particles are constantly performing irregular thermal motion, the velocity of microscopic particles is continuous, and its velocity is infinite-dimensional in the phase space. However, the details of the particle motion do not significantly determine the macroscopic motion of the fluid. Therefore, the particle velocity is simplified to a finite-dimensional velocity space in the phase space, and then the continuous distribution function is also discretized into a discrete distribution function with the same number of discrete velocities. Then the equilibrium distribution function of each lattice point will also be discretized into the same number of discrete equilibrium distribution functions as the discrete velocity The initial discrete distribution functions of oxygen and water are The macroscopic density ρ of the different phase fluids in step (2) 水 , ρ 氧 , macro speed Substitute the equilibrium distribution function Find the formula:

[0132]

[0133] Where k represents the serial number of the discrete quantity (for example, if it is discretized into 9 quantities, k is 0 to 8), a represents the serial number of different fluids in the computational domain (since this model studies oxygen and water, the range of a is 0 to 1); w k It is the weight factor of different discrete quantities. The specific value is different for different discrete quantities, but the sum of the weight factors of all discrete quantities must be equal to 1; c 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, c = 1; Represents the spatial distribution of discrete quantities in the form of vector coordinates. The coordinate values ​​of are also different, but must satisfy The modulus is equal to 1, ρ a and They represent the macroscopic density of different phase fluids and the macroscopic velocity of the computational domain respectively;

[0134] (2) According to different working conditions, the boundary conditions of the reaction inside the fuel cell diffusion layer, flow layer and catalytic layer are set; the top of the calculation domain is set as the inlet, and a constant oxygen concentration is injected through the flow layer, and diffuses in the diffusion layer and finally reaches the catalytic layer to undergo redox reaction. The bottom of the calculation domain is the catalytic layer. Since the catalytic sites are distributed at the nanometer scale, and the diffusion layer is distributed at the micrometer scale, the size of the catalytic sites is very small in comparison, so the structure of the catalytic layer can be simplified, and the distribution area of ​​the catalytic reaction sites can be set instead according to the size of the actual cracks in the catalytic layer. Cracks in the catalytic layer are distributed at the bottom of the calculation domain, and dynamic electrochemical reactions occur at the cracks. The left and right sides of the calculation domain are set as periodic boundary conditions. The solid area in the calculation domain is subjected to boundary processing. The present invention adopts a no-slip rebound boundary, that is, the cell at a certain position is determined as a boundary entity, and the normal collision process needs to be omitted, and the density of the fluid at this location will be rebounded. According to the above different boundary conditions, the values ​​of the different discrete velocity distribution functions of the particles need to be modified accordingly;

[0135] (3) By calculating the evolution equation during the migration and collision process, all grid points are traversed to obtain a new discrete distribution function for each grid point. Since the discreteness of time and space in the model is not independent, but is linked by the discrete velocity of the particle, the movement of the particle is divided into two parts: migration and collision. That is, between two time steps, the particle moves from a grid node to the corresponding adjacent grid node and collides with other particles at the grid node. This forms the evolution equation of the distribution function during the migration and collision process, namely:

[0136]

[0137] In the formula, the distribution function at point x at time t and the equilibrium distribution function are solved to obtain the adjacent points x: The distribution function at position t+Δt time; [(M a ) -1 Λ a M a ] represents the relaxation time, that is, the time required for the fluid to transition from the current state to the equilibrium state. This method is based on Λ a Parameter settings can ensure the stable migration of oxygen and water with different viscosities under different working conditions; represents the force terms of different phase fluids, represents the source term of the electrochemical reaction of different phase fluids, where M a Represented as a matrix of specific values:

[0138]

[0139] Λ a is a diagonal relaxation matrix, expressed as:

[0140]

[0141] Furthermore, in one embodiment of the present invention, in this study, the parameter The value is Relaxation time and fluid viscosity υ a Has the following relationship:

[0142]

[0143] The present invention is through a The diagonal relaxation matrix parameters are set to ensure the stable migration of oxygen and water in the fuel cell under different operating conditions and different viscosity fluid parameters. is the force term related to the interfacial tension between the fluid phase interface between oxygen and water and the solid wall, expressed as:

[0144]

[0145] in A potential function that depends on the local density and interaction strength is used to represent the interaction force between oxygen and water and the interaction force with the solid wall. It consists of two parts. The first is the interaction force between oxygen and water. The second is the force between oxygen, water and the wall The formula is as follows:

[0146]

[0147] The evolution equation in step 3 (3) represents the interaction force between oxygen and water and the interaction force between oxygen and the wall: By adding an appropriate potential function to the force equation, the fluid will automatically be separated into different phases. The potential function is also called the interparticle potential energy ψ a (x, t). The interaction force between water and oxygen is The solution process is as follows:

[0148]

[0149] In the subscripts, a and β represent different fluids; for example, if a represents water, β represents oxygen; if a represents oxygen, β represents water;

[0150] The potential function ψ a The formula for (x, t) is:

[0151] ψ a (x,t)=1-exp(-ρ a (x,t))

[0152] The specific value of the interaction force parameter between water and oxygen is determined by multiple tests based on the changes in the surface tension parameters of water and oxygen under different working conditions and the changes in the contact angles when different phase fluids come into contact. 氧-水 =1.40;

[0153] The interaction force between oxygen and water and the wall in step 3 (3) The solution process is as follows:

[0154]

[0155] In the formula represents the interaction strength between the fluid phase a and the solid phase s, when When , it indicates the non-wetting phase; When , it is represented as the wetting phase. s(x) is the indicator function, s(x) = 1 represents solid, s(x) = 0 represents fluid. Based on the wettability relationship between water, oxygen and solid wall in the fuel cell diffusion layer, Set to:

[0156] In the evolution equation of step 3 (3) The term represents the source term of the electrochemical reaction of different phase fluids, where The calculation process is as follows:

[0157]

[0158] where R a (x, t) represents the electrochemical reaction rate, including the water generation rate R 水 (x, t) and the oxygen consumption rate R 氧(x, t), solved by (3) in step 2.

[0159] Furthermore, the specific steps of step 4 are as follows:

[0160] The new discrete distribution function after the iterative evolution of water and oxygen at each grid point Substitute into the following formula, combined with the spatial distribution of discrete quantities The macroscopic density of water and oxygen in the fuel cell diffusion layer is obtained as follows: a (x,t) and macroscopic velocity

[0161]

[0162] The velocity vector in the above formula Expressed as:

[0163]

[0164] Sat 氧 (x, t) represents oxygen saturation, Sat 水 (x, t) represents water saturation, and the specific calculation formula is as follows:

[0165]

[0166] where m 氧 =32g / mol represents the molar mass of oxygen, m 水 =18g / mol represents the molar mass of water.

[0167] The macroscopic pressure P(x,t) can be obtained by the following formula ψ 氧 :

[0168]

[0169] where c s , G 氧-水 , ψ 氧 (x,t) and ψ 水 The value of (x,t) is determined by step 3.

[0170] Furthermore, the specific steps of step 5 are as follows:

[0171] Since the diffusion rate of water in the fuel cell is much smaller than that of oxygen, and the concentration factor affecting the electrochemical reaction rate in step 2 is mainly reflected in the change of the local oxygen concentration Coxygen(x,t) near the catalyst layer, the mass transfer process is only solved in oxygen, where the control equation affecting the mass transfer process is expressed as:

[0172]

[0173] Where R oxygen (x, t) is the source term of the mass transfer process control equation obtained by the basic physical parameters of step 2 and the oxygen consumption rate in the electrochemical reaction, and then the macroscopic velocity calculated in step 4 is obtained. Substitute into the mass transfer process control equation, D 氧 is the diffusion coefficient of oxygen, and the new oxygen concentration Coxygen(x,t) at each lattice point is obtained by iterative calculation of all lattice points through the migration collision process.

[0174] The boundary conditions for the mass transfer process include: a constant oxygen concentration at the top inlet of the computational domain; the electrochemical reaction in the catalytic layer in the lower part of the computational domain, which involves dynamic changes in the oxygen concentration; periodic boundary conditions on the left and right sides of the computational domain; and a no-slip mirror-bounce boundary for the solid walls within the computational domain.

[0175] Furthermore, the specific steps of step 6 are as follows:

[0176] Because the electrochemical reaction in the fuel cell will have a great influence on the temperature field due to heat release, an internal heat source is defined in the reaction area of ​​the catalyst layer in the simulation method to simulate the heat release problem of the reaction, and the heat generated by the catalyst layer is coupled to the heat conduction and heat convection process gradually transferred to the diffusion layer, affecting the entire electrode temperature evolution and the electrochemical reaction efficiency of the fuel cell. The source term added to the heat transfer control equation is:

[0177] S T (x, t) = ΔH·R 水 (x,t) / (ρ all (x,t)c p (x,t))

[0178] Where ΔH represents the heat of reaction, R 水 (x, t) is the electrochemical reaction rate of water calculated in step 2, ρ all (x, t) represents the density of different grid points in the computational domain. Due to different expressions for the flow region and the solid region, the density in the flow region can be expressed as:

[0179] ρ all (x,t)=ρ 水 (x,t)+ρ 氧 (x,t)

[0180] The parameter ρ in the above formula 水 (x,t) and ρ 氧 (x, t) can be calculated from step 4 and can be expressed in the solid region as:

[0181] ρ all (x,t)=ρ 固体 (x,t)

[0182] ρ 固体 (x, t) represents the density of solid particles, which needs to be set as a constant value. p (x, t) represents the constant-pressure specific heat capacity of multi-component fluids and solids at different lattice points, which can be expressed in the flow region as:

[0183]

[0184] In the solid region it can be expressed as:

[0185]

[0186] The above formula parameters and The value of n can be obtained in step 2, v (x, t) represents the proportion of water and oxygen at different lattice points. The specific formula is as follows:

[0187]

[0188] The oxygen saturation Sat in the above formula 氧 (x, t) and water saturation Sat 水 (x, t) can be obtained in step 4. Then, the macroscopic velocity calculated in step 4 is brought into the fluid-solid coupling heat transfer process control equation for specific solution. The specific formula is expressed as follows:

[0189]

[0190] where λ(x,t) represents the thermal conductivity of the multi-component fluid and solid at different lattice points, which can be expressed in the flow region as:

[0191] λ(x,t)=n v (x,t)*λ 氧 +(1.0-n v (x,t))*λ 水

[0192] In the solid region it can be expressed as:

[0193] λ(x,t)=λ 固体

[0194] The parameter λ in the above formula 水 ,λ 氧 and λ 固体 The value of can be obtained in step 2. By traversing all the grid points through the migration collision process, the evolution calculation of the fluid-solid coupling heat transfer process control equation is performed to obtain the new macroscopic temperature T(x,t) at each grid point.

[0195] The boundary conditions for the heat transfer process include: a constant temperature distribution at the top inlet of the computational domain, an internal heat source in the lower part of the computational domain that changes with the electrochemical reaction, and periodic boundary conditions on the left and right sides of the computational domain, ultimately achieving dynamic temperature changes in the entire computational domain.

[0196] Furthermore, the specific steps of step 7 are as follows:

[0197] When the process reaches the nth iteration and traverses all grid points, the equation is satisfied:

[0198] |Sat 水 (x,n)-Sat 水 (x,n-1)|<ε1

[0199]

[0200] Where ε1 is the steady-state threshold for water saturation, and ε2 is the steady-state threshold for macroscopic velocity. At this point, the overall electrochemical reaction in the fuel cell has reached a steady state, and the iteration cycle ends. The output includes: water and oxygen saturation values ​​for all grid points within the computational domain at different time steps, velocity gradient values ​​within the computational domain in different spatial directions, temperature distribution within the computational domain in different spatial directions, and pressure distribution within the computational domain in different spatial directions. By statistically analyzing these results, we can generate water saturation curves for the fuel cell's diffusion layer and flow layer, as well as images of the oxygen and water transport distribution within the electrode at steady state, velocity gradient images, and temperature distribution images.

[0201] Example 1

[0202] like Figure 1 The figure shows a flow chart of the method for fine-scale simulation of the diffusion layer of a fuel cell. The specific steps are as follows:

[0203] Step 1: Using the structural characteristics of the diffusion layer of the real electrode material, combined with the simplified flow layer and catalyst layer, an RGB vector image of the fuel cell structure is formed. This image is then binarized into a binary matrix of corresponding pixels and an iterative step size is set to convert the pixels into corresponding grid points of equal proportions. This serves as the basic computational domain, mainly including:

[0204] (1) Collect the porous medium structure slices of the actual electrode material diffusion layer, use image representation and processing methods to scan the electrode slices, and generate the flow layer and catalytic layer that match the size of the diffusion layer to form an RGB vector image of the overall electrode structure;

[0205] (2) Distinguish the structural features of different levels of electrodes and use image processing methods such as Matlab and OpenCV to convert the RGB vector image of the electrode structure into a binary image, such as Figure 2As shown, black represents particles, white represents the catalytic sites of the reaction process, and gray represents the flow channel. The specific parameters of the calculation domain are shown in Table 1 below:

[0206] Table 1 Main parameters of the computational domain

[0207]

[0208] (3) Figure 2 The binary image shown is converted into a corresponding binary matrix array of 0 and 1 according to the arrangement of the pixel points, and the iterative step size is set to convert the pixel points into a matrix array with the corresponding grid point distribution in equal proportion. This serves as the basic calculation domain of this simulation method, where each grid point serves as a basic calculation unit. The total number of grid points in the calculation domain of this example is: 1000×252.

[0209] Step 2: In the basic calculation domain, based on the initial distribution of oxygen and water in the calculation domain, macroscopic parameters (macroscopic density and macroscopic viscosity of water and oxygen, oxygen concentration, macroscopic velocity, macroscopic temperature, macroscopic pressure), and basic physical property parameters (interaction strength coefficient between oxygen and water, interaction strength coefficient between different phase fluids and wall, Faraday constant, exchange current density, universal gas constant, thermal conductivity, constant pressure specific heat capacity, and overvoltage), calculate the electrochemical reaction rate of oxygen and water under the target working conditions. This mainly includes:

[0210] (1) Determine the initial distribution state of oxygen and water in the basic calculation domain of the electrode. Specifically, the electrode flow layer, diffusion layer, and catalytic layer are initially set to an oxygen saturation state, and a constant oxygen concentration and temperature are injected at the top inlet of the flow layer. The specific settings are:

[0211] C oxygen = 1.0

[0212] T=343K

[0213] (2) Set the macroscopic parameters and basic physical parameters selected in the calculation domain, including the density of water and oxygen ρ 水 , ρ 氧 and viscosity υ 水 、υ 氧 , oxygen concentration C 氧 , temperature T; and basic physical parameters: interaction strength coefficient G between oxygen and water 氧-水 , the interaction strength coefficient between oxygen, water and the wall Faraday constant F, exchange current density i0, universal gas constant R, transfer coefficient α, reference oxygen concentration C 氧,ref , thermal conductivity λ 水 ,λ 氧 ,λ 固体 , specific heat capacity at constant pressure The overpotential η and diffusion coefficient are shown in Table 2:

[0214] Table 2 Basic physical properties

[0215]

[0216]

[0217] (3) Calculate the oxygen consumption rate R in the electrochemical reaction based on the set macroscopic parameters and basic physical parameters 氧 (x, t) and the water generation rate R 水 (x,t). R 氧 (x, t) can be calculated based on the local oxygen concentration C at the catalytic site. 氧 and the reaction rate constant k elec (x, t) is determined, and the specific calculation formula is as follows:

[0218]

[0219] Reaction rate constant k elec (x,t) can be calculated using the following equation:

[0220]

[0221] Where i0 is the exchange current density, F is the Faraday constant, C 氧,ref is the reference oxygen concentration, α is the transfer coefficient, R is the universal gas constant, T(x,t) is the macroscopic temperature distribution, and η is the overpotential. Based on the above parameter values, the oxygen reaction consumption rate R can be determined. 氧 (x, t). The reaction rate of water production R 水 (x, t) can be calculated based on the oxygen consumption rate R 氧 (x, t) is obtained, and the specific calculation formula is as follows:

[0222] R 水 (x,t)=-2.0*R 氧 (x,t)

[0223] Step 3: In the basic computational domain, based on the initial distribution, macroscopic parameters, basic physical properties, and electrochemical reaction rates of oxygen and water determined in Step 2, determine the initial discrete velocity distribution functions of oxygen and water under the target operating conditions. By performing migration and collision on the evolution equation of the discrete velocity distribution function, the calculation traverses all grid points to obtain a new discrete velocity distribution function for each grid point. This mainly includes:

[0224] (1) According to the initial distribution, macroscopic parameters, basic physical parameters and electrochemical reaction rate of oxygen and water determined in step 2, the initial discrete distribution function of oxygen and water under the target working condition is determined, where the initial discrete distribution function of oxygen and water is Substitute the macroscopic densities of water and oxygen ρwater, ρoxygen in step (2) into the equilibrium distribution function, where the initial velocity If set to 0, and It can be calculated by the following formula:

[0225]

[0226] Then, in the iterative process at different times, the newly calculated macroscopic density of water and oxygen ρ 水 ,ρ 氧 , macroscopic velocity is brought into the discrete equilibrium distribution function Solve in the calculation formula:

[0227]

[0228] Where k is set from 0 to 8, including 9 discrete quantities in total, a represents the serial number of oxygen and water in the calculation domain; w k is the weight factor of different discrete quantities, and its specific value is:

[0229]

[0230] Discrete equilibrium distribution function c represents the grid velocity, in this example c = 1; Represents the spatial distribution of discrete quantities, such as Figure 3 As shown, in this example, since there are 9 discrete quantities, It can be expressed as:

[0231]

[0232] (2) In the boundary condition setting of the computational domain, the left and right boundaries are set as periodic boundaries, so the specific expression is:

[0233]

[0234] The upper boundary is set to the same constant oxygen density as the initial value, and the lower boundary contains the catalytic reaction sites, involving electrochemical reactions. The solid boundary setting within the computational domain uses mirror reflection of the no-slip wall. That is, if a cell at a certain location is determined to be a boundary entity, the normal collision process is omitted, and the density is rebounded. Specifically, it can be expressed as:

[0235]

[0236] (3) By calculating the evolution equation during the migration and collision process, all grid points are traversed to obtain a new discrete distribution function for each grid point. Since the discreteness of time and space in the model is not independent, but is linked by the discrete velocity of the particle, the movement of the particle is divided into two parts: migration and collision. That is, between two time steps, the particle moves from a grid node to the corresponding adjacent grid node and collides with other particles at the grid node. This forms the evolution equation of the distribution function during the migration and collision process, namely:

[0237]

[0238] In the formula, the distribution function at point x at time t and the equilibrium distribution function are solved to obtain the adjacent points x: The distribution function at position t+Δt time; [(M a ) -1 Λ a M a ] represents the relaxation time. This method is based on Λ a The parameter setting satisfies the stable migration of oxygen and water under different working conditions when dealing with fluid parameters of different densities and viscosities; represents the force terms of different phase fluids, represents the source term of the electrochemical reaction of different phase fluids, where M a Represented as a matrix of specific values:

[0239]

[0240] Λ a is a diagonal relaxation matrix, expressed as:

[0241]

[0242] Furthermore, in one embodiment of the present invention, in this study, the parameter The value is Relaxation time and fluid viscosity υ a Has the following relationship:

[0243]

[0244] Through the a The diagonal relaxation matrix parameters are set to meet the requirements of stable migration of oxygen and water in fuel cells under different working conditions and different viscosity fluid parameters. is the force term related to the interfacial tension between the fluid phase interface between oxygen and water and the solid wall, expressed as:

[0245]

[0246] in A potential function that depends on the local density and interaction strength is used to represent the interaction force between oxygen and water and the interaction force with the solid wall. It consists of two parts. The first is the interaction force between oxygen and water. The second is the force between oxygen, water and the wall The formula is as follows:

[0247]

[0248] The evolution equation in step 3 (3) represents the interaction force between oxygen and water and the interaction force between oxygen and the wall: By adding an appropriate potential function to the force equation, the fluid will automatically be separated into different phases. The potential function is also called the interparticle potential energy ψ a (x, t). The interaction force between water and oxygen is The solution process is as follows:

[0249]

[0250] In the subscripts, a and β represent different fluids; for example, if a represents water, β represents oxygen; if a represents oxygen, β represents water;

[0251] The potential function ψ a The formula for (x, t) is:

[0252] ψ a (x,t)=1-exp(-ρ a (x,t))

[0253] The specific values ​​of the interaction force parameters between water and oxygen are determined by multiple tests based on the changes in the surface tension parameters of water and oxygen under different working conditions and the changes in the contact angles of different phase fluids. 氧-水 =1.40;

[0254] The interaction force between oxygen and water and the wall in step 3 (3) The solution process is as follows:

[0255]

[0256] In the formula represents the interaction strength between the fluid phase a and the solid phase s, when When , it indicates the non-wetting phase; When , it is represented as the wetting phase. s(x) is the indicator function, s(x) = 1 represents solid, s(x) = 0 represents fluid. Based on the wettability relationship between water and oxygen and the solid wall in the fuel cell diffusion layer, Set to:

[0257] In the evolution equation of step 3 (3) The term represents the source term of the electrochemical reaction of different phase fluids, where The calculation process is as follows:

[0258]

[0259] where R a (x, t) represents the electrochemical reaction rate, including the water generation rate r 水 (x, t) and the oxygen consumption rate R 氧 (x, t), solved by (3) in step 2.

[0260] Step 4: Calculate the water and oxygen saturation, macroscopic pressure, and macroscopic velocity at each grid point using the new discrete velocity distribution function calculated in step 3, which mainly includes:

[0261] The new discrete distribution function after the iterative evolution of water and oxygen at each grid point Substitute into the following formula, combined with the spatial distribution of discrete quantities The macroscopic density of water and oxygen in the fuel cell diffusion layer is obtained as follows: a (x,t) and macroscopic velocity

[0262]

[0263] The velocity vector in the above formula Expressed as:

[0264]

[0265] Sat 氧 (x, t) represents oxygen saturation, Sat 水 (x, t) represents water saturation, and the specific calculation formula is as follows:

[0266]

[0267] where m 氧 =32g / mol represents the molar mass of oxygen, m 水 =18g / mol represents the molar mass of water.

[0268] The macroscopic pressure P(x,t) can be obtained by the following formula:

[0269]

[0270] Step 5: Substitute the basic physical parameters and the calculated electrochemical reaction rate of water from step 2 into the source term of the mass transfer process control equation. Substitute the macroscopic velocity calculated from step 4 into the mass transfer process control equation, and traverse all grid points through the migration collision process to obtain the new oxygen concentration at each grid point. This mainly includes:

[0271] Since the diffusion rate of water in a fuel cell is much smaller than that of oxygen, the main calculation is the mass transfer and diffusion process of oxygen in the computational domain. The governing equation of oxygen in the mass transfer process is expressed as:

[0272]

[0273] The electrochemical reaction process consumes oxygen at a rate R 氧 (x, t) and the oxygen diffusion rate D 氧 It can be calculated by step 2, and then the macroscopic velocity calculated by step 4 Substitute the mass transfer process control equation and iterate through all the grid points through the migration collision process to obtain the new oxygen concentration C at each grid point. 氧 (x,t).

[0274] The boundary conditions for the mass transfer process include: a constant oxygen concentration at the top inlet of the computational domain; the electrochemical reaction in the catalytic layer in the lower part of the computational domain, which involves dynamic changes in oxygen concentration; and periodic boundary conditions on the left and right sides of the computational domain. The specific formula is as follows:

[0275] C 氧 (x max +1,t)=C 氧 (x min ,t)

[0276] C 氧 (x min -1,t)=C 氧 (x max ,t)

[0277] Step 6: Use the basic physical properties and the calculated electrochemical reaction rate of water from step 2 to bring them into the source term of the heat transfer process control equation. Use the macroscopic velocity calculated from step 4 to bring it into the fluid-solid coupling heat transfer process control equation and traverse all grid points through the migration collision process to obtain the new macroscopic temperature of each grid point. This mainly includes:

[0278] Since the electrochemical reaction in the fuel cell generates a certain amount of heat, it is necessary to define an internal heat source in the reaction area of ​​the catalyst layer to simulate the exothermic reaction. The source term added to the heat transfer control equation is:

[0279] S T (x, t) = ΔH·R 水(x,t) / (ρ all (x,t)c p (x,t))

[0280] Where ΔH represents the heat of reaction, which is set to 286.0 kJ / mol in this example. Rwater(x,t) is the electrochemical reaction rate of water calculated in step 2. all (x, t) represents the density of different grid points in the computational domain. Due to different expressions for the flow region and the solid region, the density in the flow region can be expressed as:

[0281] ρ all (x,t)=ρ 水 (x,t)+ρ 氧 (x,t)

[0282] The parameter ρ in the above formula 水 (x,t) and ρ 氧 (x, t) can be calculated from step 4 and can be expressed in the solid region as:

[0283] ρ all (x,t)=ρ 固体 (x,t)

[0284] ρ 固体 (x,t) represents the density of the solid particles, which is set to 2.0 kg / m in this example. 3 . c p (x, t) represents the constant-pressure specific heat capacity of multi-component fluids and solids at different lattice points, which can be expressed in the flow region as:

[0285]

[0286] In the solid region it can be expressed as:

[0287]

[0288] The above formula parameters and The value of n can be obtained in step 2, v (x, t) represents the proportion of water phase and oxygen phase at different lattice points. The specific formula is as follows:

[0289]

[0290] The oxygen saturation Sat in the above formula 氧 (x, t) and water saturation Sat 水 (x, t) can be obtained in step 4. Then, the macroscopic velocity calculated in step 4 is brought into the fluid-solid coupling heat transfer process control equation for specific solution. The specific formula is expressed as follows:

[0291]

[0292] where λ(x,t) represents the thermal conductivity of the multi-component fluid and solid at different lattice points, which can be expressed in the flow region as:

[0293] λ(x,t)=n v (x,t)*λ 氧 +(1.0-n v (x,t))*λ 水

[0294] In the solid region it can be expressed as:

[0295] λ(x,t)=λ 固体

[0296] The parameter λ in the above formula 水 ,λ 氧 and λ 固体 The value of can be obtained in step 2. By traversing all the grid points through the migration collision process, the evolution calculation of the fluid-solid coupling heat transfer process control equation is performed to obtain the new macroscopic temperature T(x,t) at each grid point.

[0297] The boundary conditions for the heat transfer process include: a constant temperature distribution at the top inlet of the computational domain, an internal heat source in the lower part of the computational domain that changes with the electrochemical reaction, and periodic boundary conditions on the left and right sides of the computational domain. The specific formula can be expressed as:

[0298] T(x max +1,t)=T(x min ,t)

[0299] T(x min -1,t)=T(x max ,t)

[0300] Step 7: Substitute the new macroscopic temperature and oxygen concentration calculated in steps 5 and 6 into step 2 to obtain a new electrochemical reaction rate, and repeat steps 2, 3, 4, 5, and 6 to perform a new round of iterative calculations. When the water saturation and velocity gradient values ​​do not exceed the set threshold, it is determined that the overall electrochemical reaction of the fuel cell has reached a stable state, and the iterative cycle ends. It mainly includes:

[0301] When the process reaches the nth iteration and traverses all grid points, the equation is satisfied:

[0302] |Sat 水 (x,n)-Sat 水 (x,n-1)|<1×10 -5

[0303]

[0304] At this point, the overall electrochemical reaction of the fuel cell has reached a stable state, and the loop iteration is ended. The output content at this time includes: water saturation and oxygen saturation of all grid points in the calculation domain at different time steps, velocity gradient in the calculation domain at different spatial directions, temperature distribution of the calculation domain at different spatial directions, and pressure distribution of the calculation domain at different spatial directions. By statistically analyzing these result values, the water saturation change curve of the fuel cell diffusion layer and flow layer, the migration distribution image of oxygen and water in the electrode at steady state, the velocity gradient image, and the temperature distribution image are made ( Figure 4 、 Figure 5 、 Figure 6 、 Figure 7 ). Figure 4-7 It can express the actual reaction transport process and characteristic changes of the diffusion layer during the reaction process at the nanometer to micrometer level. The oxygen and water migration distribution images and water saturation change curves can be used to evaluate the strength of the diffusion layer structure on the oxygen diffusion capacity, the evaluation of the diffusion layer structure on the water transport capacity, and the analysis of which structure is more likely to cause concentration polarization and water flooding anomalies. Further analysis of the velocity gradient image and temperature distribution image is also specifically combined with the diffusion layer structure. By analyzing where the extreme temperature mainly occurs, the GDL structure is optimized to prevent thermal runaway of the battery caused by excessive temperature. Since these images all represent the changes in physical quantities of specific reaction processes at the nanometer to micrometer level, they have more refined characteristics.

[0305] In summary, the present invention discloses a method for fine simulation of the diffusion layer of a fuel cell at mesoscale. This method performs binarization processing and converts it into matrix data to construct a basic computational domain based on the structural characteristics of the diffusion layer of the real electrode material and the scale characteristics of the flow layer and the catalytic layer. Then, a discrete velocity evolution equation, a fluid-solid coupling heat transfer evolution equation, and a mass transfer evolution equation containing the dynamic electrochemical reaction rate are established. After determining the initial and boundary conditions, a fine simulation of the multi-physical field coupling of the fuel cell is started in the computational domain. The real-time changes of the velocity field, concentration field, and temperature field are simulated to realize the study of the overall electrochemical reaction performance of the fuel cell. After reaching a steady state, the migration characteristics and change curves of water and oxygen in the diffusion layer, and the macroscopic velocity and temperature change images are output. This method is oriented towards the heterogeneous electrode structural characteristics of the real diffusion layer material. It is suitable for studying the overall dynamic electrochemical reaction performance of the fuel cell based on the diffusion layer and in combination with other functional levels, as well as the interaction and transmission mechanism of the multi-component flow, heat transfer, and mass transfer processes in the porous electrode.

[0306] The above description of the embodiments is intended to facilitate understanding and use of the invention by those skilled in the art. It will be apparent that those skilled in the art can readily make various modifications to these embodiments and apply the general principles described herein to other embodiments without requiring inventive effort. Therefore, the present invention is not limited to the above-described embodiments. 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 present invention, should be within the scope of protection of the present invention.

Claims

1. A mesoscale simulation method for a fuel cell diffusion layer, characterized in that: The steps include: Step 1: Construct a calculation domain of a diffusion layer structure, wherein the top of the calculation domain is the flow layer and the bottom of the calculation domain is the catalytic layer; Step 2: Calculate the electrochemical reaction rate under the target working conditions within the calculation domain; Step 3: Determine the discrete velocity distribution function of each grid point in the computational domain; Step 4: Calculate the water and oxygen saturation, macroscopic pressure, and macroscopic velocity; Step 5: Determine the mass transfer process control equation for each grid point in the computational domain and calculate the oxygen concentration; Step 6: Determine the governing equations of the fluid-solid coupled heat transfer process at each grid point in the computational domain and calculate the macroscopic temperature; Step 7: Determine whether the water saturation and macro velocity meet the set thresholds: If the conditions are met, it is determined that a stable state has been reached, and the water saturation change curve, the migration distribution image of oxygen and water in the electrode at steady state, the velocity gradient image and the temperature distribution image are output; If not, it is determined that the steady state has not been reached, and the macroscopic temperature and oxygen concentration calculated in steps 5 and 6 are substituted into step 2 for iteration; Step 1 includes the following steps: (1) Collect the porous media structure slices of the diffusion layer of the electrode material in the fuel cell, scan the slices, generate the flow layer and catalyst layer that match the size of the diffusion layer, and form an RGB vector image of the overall electrode structure; (2) Distinguish the structural features of different levels in the electrode structure and convert the RGB vector image into a binary image; (3) Convert the binary image into a binary matrix array according to the arrangement of pixel points as the calculation domain; Step 2 includes the following steps: (1) Determine the initial distribution of oxygen and water in the computational domain; (2) Set the macroscopic parameters and basic physical property parameters of the computational domain; (3) Calculate the oxygen consumption rate and water generation rate in the electrochemical reaction under the target operating conditions: Oxygen consumption rate : ； Reaction rate constant Calculated according to the following equation: ； in, is the exchange current density, is the Faraday constant, is the reference oxygen concentration, is the transfer coefficient, is the universal gas constant, is the macroscopic temperature, Overpotential Water production rate : 。 2. The mesoscale simulation method of a fuel cell diffusion layer according to claim 1, characterized in that: Step 3 includes the following steps: (1) Determine the initial discrete velocity distribution functions of oxygen and water : ； Where k represents the serial number of the discrete quantity, and a represents the serial number of different fluids in the computational domain; are weight factors for different discrete quantities; , c represents the grid velocity; Represents the spatial distribution of discrete quantities; Indicates the macroscopic density of different phase fluids; Indicates the macroscopic velocity of different phase fluids; (2) Constructing the interaction forces between fluids and between fluids and walls : ； in, ； , in the subscript, and Represent different fluids respectively; in, , Indicates the strength of oxygen-water interaction; ； in, represents the interaction strength between the fluid phase a and the solid phase s, is the indicator function, Represents a solid, Represents fluid; ； (3) Establish the evolution equation of the discrete distribution function considering the dynamic electrochemical reaction rate: in, represents the relaxation time, represents the force terms of different phase fluids, represents the source term of the electrochemical reaction of different phase fluids, is the time step; ； ； ； ； Represents the electrochemical reaction rate, including the oxygen consumption rate and water production rate.

3. The mesoscale simulation method of a fuel cell diffusion layer according to claim 2, characterized in that: At least one of ; Among them, relaxation time and fluid viscosity satisfy: 。 4. The mesoscale simulation method of a fuel cell diffusion layer according to claim 1, characterized in that: In step 4: Macro speed : ； in, represents the initial relaxation time, and a represents the sequence number of different fluids in the computational domain; , represents the discrete velocity distribution function, k represents the sequence number of the discrete quantity; , Represents the spatial distribution of discrete quantities; ； Water saturation : ； in, represents the molar mass of water; Oxygen saturation: : ； in, represents the molar mass of oxygen; Macro pressure : ； in, , c represents the grid velocity; Indicates the macroscopic density of oxygen; Indicates the macroscopic density of water; Indicates the strength of oxygen-water interaction; represents the potential function of oxygen; represents the potential function of water.

5. The mesoscale simulation method of a fuel cell diffusion layer according to claim 1, characterized in that: In step 5: The mass transfer process control equation is: ； in, is the oxygen consumption rate in the electrochemical reaction rate obtained in step 2, is the macroscopic velocity obtained in step 4, is the diffusion coefficient of oxygen, is the oxygen concentration.

6. The mesoscale simulation method of a fuel cell diffusion layer according to claim 1, characterized in that: In step 6: The governing equation for the fluid-solid coupled heat transfer process is: ； in, represents macroscopic temperature; , represents the heat of reaction, is the water generation rate in the electrochemical reaction rate obtained in step 2, represents the density of different grid points in the computational domain, It represents the constant-pressure specific heat capacity of multi-component fluids and solids at different lattice points; In the flow area: , represents the macroscopic density of water, Indicates the macroscopic density of oxygen; , represents the specific heat capacity of water at constant pressure, represents the specific heat capacity of oxygen at constant pressure; , Indicates oxygen saturation, Indicates water saturation; , represents the thermal conductivity of water, represents the oxygen thermal conductivity; In solid regions: , Indicates the density value of solid particles; , represents the specific heat capacity of solid at constant pressure; , represents the thermal conductivity of solid.

7. The mesoscale simulation method of a fuel cell diffusion layer according to claim 1, characterized in that: In step 7: judge: ； ； in, Indicates water saturation, represents the macroscopic speed, is the steady-state threshold of water saturation, is the steady-state threshold of macroscopic velocity.

8. The mesoscale simulation method of a fuel cell diffusion layer according to claim 1, characterized in that: In step 7, the output results are obtained through the water saturation and oxygen saturation values ​​of all grid points in the calculation domain, the velocity gradient values ​​in the calculation domain in different spatial directions, the temperature distribution of the calculation domain in different spatial directions, and the pressure distribution of the calculation domain in different spatial directions.

Images