Mesoscale simulation method of fuel cell diffusion layer

Abstract

The invention belongs to the technical field of fuel cells, and particularly relates to a mesoscale simulation method of a fuel cell diffusion layer, which comprises the following steps of: constructing a computational domain of a diffusion layer structure; calculating the electrochemical reaction rate under the target working condition in the computational domain; determining a discrete velocity distribution function of each grid point in the computational domain; calculating saturation, macroscopic pressure and macroscopic speed of water and oxygen; a mass transfer process control equation of each grid point in the computational domain is determined, and the oxygen concentration is calculated; determining a fluid-structure interaction heat transfer process control equation of each grid point in the computational domain and calculating a macroscopic temperature; and judging whether the water saturation and the macroscopic velocity meet threshold values or not. Compared with the prior art, the method solves the problem that in the prior art, a simulation method for accurate coupling and mutual feedback between mesoscale multi-physics fields of a battery diffusion layer is difficult. According to the scheme, a mesoscale fine simulation method for multiphase flow mass transfer, dynamic electrochemical reaction and heat transfer in a flow layer, a diffusion layer and a catalyst layer of the fuel cell is constructed.

Description

Technical Field

[0001] The 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. In particular, fuel cells, as an electrochemical energy conversion device, have the advantages of high efficiency, fast response speed, high power density, and zero pollution to the environment. However, the performance of fuel cells depends largely 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 can maintain high proton conductivity, but water under high load may condense into droplets in the porous electrodes. If not adequately removed, liquid water will eventually accumulate in the diffusion layer and catalyst layer, hindering the diffusion of reactants and causing flooding, ultimately reducing the overall performance of the fuel cell. Some studies have shown that optimizing the electrode diffusion layer can effectively alleviate the flooding phenomenon. The functions of the fuel cell diffusion layer include: removing accumulated water in the fuel cell, providing a pathway for reactants to diffuse to the catalyst layer, providing structural support for the reaction area, and promoting electron and heat transfer. Flooding can be further alleviated by designing different structural features of the diffusion layer or changing the surface properties of the diffusion layer.

[0003] However, the existing simulation methods mainly focus on macroscopic or system research. For example, in the numerical simulation method of a solid oxide fuel cell disclosed in CN116467913A, a macroscopic simulation method of coupled heat transfer, mass transfer, and electrochemical reaction in the stack is established, which realizes the rapid simulation of solid oxide fuel cells; in another example, in the computer simulation method of mass, momentum, energy, and charge transfer in proton exchange membrane fuel cells disclosed in CN118472326A, a physical transmission method between different layers in the electrode is established, which realizes the simulation of various fluid transmission and migration characteristics in the electrode. The above schemes can realize the performance analysis and parameter control of fuel cells, but they are mainly based on constantly adjusting parameters to find rules and then determine a more reasonable electrode design scheme. This method does not have a deep understanding of the interaction between multiphase flow and complex electrode structure, the coupling process and microscopic mechanism between heat transfer, mass transfer and electrochemical reaction, so it is often affected by experience in determining the scheme. Mesoscale research needs to consider the interaction under multiple dynamic conditions such as fluid dynamics, mass transfer, heat transfer, and electrochemical reaction at the same time, which can further explore the mechanism and influence mechanism 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 relationship between the physical behaviors from the microscopic to the macroscopic scale. (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 the actual fuel cell operates 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 scheme, based on the structural characteristics of the fuel cell diffusion layer and combined with the simplified scale characteristics of other functional levels (flow layer and catalyst layer), constructs a mesoscale fine simulation method for multiphase flow mass transfer, dynamic electrochemical reaction and heat transfer inside the flow layer, diffusion layer and catalyst layer of the fuel cell. This method is oriented to 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 levels, as well as the interaction and transmission mechanism of the multi-component flow, heat transfer and mass transfer processes in the porous electrode.

[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 condition in 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, macro pressure and macro velocity;

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

[0014] Step 6: Determine the control equation of the fluid-solid coupling heat transfer process at each grid point in the calculation 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 a 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 in the calculation domain adopts a no-slip rebound boundary.

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

[0020] (1) Collecting porous media structure slices of the diffusion layer of the electrode material in the fuel cell, scanning the slices, generating a flow layer and a catalyst layer that match the size of the diffusion layer, and 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) The binary image is converted into a binary matrix array according to the arrangement position of the pixels as the calculation domain.

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

[0024] (1) Determine the initial distribution state 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] Among them, i 0 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 generation 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 speed; 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 term 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 generation rate.

[0054] Preferably, a At least one of

[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 constant value related to the fluid viscosity υ a Related values ​​(corresponding to different working conditions), where at least one Need to be consistent with fluid viscosity a Related.

[0058] Preferably, in step 4:

[0059] Macro speed

[0060]

[0061] in, represents the initial relaxation time, 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] Macro 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 macro 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 control equation of fluid-solid coupling 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 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, It 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 the solid region:

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

[0089] It 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 macroscopic velocity, ε 1 is the steady-state threshold of water saturation, ε 2 is the steady-state threshold of the macroscopic 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 the real electrode material, combined with the flow layer and catalyst layer of the same scale characteristics to carry out simulation research, which 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 can simulate the generation and migration of water and its actual impact on the dynamic electrochemical reaction of the fuel cell by establishing a set of fuel cell dynamic water management simulation methods, which can help optimize the design of the diffusion layer, alleviate the flooding phenomenon, and improve the overall performance of the fuel cell. BRIEF DESCRIPTION OF THE DRAWINGS

[0101] Figure 1 It is a schematic flow chart of the method for fine simulation of the fuel cell diffusion layer mesoscale 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 is 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 the steady state in Example 1;

[0106] Figure 6 The velocity gradient image in the electrode of the fuel cell in the 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 a steady state. DETAILED DESCRIPTION

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

[0109] In order to overcome the shortcomings of the prior art, the present invention constructs a mesoscale fine simulation method for multiphase flow mass transfer, dynamic electrochemical reaction and heat transfer in the flow layer, diffusion layer and catalyst layer of the fuel cell based on the structural characteristics of the fuel cell diffusion layer and combined with the simplified scale characteristics of other functional layers (flow layer and catalyst layer). This method is oriented to 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, 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, and specifically 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, the RGB vector image of the fuel cell structure is formed, which is binarized into a binary matrix of corresponding pixels and the iteration step is set to convert the pixels into corresponding grid points in the same proportion 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 according to 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 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, diffusion coefficient, constant pressure specific heat capacity and overvoltage);

[0112] Step 3: In the basic calculation domain, according to the initial distribution, macroscopic parameters, basic physical parameters and electrochemical reaction rate of oxygen and water determined in step 2, determine the initial discrete velocity distribution function of oxygen and water under the target working condition, and calculate and traverse all grid points by performing migration and collision on the evolution equation of the discrete velocity distribution function to obtain a new discrete velocity distribution function for each grid point;

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

[0114] Step 5: Substitute the basic physical property parameters in step 2 and the calculated electrochemical reaction rate of water into the source term of the mass transfer process control equation, and substitute the macroscopic velocity calculated in 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;

[0115] Step 6: The basic physical property parameters and the calculated electrochemical reaction rate of water in step 2 are brought into the source term of the heat transfer process control equation. The macroscopic velocity calculated in step 4 is brought into the fluid-solid coupling heat transfer process control equation and all grid points are traversed 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 water saturation and velocity gradient numerical changes 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 cycle iteration 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 size of the diffusion layer, 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 iteration step is set to convert the pixel points into a matrix array of corresponding grid point distribution in equal proportion, which serves as the basic calculation domain of this simulation method, in which each grid point is used 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 fluid layer, diffusion layer and catalyst layer are initially set to an oxygen saturation state. A constant oxygen concentration is injected through the fluid layer at the top inlet of the fluid layer, and diffuses into the diffusion layer, and finally reaches the catalyst layer for oxidation-reduction reaction, consuming oxygen and generating water. The temperature at the top inlet of the fluid layer is also kept constant.

[0123] (2) Setting the macroscopic parameters and basic physical property parameters selected in the calculation domain, the macroscopic parameters include: 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 and water and the wall Faraday constant F, exchange current density i 0 , 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 property 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 i 0 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 generation 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 this 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 performing irregular thermal motion all the time, the speed of microscopic particles is continuous, and its speed is infinite in the phase space. However, the details of the particle movement do not significantly determine the macroscopic movement of the fluid. Therefore, the particle speed is simplified to a finite-dimensional speed space in the phase space. Then, the continuous distribution function is also discretized into a discrete distribution function with the same number of discrete speeds. Then the equilibrium distribution function of each grid point will also be discretized into discrete equilibrium distribution functions with the same number of discrete velocities 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] In the formula, k represents the serial number of the discrete quantity (for example, if it is discretized into 9 quantities, k represents 0 to 8), a represents the serial number of different fluids in the calculation 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 grid speed, and its value needs to be converted to the actual sound speed. When the actual sound speed 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 they 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) Boundary conditions for the reaction inside the fuel cell diffusion layer combined with the flow layer and the catalyst layer are set according to different working conditions; the top of the calculation domain is set as an inlet, and a constant oxygen concentration is injected through the flow layer, and diffuses in the diffusion layer and finally reaches the catalyst layer to undergo a redox reaction. The bottom of the calculation domain is the catalyst 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. Therefore, the structure of the catalyst layer can be simplified, and the distribution area of ​​the catalytic reaction sites can be set based on the size of the actual cracks in the catalyst layer. Cracks in the catalyst 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 within the calculation domain is subjected to boundary processing. The present invention adopts a no-slip rebound boundary, that is, a 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-mentioned different boundary conditions, the values ​​of the different discrete velocity distribution functions of the particles need to be modified accordingly;

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

[0136]

[0137] In the formula, the distribution function at point x at time t and the equilibrium distribution function are solved to obtain the distribution function of 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 The parameter setting can meet the stable migration of oxygen and water with different viscosity values ​​under different working conditions; represents the force term 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] Further, 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 setting of the diagonal relaxation matrix parameters satisfies the stable migration of oxygen and water in the fuel cell under different operating conditions when dealing with 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 interaction between oxygen, water and the wall The formula is as follows:

[0146]

[0147] The terms representing the interaction force between oxygen and water and the interaction force between oxygen and the wall in the evolution equation of (3) in step 3 are: By adding an appropriate potential function to the force equation, the fluid will automatically separate into different phases. The potential function is also called the interparticle potential energy ψ a (x, t). The interaction term 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 calculation 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 of different phase fluids in 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 , 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 iteration of water and oxygen evolution 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 It is 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 in step 2 and the oxygen consumption rate in the electrochemical reaction calculated, and then the macroscopic velocity calculated in step 4 is Substitute into the mass transfer process control equation, D 氧 is the diffusion coefficient of oxygen. 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 of the mass transfer process include: a constant oxygen concentration is set at the top inlet of the computational domain, the electrochemical reaction of the catalytic layer in the lower part of the computational domain involves the dynamic change of oxygen concentration, periodic boundary conditions are maintained on the left and right sides of the computational domain, and a no-slip mirror rebound boundary is used for the solid wall in 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 impact 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 convection process gradually to the diffusion layer, which affects the evolution of the entire electrode temperature 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, it can be expressed in the flow region as follows:

[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 fluid and solid at different grid points, which can be expressed as:

[0183]

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

[0185]

[0186] The above formula parameters and The value of can be obtained in step 2, n v (x, t) represents the proportion of water and oxygen at different grid 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 control equation of the fluid-solid coupling heat transfer process for specific solution. The specific formula is expressed as follows:

[0189]

[0190] Where λ(x,t) represents the thermal conductivity of multi-component fluid and solid at different grid points, which can be expressed 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 above formula parameter λ 水 ,λ 氧 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 control equation of the fluid-solid coupling heat transfer process is performed to obtain the new macroscopic temperature T(x, t) at each grid point.

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

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

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

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

[0199]

[0200] In the formula, ε 1 is the steady-state threshold of water saturation, ε 2 is the steady-state threshold of the macroscopic velocity; at this point, the overall electrochemical reaction of the fuel cell has reached a steady state, and the loop iteration ends. At this point, the output includes: the water saturation and oxygen saturation values ​​of all grid points in the computational domain at different time steps, the velocity gradient values ​​in the computational domain at different spatial directions, the temperature distribution of the computational domain at different spatial directions, and the pressure distribution of the computational domain at different spatial directions. By statistically analyzing these result values, the water saturation change curves of the fuel cell diffusion layer and flow layer, the migration distribution images of oxygen and water in the electrode at steady state, the velocity gradient images, and the temperature distribution images are made.

[0201] Example 1

[0202] like Figure 1 The flowchart of the fuel cell diffusion layer mesoscale fine simulation method is shown, and 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, the RGB vector image of the fuel cell structure is formed, which is binarized into a binary matrix of corresponding pixels and the iteration step is set to convert the pixels into corresponding grid points in the same proportion as the basic calculation domain, which mainly includes:

[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 catalyst 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 2 As 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 corresponding 0, 1 binary matrix arrays according to the arrangement position of the pixel points, and the iteration step is set to convert the pixel points into the matrix array of the corresponding grid point distribution in equal proportion, which serves as the basic calculation domain of this simulation method. Each grid point is used 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, according to 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 intensity coefficient between oxygen and water, interaction intensity 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 condition, mainly including:

[0210] (1) Determine the initial distribution state of oxygen and water in the basic calculation domain of the electrode, which is specifically expressed as follows: the electrode flow layer, diffusion layer and catalyst 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, which are specifically set as follows:

[0211] C oxygen = 1.0

[0212] T=343K

[0213] (2) Setting the macroscopic parameters and basic physical property parameters selected in the calculation domain, the macroscopic parameters include: 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 and water and the wall Faraday constant F, exchange current density i 0 , 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 property 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 i 0 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 generation 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 calculation domain, according to the initial distribution, macroscopic parameters, basic physical parameters and electrochemical reaction rate of oxygen and water determined in step 2, determine the initial discrete velocity distribution function of oxygen and water under the target working condition, and calculate and traverse all grid points by migration and collision of the evolution equation of the discrete velocity distribution function to obtain a new discrete velocity distribution function for each grid point, which 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 and ρ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] In the formula, 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 speed, 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 one, and the lower boundary contains the catalytic reaction site, involving electrochemical reactions; and the solid boundary setting inside the computational domain uses the mirror reflection of the non-slip wall. That is, if a cell at a certain position is determined as a boundary entity, the normal collision process is omitted, and the density will be rebounded, which can be specifically expressed as:

[0235]

[0236] (3) Through the calculation of the evolution equation in the migration and collision process, all grid points are traversed to obtain the new discrete distribution function of 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, the particle moves from a grid node to the corresponding adjacent grid node between two time steps and collides with other particles at the grid node. This forms the evolution equation of the distribution function in 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 distribution function of 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 term 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] Further, 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 setting of the diagonal relaxation matrix parameters satisfies 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:

[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 interaction between oxygen, water and the wall The formula is as follows:

[0247]

[0248] The terms representing the interaction force between oxygen and water and the interaction force between oxygen and the wall in the evolution equation of (3) in step 3 are: By adding an appropriate potential function to the force equation, the fluid will automatically separate into different phases. The potential function is also called the interparticle potential energy ψ a (x, t). The interaction term 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 calculation 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 , 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 in the fuel cell diffusion layer and the solid wall, 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 saturation, macro pressure and macro velocity of water and oxygen at this time through the new discrete velocity distribution function of each grid point calculated in step 3, mainly including:

[0261] The new discrete distribution function after the iteration of water and oxygen evolution 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 It is 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: The basic physical parameters in step 2 and the calculated electrochemical reaction rate of water are introduced into the source term of the mass transfer process control equation. The macroscopic velocity calculated in step 4 is introduced into the mass transfer process control equation, and all grid points are traversed through the migration collision process to obtain the new oxygen concentration at each grid point, which mainly includes:

[0271] Since the diffusion rate of water in the fuel cell is much smaller than that of oxygen, the mass transfer diffusion process of oxygen in the computational domain is mainly calculated, where the control 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 diffusion rate of oxygen D 氧 It can be calculated by step 2, and then the macro velocity calculated by step 4 Substitute the mass transfer process control equation into it, 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 of the mass transfer process include: a constant oxygen concentration is set at the top inlet of the computational domain, the electrochemical reaction of the catalytic layer in the lower part of the computational domain involves the dynamic change of oxygen concentration, and periodic boundary conditions are maintained on the left and right sides of the computational domain. The specific formula is expressed 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: The basic physical parameters in step 2 and the electrochemical reaction rate of water calculated are brought into the source term of the heat transfer process control equation. The macroscopic velocity calculated in step 4 is brought into the fluid-solid coupling heat transfer process control equation and all grid points are traversed through the migration collision process to obtain the new macroscopic temperature of each grid point, which mainly includes:

[0278] Since the electrochemical reaction of 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 heat release problem. 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, it can be expressed in the flow region as follows:

[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 fluid and solid at different grid points, which can be expressed as:

[0285]

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

[0287]

[0288] The above formula parameters and The value of can be obtained in step 2, n 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 control equation of the fluid-solid coupling heat transfer process for specific solution. The specific formula is expressed as follows:

[0291]

[0292] Where λ(x,t) represents the thermal conductivity of multi-component fluid and solid at different grid points, which can be expressed 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 above formula parameter λ 水 ,λ 氧 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 control equation of the fluid-solid coupling heat transfer process is performed to obtain the new macroscopic temperature T(x, t) at each grid point.

[0297] The boundary conditions of the heat transfer process include: a constant temperature distribution is set at the top inlet of the computational domain, the internal heat source in the lower part of the computational domain changes with the electrochemical reaction, and the periodic boundary conditions are maintained 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: Bring 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 cycle iteration ends. It mainly includes:

[0301] When the process reaches the nth iteration, all grid points are traversed and 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 ends. 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 at the nanometer to micrometer level during the reaction process. 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 diffusion layer structure on the water transport capacity, and analyze under which structure concentration polarization and flooding anomalies are more likely to occur. 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 are more refined.

[0305] In summary, the present invention discloses a method for fine simulation of the diffusion layer of a fuel cell at a mesoscale. This method performs binarization processing and converts it into matrix data to construct a basic calculation 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 catalyst 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 calculation domain. The overall electrochemical reaction performance of the fuel cell is studied by simulating the real-time changes of the velocity field, concentration field, and temperature field. 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 to 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 levels, as well as the interaction and transmission mechanism of the multi-component flow, heat transfer, and mass transfer processes contained in the porous electrode.

[0306] The above description of the embodiments is to facilitate the understanding and use of the invention by those skilled in the art. It is obvious that those skilled in the art can easily make various modifications to these embodiments and apply the general principles described herein to other embodiments without creative work. Therefore, the present invention is not limited to the above embodiments, and 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 condition in 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, macro pressure and macro velocity; Step 5: Determine the mass transfer process control equation for each grid point in the calculation domain and calculate the oxygen concentration; Step 6: Determine the control equation of the fluid-solid coupling heat transfer process at each grid point in the calculation 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 a 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.

2. The mesoscale simulation method of a fuel cell diffusion layer according to claim 1, characterized in that: Step 1 includes the following steps: (1) Collecting porous media structure slices of the diffusion layer of the electrode material in the fuel cell, scanning the slices, generating a flow layer and a catalyst layer that match the size of the diffusion layer, and forming 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) The binary image is converted into a binary matrix array according to the arrangement position of the pixels as the calculation domain.

3. The mesoscale simulation method of a fuel cell diffusion layer according to claim 1, characterized in that: Step 2 includes the following steps: (1) Determine the initial distribution state of oxygen and water in the computational domain; (2) Setting 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 R 氧 (x,t): Reaction rate constant k elec (x,t) is calculated according to the following equation: 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 Water generation rate R 水 (x,t): R 水 (x,t)=-2.0*R 氧 (x,t)。 4. 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, 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 speed; represents the spatial distribution of discrete quantities; ρ a 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 subscripts, a and β represent different fluids; Among them, ψ a (x,t)=1-exp(-ρ a (x,t)), G 氧-水 Indicates the strength of oxygen-water interaction; 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; (3) Establish the evolution equation of the discrete distribution function considering the dynamic electrochemical reaction rate: Among them, [(M a ) -1 Λ a M a ] represents the relaxation time, represents the force term of different phase fluids, represents the source term of the electrochemical reaction of different phase fluids, and Δt is the time step; Λ a is the diagonal relaxation matrix; R a (x, t) represents the electrochemical reaction rate, including the oxygen consumption rate and water generation rate.

5. The mesoscale simulation method of a fuel cell diffusion layer according to claim 4, characterized in that: Λ a At least one of Among them, relaxation time and fluid viscosity υ a satisfy:

6. 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, 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 Sat 水 (x,t): Among them, m 水 represents the molar mass of water; Oxygen saturation: Sat 氧 (x,t): Among them, m 氧 represents the molar mass of oxygen; Macro pressure P(x,t): 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.

7. 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: Among them, r 氧 (x, t) is the oxygen consumption rate in the electrochemical reaction rate obtained in step 2, is the macro velocity obtained in step 4, D 氧 is the diffusion coefficient of oxygen, C 氧 (x, t) is the oxygen concentration.

8. The mesoscale simulation method of a fuel cell diffusion layer according to claim 1, characterized in that: In step 6: The control equation of fluid-solid coupling heat transfer process is: Where T(x,t) represents the macroscopic temperature; 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 multi-component fluid and solid at different lattice points; In the flow area: ρ all (x,t)=ρ 水 (x,t)+ρ 氧 (x,t),ρ 水 (x, t) represents the macroscopic density of water, ρ 氧 (x, t) represents the macroscopic density of oxygen; represents the specific heat capacity of water at constant pressure, It represents the specific heat capacity of oxygen at constant pressure; Sat 氧 (x, t) represents oxygen saturation, Sat 水 (x, t) represents water saturation; λ(x,t)=n v (x,t)*λ 氧 +(1.0-n v (x,t))*λ 水 ,λ 水 represents the thermal conductivity of water, λ 氧 represents the oxygen thermal conductivity; In the solid region: ρ all (x,t)=ρ 固体 (x,t),ρ 固体 (x, t) represents the density value of solid particles; It represents the specific heat capacity of solid at constant pressure; λ(x,t)=λ 固体 ,λ 固体 Represents the thermal conductivity of solid.

9. The mesoscale simulation method of a fuel cell diffusion layer according to claim 1, characterized in that: In step 7: judge: |Sat 水 (x,n)-Sat 水 (x,n-1)|<ε1; 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.

10. 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 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.

Images