Method and device for predicting the operating performance of a proton exchange membrane fuel cell

Abstract

The application relates to the technical field of fuel cells, in particular to a working performance prediction method and device for a proton exchange membrane fuel cell, wherein the method comprises the following steps: dividing a target single-piece fuel cell into multiple network units, each network unit comprising a bipolar plate, a gas channel, a membrane electrode and cooling liquid on the upper and lower sides for heat dissipation; establishing a relational expression among the multiple network units, combining a physical boundary condition and a second preset boundary condition to construct a two-dimensional network model in a plane direction, and establishing a one-dimensional differential model in the membrane electrode of each network unit; coupling the two-dimensional network model and the one-dimensional differential model to establish a two-dimensional+one-dimensional mathematical model of the proton exchange membrane fuel cell with multi-physical field coupling; and solving the distribution of physical quantities such as current density and temperature by using the mathematical model, and then predicting the performance of the target single-piece fuel cell. Therefore, the problems that the simulation model for the fuel cell cannot meet the membrane electrode rapid iterative design requirement and has high cost and the like are solved.

Description

Technical Field

[0001] This invention relates to the field of fuel cell technology, and in particular to a calculation method and apparatus for rapid prediction of the performance of a proton exchange membrane fuel cell. Background Technology

[0002] Proton exchange membrane fuel cell (PEMFC) technology is a green and efficient new energy power generation technology. It boasts high operating efficiency, emits no high-temperature harmful gases, and has broad application prospects. The core component of a fuel cell is the membrane electrode assembly (MEA). During operation, a complex physical process involving electrochemical reactions, flow, mass transfer, and heat transfer coupling occurs within the MEA. Revealing the coupling mechanisms between these physical fields is one of the key technologies urgently needing breakthroughs in high-performance MEA design. Constrained by the geometric dimensions of the MEA thickness, current technologies struggle to measure the distribution characteristics of various physical quantities within the MEA during operation online. Therefore, numerical simulation technology can be used to deeply analyze the multi-field coupling physical mechanisms within the MEA, rapidly evaluate the fuel cell's performance, and accelerate the iterative design process of the fuel cell MEA.

[0003] The operation of a proton exchange membrane fuel cell (PEMFC) is a complex multiphysics coupling problem involving electrochemical reactions, multiphase and multi-component fluid flow, and various forms of heat transfer. In practical applications, large-scale fuel cell systems exhibit significant non-uniformity in internal gas pressure, temperature, humidity, gas concentration, and water distribution, leading to substantial variations in current density. Therefore, mathematical models used in some related technologies that simply assume uniform distribution of physical quantities within the membrane electrode assembly (MEA) structure (e.g., gas diffusion layer, catalyst layer, proton exchange membrane) or throughout the entire fuel cell will result in distorted calculations and hinder in-depth analysis of the causes of low efficiency in localized areas of the fuel cell, thus impeding MEA design improvements.

[0004] To identify the spatial distribution of various physical quantities, three-dimensional models are also used in related technologies to carry out numerical simulation analysis. However, simulating the operating state of large-size fuel cells requires huge computing resources, has a long calculation cycle, and is costly, which cannot meet the needs of rapid iterative design of membrane electrodes.

[0005] Therefore, a simplified simulation model tailored to the characteristics of fuel cells is needed to improve computational efficiency. Summary of the Invention

[0006] This invention provides a method and apparatus for predicting the performance of a proton exchange membrane fuel cell, which solves the problems of existing three-dimensional models for fuel cells having long calculation cycles, not meeting the requirements for rapid iteration, and overly simplified models failing to identify potential problems in the internal structure of the membrane electrode, thus failing to support the design and development of the membrane electrode.

[0007] The first aspect of this invention provides a method for predicting the performance of a proton exchange membrane fuel cell, comprising the following steps: dividing a target monolithic fuel cell into multiple network units, each network unit including a bipolar plate, a gas channel, a membrane electrode assembly (MEA), and coolant on the upper and lower sides for heat dissipation; establishing the relationship between the multiple network units, constructing a two-dimensional network model in a two-dimensional plane direction in combination with a first preset physical boundary condition, and establishing a one-dimensional differential model in the MEA of each network unit; coupling the two-dimensional network model and the one-dimensional differential model to establish a two-dimensional + one-dimensional mathematical model of the proton exchange membrane fuel cell with multi-physics coupling; and using the two-dimensional + one-dimensional mathematical model of the proton exchange membrane fuel cell with multi-physics coupling to solve for the distribution of current density, pressure, flow rate, temperature, humidity, and gas concentration to predict the performance of the target monolithic fuel cell.

[0008] Optionally, dividing the target monolithic fuel cell into multiple network units includes: setting the two-dimensional plane containing the target monolithic fuel cell as the xy plane, with the positive x-axis as the anode gas flow direction and the z-axis as the thickness direction of the target monolithic fuel cell; dividing the target monolithic fuel cell into N network units along the y-axis direction. y N network units, each with a size of Δy, are uniformly divided along the x-axis. x A plurality of network units with a size of Δx are obtained, and the unit numbers are (i,j).

[0009] Optionally, dividing the target monolithic fuel cell into multiple network units further includes: setting a layer of boundary virtual point grid units with the same size as the network units on the outer wall of the divided target monolithic fuel cell to set the equations corresponding to the boundary conditions.

[0010] Optionally, the two-dimensional network model is:

[0011]

[0012] in, denoted as r, the flow rates of component r in the mixed gas at the center of network cell (i,j) and the center of network cell (i,j+1), respectively. r represents hydrogen and water at the anode and oxygen at the cathode. The mass flux density of gas component r through the airway-membrane electrode interface in network units (i,j) and (i,j+1) are respectively. gc,z α represents the contact area between the airway and the membrane electrode in a single network unit. i,j α i,j-1 These are the pressure loss factors in network units (i,j) and (i,j-1), respectively, p i,j+1 p i,j p i,j-1Let β represent the pressure at the center of network cell (i,j+1), the center of (i,j), and the center of (i,j-1), respectively. i,j β i,j-1 These are the temperature factors at the center of network cell (i,j) and the center of network cell (i,j-1), respectively. h represents the total mass flux density of the mixed gas passing through the gas channel-membrane electrode interface at the centers of network units (i,j), (i,j+1), and (i,j-1), respectively. gb A is the convective heat transfer coefficient between the airflow and the bipolar plates in the air passage. gb h is the contact area between the airway and the bipolar plate in a single network unit. gc,z C is the convective heat transfer coefficient between the airflow in the gas channel and the boundary of the gas diffusion layer of the membrane electrode. p,k Let K be the specific heat capacity of gaseous component k, where k represents only hydrogen at the anode and only oxygen at the cathode. C represents the mass flux density of gas component k through the gas channel-membrane electrode interface in network units (i,j) and (i,j+1), respectively. p,r Let r be the specific heat capacity of the gas component. C represents the gas temperature at the center of network cell (i,j) and the center of network cell (i,j+1), respectively. p,v The specific heat capacity of water vapor. C represents the mass flux density of water vapor through the gas channel-membrane electrode interface in network units (i,j) and (i,j-1), respectively. p,l The specific heat capacity of liquid water, These represent the mass flux density of liquid water through the gas channel-membrane electrode interface in network units (i,j) and (i,j+1), respectively. These are the boundary temperatures of the membrane electrode gas diffusion layer in network cell (i,j) and (i,j+1), respectively. These represent the temperatures of the bipolar plates in network units (i,j), (i,j+1), and (i,j-1), respectively, and α is the temperature of the bipolar plate. bp k is the thermal diffusivity factor. bp Let A be the thermal conductivity of the bipolar plate. bp,z b represents the contact area between the bipolar plate and the membrane electrode in a single network unit. bp h represents the thickness of the bipolar plate along the z-axis. cool A is the convective heat transfer coefficient between the coolant and the bipolar plate. cool A represents the contact area between the bipolar plate and the coolant. bp,x Let A be the cross-sectional area of ​​the bipolar plate along the x-axis, Δx be the dimension of each network cell along the x-axis, and A be the cross-sectional area of ​​the bipolar plate along the x-axis. bp,y Let be the cross-sectional area of ​​the bipolar plate along the y-axis, and Δy be the dimension of each network cell along the y-axis. Let be the boundary temperature of the bipolar plate in network cell (i,j) at the contact surface with the membrane electrode. Let q represent the temperatures of the coolant in network cells (i,j), (i,j+1), and (i,j-1), respectively. cool C represents the coolant flow rate. p,cool This refers to the specific heat capacity of the coolant.

[0013] Optionally, the solution process for the gas pressure in the two-dimensional network model is as follows:

[0014] The expression for constructing the flow velocity in the airway is:

[0015]

[0016] The coordinate axes x', y', and z' are parallel to the x, y, and z axes, respectively. The range of x' is [0, Δx], the range of y' is [-a / 2, a / 2], and the range of z' is [-b / 2, b / 2]. x' Let x' be the gas velocity component along the x' axis, distributed across the cross-section of the airway, μ be the gas dynamic viscosity coefficient, and p be the gas pressure.

[0017] For velocity u x' Integrating over the cross-section of the DC channel and combining this with the perfect gas law pM = ρRT, the molar flow rate q is obtained. n The relationship between (x') and the pressure gradient is as follows:

[0018]

[0019] Where T is the gas temperature, Γ is an intermediate variable, and q n (x) represents the gas molar flow rate;

[0020] The molar flow rate q in the airway of network unit (i,j) n The first transformation is obtained by integrating the relationship between (x') and the pressure gradient:

[0021]

[0022] Where, p i,j p i,j+1 These represent the pressure at the center of network cell (i,j) and the center of network cell (i,j+1), respectively. Let A be the gas temperature at the center of network cell (i,j) and the center of network cell (i,j+1), respectively. gc,z This refers to the contact area between the airway and the membrane electrode in a single network unit. denoted as the total mass flux density of the mixed gas passing through the gas channel-membrane electrode interface at the center of network unit (i,j) and the center of (i,j+1), respectively. Let be the gas molar flow rate at the inlet of network cell (i,j);

[0023] The molar flow rate q in the airway of network unit (i,j-1) n The second transformation is obtained by integrating the equation relating (x') to the pressure gradient:

[0024]

[0025] Where, p i,j-1 The pressure at the center of network cell (i,j-1) is... Let (i,j-1) be the gas temperature at the center of the network cell. Let be the total mass flux density of the mixed gas at the center of network cell (i,j) through the gas channel-membrane electrode interface. Let $\frac{i}{j-1}$ be the gas molar flow rate at the inlet of network cell (i,j-1).

[0026] Subtract the first transformation from the second transformation and substitute the result into the airway flow control equation. The governing equations for airway pressure in network cell (i,j) are obtained as follows:

[0027]

[0028] Where, α i,j α i,j-1 These are the pressure loss factors in network units (i,j) and (i,j-1), respectively, β i,j β i,j-1 These are the temperature factors at the center of network unit (i,j) and the center of network unit (i,j-1), respectively.

[0029] Optionally, the one-dimensional differential model is:

[0030]

[0031] Where, σ p φ is the conductivity of the proton. p Let z be the proton potential, z be the coordinate along the thickness of the film electrode, and s be the coordinate along the thickness of the film electrode. p s e These are the source terms of the anode and cathode catalyst layers, respectively, σ e φ is the conductivity of electrons. e For the electron potential, j an and j ca These represent the electrochemical reaction rates at the anode and cathode, respectively, in seconds. l This represents the volume fraction of liquid water. These represent the anode and cathode exchange current densities, respectively; F is the Faraday constant; R is the gas constant; T is the temperature; and η is the estrus. an and η ca These are the overpotentials at the anode and cathode, respectively. and These are the four transfer coefficients in the Butler-Volmer model, V ocv ω is the open-circuit voltage of the proton exchange membrane fuel cell, ω is the volume fraction of the ion-polymer, and D is the open-circuit voltage. λ Let ρ be the diffusion coefficient of film water. mem Let be the density of the ionomer, EW be the equivalent molar mass of the ionomer, λ be the membrane water content, and i be the density of the ionomer. p For proton current density, s ad D is the source term for water adsorption / desorption by ionomers. s Let ρ be the diffusion coefficient of liquid water. l M is the density of liquid water. H2O Let be the molar mass of water, and s be the volume fraction of liquid water. ec For the liquid water source term, c represents the concentration of the mixed gas. Let x be the effective diffusion coefficient of gas component r. r Let s be the volume fraction of gaseous component r, where r represents hydrogen and water vapor at the anode and oxygen and water vapor at the cathode. r For gaseous component source terms, For the equivalent thermal conductivity, s T This is a heat absorption / heat release source term.

[0032] Optionally, the coupling of the two-dimensional network model and the one-dimensional differential model to establish a two-dimensional + one-dimensional mathematical model of multi-physics coupling in a proton exchange membrane fuel cell includes:

[0033] The two-dimensional network model is solved based on the preset physical boundary conditions and the second preset boundary conditions and the temperature field to update the gas pressure, humidity, gas components and liquid water volume fraction in the two-dimensional network model, and to serve as the first preset boundary conditions required to solve the one-dimensional differential model.

[0034] Based on the control equation and temperature field of the membrane electrode boundary temperature in the first preset boundary condition, the one-dimensional differential model is solved to update the distribution of electric potential, current density, each material component and endothermic / exothermic source terms in the membrane electrode, and to obtain the material flux density of each substance passing through the membrane electrode and the airway interface, which is used as the second preset boundary condition required to solve the two-dimensional network model.

[0035] The two-dimensional network model is solved based on the physical boundary conditions and the second preset boundary conditions to update the temperature distribution of the two-dimensional network model and obtain the temperature of the membrane electrode boundary, which serves as the first preset boundary condition required for solving the one-dimensional differential model.

[0036] Based on the control equation and energy equation of the membrane electrode boundary temperature in the first preset boundary condition, the heat absorption / heat release source term of the energy equation is used to solve the one-dimensional differential model to update the temperature distribution in the membrane electrode and obtain the heat flux density at the airway-membrane electrode interface, which serves as the second preset boundary condition for the two-dimensional network model.

[0037] The above solution process is iterated, continuously updating the current density, the electric potential, the temperature, the gas pressure, the humidity, the volume fraction of each gas component and liquid water, the mass flux density and heat flux density of the gas channel-membrane electrode interface, until the iteration converges, thus obtaining the two-dimensional + one-dimensional mathematical model of the multi-physics coupling of the proton exchange membrane fuel cell.

[0038] A second aspect of the present invention provides a device for predicting the performance of a proton exchange membrane fuel cell, comprising: a partitioning module for dividing a target monolithic fuel cell into multiple network units, each network unit including a bipolar plate, a gas channel, a membrane electrode assembly, and coolant on the upper and lower sides for heat dissipation; a construction module for establishing the relationship between the multiple network units, constructing a two-dimensional network model in a two-dimensional plane direction in combination with a first preset physical boundary condition, and establishing a one-dimensional differential model in the membrane electrode assembly of each network unit; a coupling module for coupling the two-dimensional network model and the one-dimensional differential model to establish a two-dimensional + one-dimensional mathematical model of multi-physics coupling of the proton exchange membrane fuel cell; and a prediction module for using the two-dimensional + one-dimensional mathematical model of multi-physics coupling of the proton exchange membrane fuel cell to solve for the distribution of current density, pressure, flow rate, temperature, humidity, and gas concentration to predict the performance of the target monolithic fuel cell.

[0039] A third aspect of the present invention provides an electronic device, including: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the proton exchange membrane fuel cell performance prediction method as described in the above embodiments.

[0040] A fourth aspect of the present invention provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the above-described method for predicting the performance of a proton exchange membrane fuel cell.

[0041] The method and apparatus for predicting the performance of a proton exchange membrane fuel cell according to embodiments of the present invention have the following beneficial effects:

[0042] A single fuel cell is divided into a series of network units, and a two-dimensional network model is established. For each network unit, a one-dimensional differential model of the membrane electrode structure is established, and a multi-physics coupled "two-dimensional + one-dimensional" mathematical model is established, which can simulate and calculate the electrochemical reaction process, the flow and mass transfer process in the gas channel and porous medium, and the heat transfer (heat dissipation) process in the fuel cell.

[0043] The numerical calculation process of bipolar plates is simplified, so that computational resources are relatively concentrated on the numerical calculation part of membrane electrode, thereby obtaining high-resolution physical quantity distribution characteristics, while effectively limiting the calculation time, which meets the design requirements of fuel cell membrane electrode.

[0044] The computational process of coupling the two-dimensional network model and the one-dimensional differential model has been verified through actual operation and can stably achieve the computational function.

[0045] Additional aspects and advantages of the invention will be set forth in part in the description which follows, and in part will be obvious from the description, or may be learned by practice of the invention. Attached Figure Description

[0046] The above and / or additional aspects and advantages of the present invention will become apparent and readily understood from the following description of the embodiments taken in conjunction with the accompanying drawings, wherein:

[0047] Figure 1 A flowchart illustrating a method for predicting the performance of a proton exchange membrane fuel cell according to an embodiment of the present invention;

[0048] Figure 2 This is a schematic diagram of a two-dimensional network model provided according to an embodiment of the present invention;

[0049] Figure 3 This is a schematic diagram of the network unit provided in the embodiment of the present invention in the z-axis direction;

[0050] Figure 4 This is a heat transfer path analysis diagram of each component in a two-dimensional network model provided according to an embodiment of the present invention;

[0051] Figure 5 This is a schematic diagram of a one-dimensional differential model provided according to an embodiment of the present invention;

[0052] Figure 6 This is a schematic diagram illustrating the coupling between a two-dimensional network model and a one-dimensional differential model provided according to an embodiment of the present invention.

[0053] Figure 7 This is a flowchart illustrating the coupled computation of a two-dimensional network model and a one-dimensional differential model according to an embodiment of the present invention.

[0054] Figure 8 This is a schematic diagram of the average current density curves of a fuel cell under different voltage states according to an embodiment of the present invention;

[0055] Figure 9 This is a schematic diagram of the current density distribution inside a fuel cell when the voltage is 0.7V, according to an embodiment of the present invention.

[0056] Figure 10 This is a schematic diagram showing the distribution of relative humidity in the cathode gas channel at a voltage of 0.7V, according to an embodiment of the present invention.

[0057] Figure 11 This is a block diagram of a proton exchange membrane fuel cell performance prediction device according to an embodiment of the present invention;

[0058] Figure 12 This is a schematic diagram of the structure of an electronic device according to an embodiment of the present invention. Detailed Implementation

[0059] Embodiments of the present invention are described in detail below. Examples of these embodiments are illustrated in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and intended to explain the present invention, and should not be construed as limiting the present invention.

[0060] The following describes, with reference to the accompanying drawings, a method and apparatus for predicting the performance of a proton exchange membrane fuel cell according to an embodiment of the present invention.

[0061] Figure 1 This is a flowchart illustrating a method for predicting the performance of a proton exchange membrane fuel cell, as provided in an embodiment of the present invention.

[0062] like Figure 1 As shown, the method for predicting the performance of a proton exchange membrane fuel cell includes the following steps:

[0063] In step S101, the target monolithic fuel cell is divided into multiple network units. Each network unit includes a bipolar plate, a gas channel, a membrane electrode assembly, and coolant on the upper and lower sides for heat dissipation.

[0064] Specifically, the two-dimensional plane containing the target monolithic fuel cell is set as the xy plane, with the positive x-axis representing the anode gas flow direction and the z-axis representing the thickness direction of the target monolithic fuel cell; for example... Figure 2 As shown, the two-dimensional plane of the target monolithic fuel cell is divided. Let N be the number of gas channels with the bipolar plates of the fuel cell arranged side by side along the y-axis. y Then divide N along the y-axis direction yThere are 1 network element with a size of Δy; the network element size is set to Δx in the x-axis direction based on the simulation accuracy requirements, and the number of network elements N is given in combination with the total length of the fuel cell. x Among them, such as Figure 3 As shown, the structure of each network unit includes a bipolar plate, a gas channel, a membrane electrode, and coolant on the upper and lower sides for heat dissipation.

[0065] In some embodiments of the present invention, a layer of boundary virtual point grid cells with the same size as the network cells is set on the outer wall surface of the divided target monolithic fuel cell to set the equations corresponding to the boundary conditions. Specifically, as shown in the figure... Figure 2 As shown, the actual area occupied by the fuel cell is the white grid area, and the thick black solid line is the side wall of the fuel cell, that is, the boundary of the two-dimensional network model.

[0066] To standardize the governing equations within the computational domain and reduce program modifications when boundary conditions are altered in the future, a layer of identically sized boundary virtual point network elements was added around the fuel cell computational domain, as shown in the shaded elements in the figure. The boundary conditions defined on the bold black solid lines were then transformed into governing equations for the corresponding variables on the boundary virtual point network elements. At this point, for the computational domain of the two-dimensional network model, N exists in the x-axis direction. x +2 units, N in the y-axis direction y +2 units. i = 1 or N y +2 and j = 1 or N x +2 represents the boundary virtual point network unit.

[0067] In step S102, the relationship between multiple network units is established, a two-dimensional network model in the two-dimensional plane direction is constructed in combination with the first preset physical boundary conditions, and a one-dimensional differential model is established in the membrane electrode of each network unit.

[0068] Specifically, in this embodiment of the invention, based on the principles of conservation of mass, momentum, and energy, relationships between network units are established in the xy-plane, and a two-dimensional network model in the xy-plane direction is constructed by combining physical boundary conditions. Then, in each network unit, a one-dimensional differential model is established in the membrane electrode. For example... Figure 4 As shown, the membrane electrode consists of a gas diffusion layer, a catalyst layer, and a proton exchange membrane. One-dimensional differential conservation equations and electrochemical reaction equations are established for each structure to solve for the distribution of each physical quantity along the z-axis.

[0069] The two-dimensional network model in the xy plane is as follows:

[0070] First, assuming that both the anode and cathode gases flow along the positive x-axis, the flow control equations for the anode / cathode gas channels are constructed as follows:

[0071]

[0072] In the formula, denoted as r, the flow rates of component r in the mixed gas at the center of network cell (i,j) and the center of network cell (i,j+1), respectively. r represents hydrogen and water at the anode and oxygen at the cathode. The mass flux density of gas component r through the airway-membrane electrode interface in network units (i,j) and (i,j+1) are respectively. gc,z This represents the contact area between the airway and the membrane electrode in a single network unit. It should be noted that to solve this system of equations, appropriate boundary conditions need to be set at the airway inlet / outlet. The method for setting boundary conditions for the two-dimensional network model will be described in a unified manner later.

[0073] After determining the flow distribution within the airway, the proportion of each component can be calculated. The volume fraction of each gas component in the airway of network cell (i,j) is the component flow rate. Total flow rate of the mixed gas The proportion of.

[0074] The volume fraction of liquid water s in the airway of network unit (i,j) i,j The expression is:

[0075]

[0076] The relative humidity in the airway of network unit (i,j) is:

[0077]

[0078] Among them, in formulas (2) and (3) The expression is:

[0079]

[0080] In the formula, Let be the volume fraction of water in the airway of network cell (i,j). p represents the flow rates of water vapor and liquid water in the air passage of network unit (i,j). i,j Let be the gas pressure in the airway of network cell (i,j). The saturated vapor pressure in the airway of network unit (i,j) is obtained by solving the subsequent formula (30), ρ. l Let R be the density of liquid water, and R be the gas constant. Let be the temperature in the airway of network cell (i,j). Let be the molar mass of water. It should be noted that the volume fractions of each gas component, the volume fraction of liquid water, and the relative humidity in the aforementioned airway are the first preset boundary conditions provided to the one-dimensional network model in this paper.

[0081] Then, the pressure control equations for the anode / cathode gas passages are constructed. For a rectangular cross-section direct-flow passage with width *a*, height *b*, and length *Δx*, the internal gas flow process is analyzed. Under the constraints of the operating parameters within the fuel cell gas passage, the relationship between gas flow rate and pressure can be calculated from the analytical solution of the fluid control equations, the Navier-Stokes equations. The expression for the flow velocity in the gas passage is derived mathematically as follows:

[0082]

[0083] The coordinate axes x', y', and z' are parallel to the x, y, and z axes, respectively. The range of x' is [0, Δx], the range of y' is [-a / 2, a / 2], and the range of z' is [-b / 2, b / 2]. x' Let be the gas velocity component along the x' axis, which has a certain distribution across the cross-section of the gas passage, and its specific value depends on the coordinates y' and z'. μ is the gas dynamic viscosity coefficient, and p is the gas pressure. For velocity u... x' Integrating over the cross-section of the DC channel and combining this with the perfect gas law pM = ρRT, the molar flow rate q can be obtained. n The relationship between (x') and the pressure gradient is as follows:

[0084]

[0085] In the formula, T is the gas temperature, Γ is an intermediate variable, and q n (x) represents the gas molar flow rate.

[0086] The expression for the intermediate variable Γ is:

[0087]

[0088] function The expression is:

[0089]

[0090] In the two-dimensional network model, the governing equation (6) for gas pressure p includes two sets of unknown variables—gas pressure p and flow rate q. nThis requires solving for these two sets of variables in the flow field using the airflow control equation. However, the boundary conditions of the airflow control equation are usually set at the upstream inlet, while the boundary conditions of the gas pressure control equation are usually set at the downstream outlet. Therefore, iterative algorithms need to be designed to advance the calculation process that satisfies the boundary conditions, which carries the risk of slow iterative convergence or poor computational stability.

[0091] Equation (6) is transformed in the following way to obtain the governing equation for independently solving the gas pressure p, thus decoupling the calculation process of gas flow rate and pressure. The integral in Equation (2) is calculated in the gas passage of network unit (i,j), and then rearranged to obtain:

[0092]

[0093] Similarly, the integral in formula (6) is calculated in the airway of network unit (i,j-1), and the result is:

[0094]

[0095] Subtracting equations (9) and (10) and substituting the result into the airway flow control equation, and combining this with the mass conservation equation, eliminates the variable flow rate q in the pressure control equation. n Thus, the governing equation for airway pressure in network cell (i,j) can be obtained:

[0096]

[0097] In the formula, the superscript (i,j) in the variable symbol indicates that the variable belongs to network unit (i,j), α i,j p is the pressure loss factor in network cell (i,j). i,j β is the pressure at the center of network cell (i,j). i,j Temperature factor Let be the total mass flux density of the mixed gas in network cell (i,j) through the gas channel-membrane electrode interface. Let be the gas molar flow rate at the inlet of network cell (i,j-1).

[0098] Formula (11) in each network unit forms the governing equations for solving the gas pressure within the airway. It should be noted that formula (11) couples the mass and momentum conservation equations within a single grid network unit, eliminating the flow rate variable to obtain the governing equations for gas pressure. This differs from the traditional approach of establishing gas pressure governing equations in a closed loop where the gas pressure drop is zero. It should also be noted that, in this embodiment of the invention, to solve this governing equation set, corresponding boundary conditions need to be set at the airway inlet / outlet. The method for setting boundary conditions in the two-dimensional network model will be described uniformly later.

[0099] The coefficients in formula (11) are expressed as follows:

[0100]

[0101] In the formula, a and b are the width and height of the airway cross section, respectively, μ is the gas dynamic viscosity coefficient, and n is the sequence number of the series summation in formula (12).

[0102] The expression for the gas dynamic viscosity coefficient μ in formula (12) is as follows:

[0103]

[0104] In the formula, r and s represent the material components, referring to hydrogen, water vapor, and liquid water at the anode, and oxygen, nitrogen, water vapor, and liquid water at the cathode, and x r or x s μ represents the volume fraction of gaseous component r or s. r or μ s M is the kinetic viscosity coefficient of bulk component r or s. s Or M r Let r be the molar mass of the gaseous component r or s.

[0105] Next, the temperature field is calculated. In the two-dimensional network model, the heat transfer paths between the membrane electrode, gas channels, bipolar plates, and coolant are as follows: Figure 4 As shown. The central membrane electrode continuously releases heat during operation, outputting heat flow to both the anode and cathode sides. Taking the anode as an example, at the boundary of the membrane electrode, part of the heat flow is transferred to the bipolar plate via solid-state thermal conduction (there is contact thermal resistance between the membrane electrode and the bipolar plate, therefore the membrane electrode boundary temperature...). Bipolar plate boundary temperature (The amounts of heat are not equal); another portion of the heat flow is transferred to the air passage via convection, changing the gas temperature. Simultaneously, the airflow in the air passage transfers heat to the bipolar plate walls via convection. Furthermore, the coolant at the top of the bipolar plate dissipates heat from the bipolar plate via forced convection. The sidewalls of the bipolar plate are exposed to the environment, where natural convection heat transfer occurs.

[0106] The temperature control equation for the gas channel is constructed below based on the energy conservation equation. Since the gas exhibits one-dimensional flow characteristics, the heat transfer process is primarily convection-based, neglecting the gas's conduction effect. Simultaneously, the convective heat transfer processes between the gas and the membrane electrode boundary and the bipolar plate wall must also be considered. Based on this, a one-dimensional temperature control equation can be constructed for network unit (i,j):

[0107]

[0108] Formula (14) in each network unit forms the governing equations for solving the gas temperature in the airway. In the formula, the subscripts r and k indicate that the variable belongs to gas component r or k (at the anode, r represents hydrogen and water, and k represents only hydrogen; at the cathode, r represents oxygen, nitrogen, and water, and k represents only oxygen). gb A is the convective heat transfer coefficient between the airflow and the bipolar plates in the air passage. gb h is the contact area between the airway and the bipolar plate in a single network unit. gc,z C is the convective heat transfer coefficient between the airflow in the gas channel and the boundary of the gas diffusion layer of the membrane electrode. p,k Let k be the specific heat capacity of the gas component. Let C be the mass flux density of gas component k through the gas channel-membrane electrode interface in network cell (i,j). p,r Let r be the specific heat capacity of the gas component. Let C be the gas temperature in the airway of network cell (i,j). p,v The specific heat capacity of water vapor. Let C be the mass flux density of water vapor through the gas channel-membrane electrode interface in network cell (i,j). p,l The specific heat capacity of liquid water, Let be the mass flux density of liquid water through the gas channel-membrane electrode interface in network cell (i,j). The boundary temperature of the gas diffusion layer of the membrane electrode in network cell (i,j) is given. Let be the temperature of the bipolar plate in network cell (i,j). It should be noted that to solve this system of equations, appropriate boundary conditions need to be set at the airway inlet / outlet. The method for setting boundary conditions for the two-dimensional network model will be described uniformly later.

[0109] Next, the temperature control equation for the bipolar plate is established. The temperature distribution within the bipolar plate is mainly controlled by solid-state heat conduction. Convective heat transfer occurs at the boundary with the airflow, cooling water, and the natural environment, while heat conduction occurs with the gas diffusion layer of the membrane electrode. Furthermore, there is contact thermal resistance between the bipolar plate and the gas diffusion layer. Since the bipolar plate is very thin, its characteristic dimension in the two-dimensional plane is much larger than its thickness. Therefore, the temperature change along the thickness direction of the bipolar plate is neglected, and the main problem to solve is the temperature distribution in the two-dimensional plane. The bipolar plate temperature is derived from the energy conservation equation in network element (i,j). The governing equations are:

[0110]

[0111] Formula (15) in each network unit forms the governing equations for solving the bipolar plate temperature. In the formula, α bp k is the thermal diffusivity factor. bp Let A be the thermal conductivity of the bipolar plate. bp,zb represents the contact area between the bipolar plate and the membrane electrode in a single network unit. bp h represents the thickness of the bipolar plate along the z-axis. cool A is the convective heat transfer coefficient between the coolant and the bipolar plate. cool A represents the contact area between the bipolar plate and the coolant. bp,x Let N be the cross-sectional area of ​​the bipolar plate along the x-axis, and Δx be N. x The size of each network unit, A bp,y Let N be the cross-sectional area of ​​the bipolar plate along the y-axis, and Δy be N. y The size of each network unit Let be the boundary temperature of the bipolar plate in network cell (i,j) at the contact surface with the membrane electrode. Let α be the temperature of the coolant in network cell (i,j). bp The expression is It should be noted that in order to solve this system of equations, it is also necessary to set corresponding boundary conditions on the bipolar plate boundary. The method of setting boundary conditions for the two-dimensional network model will be described in a unified manner later.

[0112] Next, the temperature control equation for the coolant is established. The flow direction of the coolant is determined by the actual cooling channel structure and needs to be set according to the actual situation. Here, it is assumed that the flow direction of the coolant is the same as that of the gas in the anode channel. Since the coolant absorbs heat from the bipolar plates through convection heat transfer, the diffusion term in the energy equation can be ignored, simplifying it to a one-dimensional equation, resulting in the temperature control equation for the coolant:

[0113]

[0114] Formula (16) in each network unit forms the governing equations for solving the coolant temperature. Where q cool C represents the coolant flow rate. p,cool This represents the specific heat capacity of the coolant. It should be noted that, to solve this system of equations, appropriate boundary conditions need to be set at the inlet / outlet of the coolant passage.

[0115] The following explains the method for setting boundary conditions in a two-dimensional network model. When solving the governing equations for bipolar plate temperature, airflow, pressure, and temperature using a two-dimensional network model, boundary conditions for the corresponding variables need to be set on the boundaries of the model. In this model, the expressions on the boundaries are transformed into governing equations within the virtual point network elements of the boundaries.

[0116] For the airway flow control equation, one additional boundary condition is needed to achieve a boundary value solution. Since gas convection effects are dominant, this can be directly applied to the network elements (i = 2 ~ N). y Set the inlet flow rate of each component in (j=1).

[0117]

[0118] Because gas convection is the dominant effect, the boundary virtual point network units (i = 2 ~ N) at the airway outlet are... y j=N x The value of +2) will not affect the flow rate in the inner domain. However, to maintain the integrity of the model, a fully developed boundary condition for the gas flow can be specified at the airway outlet:

[0119]

[0120] For the airway pressure governing equation, two additional boundary conditions are needed to achieve a definite solution. Typically, these are applied at the airway inlet (i = 2 ~ N). y (j=1) Determined by inlet flow Controlling the pressure at the outlet (i = 2 ~ N) y j=N x +2) Directly specify the gas pressure (back pressure).

[0121] At the airway inlet, let the formula (11) Obtain network units (i = 2 ~ N) y The governing equation at (j=1) is

[0122]

[0123] At the airway outlet (i = 2 ~ N) y j=N x +2), assuming the gas temperature has fully developed, that is And by the specified outlet gas pressure p out Replace the variable p in formula (11) i,j+1 The derivation of formula (11) is repeated to obtain the boundary virtual point network unit (i = 2 ~ N). y j=N x The governing equation for +2) is

[0124]

[0125] For the temperature control equation in the airway, one additional boundary condition is needed to achieve a boundary value solution. Since gas convection is the dominant effect, it can be directly applied to the network elements (i = 2 ~ N). y Set the inlet temperature in (j=1)

[0126]

[0127] Because gas convection is the dominant effect, the boundary virtual point network units (i = 2 ~ N) at the airway outlet are... y j=Nx The +2) value will not affect the temperature in the inner domain. However, to maintain the integrity of the model, a fully developed boundary condition for the gas temperature can be specified at the airway outlet:

[0128]

[0129] For the two-dimensional temperature control equation of a bipolar plate, four boundary conditions need to be added. Natural convection heat transfer conditions are typically set on the four sidewalls of the bipolar plate. These conditions are determined by network elements (j=2 or N). x +1, i = 2 or N y Linear extrapolation of the bipolar plate temperature and wall temperature (+1) yields boundary network elements (j=1 or N) that satisfy the natural convection heat transfer boundary conditions. x +2, i = 1 or N y +2) Bipolar plate temperature for:

[0130]

[0131] In the formula, For the network elements of the neighboring element (i,j) in the x-axis direction Belongs to network unit j=2 or N x +1, in the y-axis direction Belongs to network unit i = 2 or N y +1) temperature, T amb For ambient temperature, h bc Let Δl be the convective heat transfer coefficient between the bipolar plate sidewall and the environment. For the boundary along the x-axis, Δl = Δx, and for the boundary along the y-axis, Δl = Δy.

[0132] For the temperature control equation of the coolant, one boundary condition needs to be added. Since the coolant convection effect is dominant, it can be directly applied to the network elements (i = 2 ~ N). y Set the inlet temperature in (j=1)

[0133]

[0134] Since the coolant convection effect is dominant, the boundary virtual point network elements (i = 2 ~ N) at the coolant outlet are... y j=N x The value of +2) will not affect the temperature in the inner domain. However, to maintain the integrity of the model, a fully developed temperature boundary condition can be specified at the coolant outlet:

[0135]

[0136] In summary, a two-dimensional network model in the two-dimensional plane of the target single fuel cell is constructed according to formulas (1)-(16):

[0137]

[0138] Formulas (17)-(25) constitute the boundary conditions for the definite solution of the two-dimensional network model.

[0139] Furthermore, such as Figure 5 As shown, for any network unit, the one-dimensional differential model equation established in the membrane electrode is:

[0140] First, construct the charge conservation equations for the proton exchange membrane, catalytic layer, and gas diffusion layer:

[0141]

[0142] In the formula, σ p φ is the conductivity of the proton. p Let z be the proton potential, z be the coordinate along the thickness of the film electrode, and s be the coordinate along the thickness of the film electrode. p s e These are the electrochemical reaction source terms for the anode / cathode catalyst layers, σ e φ is the conductivity of electrons. e This is the electron potential.

[0143] In formula (26), the source term s p s e The expressions for the catalyst layer are listed in the table below:

[0144]

[0145] In the gas diffusion layer and proton exchange membrane, the source term s p s e All are 0, j an j ca These are the cation / cathode electrochemical reaction rates, expressed as:

[0146]

[0147] In the formula, s l The volume fraction of liquid water is determined by the liquid water transport equation. The anode exchange current density, and These are the four transfer coefficients in the Butler-Volmer model, where F is the Faraday constant, R is the gas constant, T is the temperature, and η is the overpotential. In the anode catalyst layer, the overpotential η = φ e -φ p In the cathode catalyst layer, η = φ e-φ p -V ocv V ocv It is the open-circuit voltage of a proton exchange membrane fuel cell. This represents the cathode exchange current density.

[0148] Among them, the anode exchange current density and cathode exchange current density The expression is:

[0149]

[0150] In the formula, A an P is the anode reverse reference current density. H2 For the partial pressure of hydrogen, P ref For reference pressure, E an R is the anodic activation barrier, R is the gas constant, and T is the anodic activation barrier. ref Here, T is the reference temperature, and A is the actual temperature. ca P is the cathode reference current density. O2 E represents the partial pressure of oxygen. ca This is the cathode activation barrier.

[0151] σ in formula (20) p σ e It is the conductivity of protons and electrons, where,

[0152]

[0153] In the formula, ω is the volume fraction of the ionomer, λ is the membrane water content, and e is the natural logarithm.

[0154] Then, the governing equations for the membrane water content λ in the catalyst layer and proton exchange membrane are constructed:

[0155]

[0156] In the formula, z is the coordinate along the thickness direction of the membrane electrode, and D λ Let ρ be the diffusion coefficient of film water. mem Let EW be the volume fraction of the ionomer, EW be the equivalent molar mass of the ionomer, and F be the Faraday constant. The proton current density is given by formula (26), and the proton potential φ is obtained by solving the formula. p Decision, s ad The conversion rate is denoted as water vapor to membrane water.

[0157] Among them, the diffusion coefficient D of film water λ for:

[0158]

[0159] In formula (30), s ad The source term for water adsorption / desorption by ionomers is expressed as follows:

[0160]

[0161] In the formula, κ ad λ represents the adsorption / desorption conversion rate. eq To balance the water content of the membrane.

[0162] Among them, the water content λ of the equilibrium membrane eq The expression is:

[0163] λ eq =0.043 + 17.81a - 39.85a 2 +36.0a 3 (33)

[0164] In the formula, a = RH + 2s l RH represents the activity of water vapor.

[0165] Among them, the adsorption / desorption conversion rate κ ad The expression is:

[0166]

[0167] In the formula, κ ad0 Intrinsic adsorption / desorption conversion rate, V m Let be the volume fraction of water in the film, h be the thickness of the catalyst layer, R be the gas constant, and T be the volume fraction of water in the film. ref T is the reference temperature, and T is the actual temperature.

[0168] Among them, the intrinsic adsorption / desorption conversion rate κ ad0 The expression is:

[0169]

[0170] Next, the governing equations for the volume fraction of liquid water in the gas diffusion layer and the catalyst layer are constructed as follows:

[0171]

[0172] In the formula, D s Let ρ be the diffusion coefficient of liquid water. l The density of liquid water, Let be the molar mass of liquid water, and s be the volume fraction of liquid water. ec This is a water vapor / liquid water phase change source term.

[0173] Among them, the diffusion coefficient D of liquid water s The expression is:

[0174]

[0175] In the formula, K0 is the inherent permeability of liquid water, and s red To correct for the liquid water volume fraction, μ l ρ is the dynamic viscosity coefficient of liquid water, and p is the gas pressure in the porous medium.

[0176] Among them, the corrected liquid water volume fraction s red The expression is:

[0177]

[0178] In the formula, s im This refers to the amount of stagnant water in the gas diffusion layer and the catalyst layer.

[0179] Capillary pressure driving force The expression is:

[0180]

[0181] In the formula, θ is the contact angle of liquid water in the gas diffusion layer and the catalyst layer.

[0182] In formula (36), the water vapor / liquid water phase change source term s ec The expression is:

[0183]

[0184] In the formula, γ ec Let be the phase transition rate, and c be the concentration of the mixture. x represents the volume fraction of water. sat This represents the volume fraction of saturated water vapor.

[0185] Among them, the phase transition rate γ ec The expression is:

[0186]

[0187] In the formula, Let R be the molar mass of water, R be the gas constant, and T be the temperature.

[0188] Wherein, saturated water vapor volume fraction x sat The expression is:

[0189] log 10 x sat =-2.1794+0.02953T-9.1837×10 -5 T 2 +1.4454×10 -7 T 3(42)

[0190] Next, the governing equations for the gas components in the gas diffusion layer and the catalyst layer are constructed:

[0191]

[0192] In the formula, z is the coordinate along the thickness direction of the membrane electrode, and c is the concentration of the mixture. Let x be the effective diffusion coefficient of gas component r. r s is the volume fraction of gaseous component r. r This represents the source term for the generation / consumption of gaseous component r during an electrochemical reaction or phase transition. The subscript r represents hydrogen and water vapor at the anode, and oxygen and water vapor at the cathode.

[0193] Wherein, the effective diffusion coefficient of gas component r The expression is:

[0194]

[0195] In the formula, Let ε be the intrinsic diffusion coefficient of gas component r, ε be the porosity of the gas diffusion layer / catalytic layer, and s be the intrinsic diffusion coefficient of gas component r. l T represents the volume fraction of liquid water. ref The reference temperature is T, the actual temperature is p. ref p is the reference pressure, where p is the gas pressure.

[0196] In formula (43), the source term s for the generation / consumption of gaseous component r during electrochemical reaction or phase transition is... r Related to electrochemical reactions and phase transition rates, the specific expressions are listed in the table below:

[0197]

[0198] Finally, the temperature control equations for each component are constructed as follows:

[0199]

[0200] In the formula, z is the coordinate along the thickness direction of the membrane electrode. The equivalent thermal conductivity is given by T, where T is the temperature and s is the temperature. T This is a heat absorption / heat release source term.

[0201] Among them, the heat absorption / heat release source item s T The expressions for each component are:

[0202]

[0203] In the formula, H ec The latent heat of the water vapor-liquid water phase transition, S FFor the electrochemical reaction rate, ΔS HOR This is the anodic hydrogen oxidation reaction, H ad The latent heat of adsorption / desorption of water by ionomers, ΔS ORR The reaction entropy is the oxygen reduction reaction at the cathode.

[0204] In summary, based on formulas (26)-(45), the one-dimensional differential model of the membrane electrode is constructed as follows:

[0205]

[0206] After constructing the two-dimensional network model and the one-dimensional differential model, the two models need to be coupled to establish a "two-dimensional + one-dimensional" mathematical model of multi-physics coupling of proton exchange membrane fuel cell. Then, the distribution of physical quantities such as current density, temperature, gas pressure, flow rate, humidity, and volume fraction of each component can be solved to realize the function of rapid prediction of fuel cell performance. The specific content is described below.

[0207] In step S103, a two-dimensional network model and a one-dimensional differential model are coupled to establish a two-dimensional + one-dimensional mathematical model of multi-physics coupling for a proton exchange membrane fuel cell.

[0208] In step S104, the distribution of current density, temperature, gas pressure, flow rate, humidity, and volume fraction of each component is solved using a two-dimensional + one-dimensional mathematical model of multi-physics coupling of proton exchange membrane fuel cell to predict the performance of the target single fuel cell.

[0209] Furthermore, in one embodiment of the present invention, a two-dimensional network model and a one-dimensional differential model are coupled to establish a two-dimensional + one-dimensional mathematical model of multi-physics coupling in a proton exchange membrane fuel cell, including:

[0210] The temperature distribution solution module provided by the temperature field module solves the two-dimensional network model based on the preset physical boundary conditions and the second preset boundary conditions and the temperature field, so as to update the gas pressure, humidity, gas components and liquid water volume fraction in the two-dimensional network model, and serve as the first preset boundary conditions required to solve the one-dimensional differential model.

[0211] Based on the control equation and temperature field of the membrane electrode boundary temperature in the first preset boundary condition, a one-dimensional differential model is solved to update the distribution of electric potential, current density, material composition and endothermic / exothermic source terms in the membrane electrode, and obtain the material flux density of each material through the membrane electrode and the gas channel interface, which is used as the second preset boundary condition required to solve the two-dimensional network model.

[0212] The module for solving the temperature field based on the endothermic / exothermic source terms of the energy equation provided by the basic module: It solves the two-dimensional network model based on physical boundary conditions and second preset boundary conditions to update the temperature distribution of the two-dimensional network model and obtain the temperature of the membrane electrode boundary, which serves as the first preset boundary condition required for solving the one-dimensional differential model.

[0213] Based on the control equation and energy equation of the membrane electrode boundary temperature in the first preset boundary condition, the one-dimensional differential model is solved to update the temperature distribution inside the membrane electrode and obtain the heat flux density at the gas channel-membrane electrode interface, which serves as the second preset boundary condition for the two-dimensional network model.

[0214] The above solution process is iterated, continuously updating the current density, electric potential, temperature, gas pressure, humidity, volume fraction of each gas component and liquid water, mass flux density and heat flux density at the gas channel-membrane electrode interface, until the iteration converges, resulting in a two-dimensional + one-dimensional mathematical model of multi-physics coupling of the proton exchange membrane fuel cell.

[0215] Specifically, solving the two-dimensional + one-dimensional mathematical model of multi-physics coupling in a proton exchange membrane fuel cell requires setting physical boundary conditions, including gas flow rate, hydrogen and oxygen ratio, relative humidity, and temperature at the anode / cathode inlet; gas pressure at the outlet; coolant flow rate and inlet temperature; ambient temperature; and operating voltage U. The one-dimensional differential model of the membrane electrode assembly (MEA) and the two-dimensional network model are coupled using the following method to form the two-dimensional + one-dimensional mathematical model of multi-physics coupling in the proton exchange membrane fuel cell.

[0216] like Figure 6 As shown, for the physical quantities in the basic module, the two-dimensional network model calculates and outputs the gas pressure, humidity, volume fraction of each gas component and liquid water, which serve as the first preset boundary conditions required for solving the one-dimensional differential model equations; the one-dimensional differential model calculates and outputs the mass flux density of each substance passing through the membrane electrode and gas channel interface, which serves as the interface mass flux density in formulas (1), (11), and (14) of the two-dimensional network model. and

[0217] For the temperature field module, the boundary temperature of the differential single-pool model is calculated. The governing equations are:

[0218]

[0219] Among them, h gc,z A is the convective heat transfer coefficient between the gas in the gas channel and the gas diffusion layer boundary of the membrane electrode. gc,z R represents the contact area between the airway and the membrane electrode in the network unit. ct The contact thermal resistance between the bipolar plate and the membrane electrode. Let (i,j) be the boundary temperature of the membrane electrode in network cell (i,j). It is the temperature of the bipolar plate in network cell (i,j) at the boundary where it contacts the membrane electrode. This represents the heat flux density through the gas channel-membrane electrode interface in network cell (i,j). It should be noted that the subscript k represents hydrogen at the anode and oxygen at the cathode.

[0220] Solving formula (46) in the two-dimensional array network model can output the boundary temperature required for solving the differential single-pool model. The differential single-cell model can output the heat flux density at the airway-membrane electrode interface. Formula (46) is used as the boundary condition for solving the temperature field of the membrane electrode using a one-dimensional differential model.

[0221] Finally, the coupling method between the flow field and the temperature field is as follows: solve the model equations involved in the basic module, and output the gas pressure, humidity, volume fraction of each component, and endothermic / exothermic source terms of the energy equation; based on multiple physical quantities, solve the model equations involved in the temperature field module, and output the various temperature values ​​required to solve the basic module.

[0222] Based on the above method, the two-dimensional network model and the one-dimensional differential model of the membrane electrode can be coupled to obtain a complete two-dimensional + one-dimensional mathematical model of the proton exchange membrane fuel cell. Simultaneously, based on the above coupling method, a specific solution program can be implemented, such as... Figure 7 After the calculation begins, 1) set the geometric dimensions of the fuel cell to be simulated, the physical properties required by the model, the physical boundary conditions, the number of two-dimensional array network elements, and the number of differential single cells; 2) set the parameters required for numerical calculation, such as the allowable residual ε0 at iteration convergence and the maximum number of iterations n. max 3) Set the initial values ​​for each physical quantity in the model for iteration. 4) Solve equations (1)-(13) in the two-dimensional network model based on boundary condition equations (17)-(20), and update the gas pressure, humidity, flow rate, and volume fraction of each component in the two-dimensional array network model; 5) Based on the updated gas pressure, humidity, flow rate, volume fraction of each component, and boundary temperature set by equation (46), solve equations (26)-(45) in the one-dimensional differential single-cell model, and update the distribution of each physical quantity in the membrane electrode, the mass flux density and heat flux density of each component passing through the gas channel-membrane electrode interface; 6) Based on the updated mass flux density, heat flux density, and boundary condition equations (21)-(25), solve equations (14)-(16) in the two-dimensional network model, and update the temperature value in the two-dimensional network model. At this point, all physical quantities in the model have been updated to... 7) Calculate the iterative residuals of each physical quantity. The calculation ends when the convergence condition ε < ε0 is met; otherwise, proceed to step 4) to begin the (n+1)th iteration of the calculation.

[0223] In some embodiments, steps 5) and 6) can also be replaced by the following method:

[0224] 5) Based on the updated gas pressure, humidity, flow rate and volume fraction of each component, solve the formulas (26)-(44) in the one-dimensional differential single-cell model, update the distribution of each physical quantity in the membrane electrode, and the source term s of the energy equation. T The mass flux density of each component through the airway-membrane electrode interface;

[0225] 6) Couple the solution of the whole field temperature equations (14)-(16), (45) and (46) and update the temperature values.

[0226] The following specific embodiment further illustrates the two-dimensional + one-dimensional mathematical model of multi-physics coupling in proton exchange membrane fuel cells proposed in this invention.

[0227] The performance of a 5cm*5cm proton exchange membrane fuel cell was simulated using a two-dimensional + one-dimensional mathematical model of multiphysics coupling. The main parameters are as follows: the gas channels within the bipolar plates adopt a parallel DC channel mode, with humidification of the gas inlet at a relative humidity of 90%. Both anode and cathode gases flow along the positive x-axis. The anode hydrogen flow rate is 0.0013 mol / s, and the cathode air flow rate is 0.0032 mol / s. The inlet gas temperature is 80°C, and the outlet pressure is 1.7 bar. The voltage is 0.95–0.5 V. The coolant temperature is set to a constant 40°C. Formulas (14)-(16) are not solved in the two-dimensional network model.

[0228] In the two-dimensional network model, Δx = 1 cm, Δy = 1 cm, and there are a total of 25 network units. A two-dimensional + one-dimensional model is established and solved. The two-dimensional + one-dimensional mathematical model of multi-physics coupling in an intrinsic exchange membrane fuel cell can be solved using various numerical calculation methods. In one embodiment, the three-stage Lobatto IIIa algorithm with fourth-order accuracy is used to solve the one-dimensional differential model, and the BICG algorithm is used to solve the two-dimensional network model.

[0229] like Figure 8 As shown, the average current density of the fuel cell under different voltage conditions is obtained through numerical calculation. Figure 9 The figure shows the current density inside the fuel cell when the voltage is 0.7V. Figure 10 As shown, the relative humidity is distributed in the cathode gas channel.

[0230] In addition, the ambient temperature, coolant inlet temperature, and flow rate can be set, and the temperature distribution can be obtained by solving formulas (14)-(16); and / or

[0231] The anode gas is set to flow in the positive x-axis direction, and the cathode gas flows in the negative x-axis direction, i.e., the anode and cathode gases flow in opposite directions; and / or Δx = 1 mm, N x =50, that is, the number of networks in the encrypted x-axis direction is used for numerical simulation; replace the previous steps to simulate the working performance of the 5cm*5cm proton exchange membrane fuel cell again.

[0232] The performance prediction method for proton exchange membrane fuel cells proposed in this invention can analyze the distribution of various physical quantities in the internal components of the membrane electrode assembly (MEA) of a large-size proton exchange membrane fuel cell, reveal the coupling mechanism of heat, electricity, and flow, and provide design ideas for improving the performance of the MEA. Compared with existing technologies, it has higher computational efficiency while maintaining the resolution of the multiphysics field inside the MEA. Furthermore, combined with experimental verification, it can deeply analyze the factors restricting the performance of fuel cells during the development and design of the MEA, quickly verify the effect of parameter optimization design, save design and development costs, and accelerate the iterative design process.

[0233] Next, the performance prediction device for a proton exchange membrane fuel cell according to an embodiment of the present invention is described with reference to the accompanying drawings.

[0234] Figure 11 This is a block diagram of a proton exchange membrane fuel cell performance prediction device according to an embodiment of the present invention.

[0235] like Figure 11 As shown, the performance prediction device 100 for the proton exchange membrane fuel cell includes: a division module 101, a construction module 102, a coupling module 103, and a prediction module 104.

[0236] The system comprises the following modules: A partitioning module 101 divides the target monolithic fuel cell into multiple network units, each including a bipolar plate, gas channels, membrane electrode assembly (MEA), and coolant on both upper and lower sides for heat dissipation. A construction module 102 establishes the relationships between the multiple network units, constructs a two-dimensional network model in a two-dimensional plane based on first preset physical boundary conditions, and establishes a one-dimensional differential model in the MEA of each network unit. A coupling module 103 couples the two-dimensional network model and the one-dimensional differential model to establish a two-dimensional + one-dimensional mathematical model of the proton exchange membrane fuel cell with multi-physics coupling. A prediction module 104 uses the two-dimensional + one-dimensional mathematical model of the proton exchange membrane fuel cell with multi-physics coupling to solve for the distributions of current density, pressure, flow rate, temperature, humidity, and gas concentration to predict the performance of the target monolithic fuel cell.

[0237] Furthermore, in one embodiment of the present invention, the dividing module 101 is specifically used for:

[0238] The two-dimensional plane in which the target monolithic fuel cell is located is set as the xy plane, with the positive x-axis being the direction of anode gas flow and the z-axis being the thickness direction of the target monolithic fuel cell.

[0239] Divide the target monolithic fuel cell into N sections along the y-axis. y N network units, each with a size of Δy, are uniformly divided along the x-axis. x One network element with size Δx is used to obtain multiple network elements, which are numbered (i,j).

[0240] Furthermore, in one embodiment of the present invention, the target monolithic fuel cell is divided into multiple network units, and further includes:

[0241] A layer of boundary virtual point grid cells with the same size as the network cells is set on the outer wall of the target single fuel cell after partitioning, so as to set the governing equations corresponding to the boundary conditions.

[0242] Furthermore, in one embodiment of the present invention, the two-dimensional network model includes:

[0243]

[0244] in, denoted as r, the flow rates of component r in the mixed gas at the center of network cell (i,j) and the center of network cell (i,j+1), respectively. r represents hydrogen and water at the anode and oxygen at the cathode. The mass flux density of gas component r through the airway-membrane electrode interface in network units (i,j) and (i,j+1) are respectively. gc,z α represents the contact area between the airway and the membrane electrode in a single network unit. i,j α i,j-1 These are the pressure loss factors in network units (i,j) and (i,j-1), respectively, p i,j+1 p i,j p i,j-1 Let β represent the pressure at the center of network cell (i,j+1), the center of (i,j), and the center of (i,j-1), respectively. i,j β i,j-1 These are the temperature factors at the center of network cell (i,j) and the center of network cell (i,j-1), respectively. h represents the total mass flux density of the mixed gas passing through the gas channel-membrane electrode interface at the centers of network units (i,j), (i,j+1), and (i,j-1), respectively. gb A is the convective heat transfer coefficient between the airflow and the bipolar plates in the air passage. gb h is the contact area between the airway and the bipolar plate in a single network unit. gc,zC is the convective heat transfer coefficient between the airflow in the gas channel and the boundary of the gas diffusion layer of the membrane electrode. p,k Let K be the specific heat capacity of gaseous component k, where k represents only hydrogen at the anode and only oxygen at the cathode. C represents the mass flux density of gas component k through the gas channel-membrane electrode interface in network units (i,j) and (i,j+1), respectively. p,r Let r be the specific heat capacity of the gas component. C represents the gas temperature at the center of network cell (i,j) and the center of network cell (i,j+1), respectively. p,v The specific heat capacity of water vapor. C represents the mass flux density of water vapor through the gas channel-membrane electrode interface in network units (i,j) and (i,j-1), respectively. p,l The specific heat capacity of liquid water, These represent the mass flux density of liquid water through the gas channel-membrane electrode interface in network units (i,j) and (i,j-1), respectively. These are the boundary temperatures of the membrane electrode gas diffusion layer in network cell (i,j) and (i,j+1), respectively. These represent the temperatures of the bipolar plates in network units (i,j), (i,j+1), and (i,j-1), respectively, and α is the temperature of the bipolar plate. bp k is the thermal diffusivity factor. bp Let A be the thermal conductivity of the bipolar plate. bp,z b represents the contact area between the bipolar plate and the membrane electrode in a single network unit. bp h represents the thickness of the bipolar plate along the z-axis. cool A is the convective heat transfer coefficient between the coolant and the bipolar plate. cool A represents the contact area between the bipolar plate and the coolant. bp,x Let A be the cross-sectional area of ​​the bipolar plate along the x-axis, Δx be the dimension of each network cell along the x-axis, and A be the cross-sectional area of ​​the bipolar plate along the x-axis. bp,y Let be the cross-sectional area of ​​the bipolar plate along the y-axis, and Δy be the dimension of each network cell along the y-axis. Let be the boundary temperature of the bipolar plate in network cell (i,j) at the contact surface with the membrane electrode. Let q represent the temperatures of the coolant in network cell (i,j), network cell (i,j+1), and network cell (i,j-1), respectively. cool C represents the coolant flow rate. p,cool This refers to the specific heat capacity of the coolant.

[0245] Furthermore, in one embodiment of the present invention, the process of solving for the gas pressure in the two-dimensional network model is as follows:

[0246] The expression for constructing the flow velocity in the airway is:

[0247]

[0248] The coordinate axes x', y', and z' are parallel to the x, y, and z axes, respectively. The range of x' is [0, Δx], the range of y' is [-a / 2, a / 2], and the range of z' is [-b / 2, b / 2]. x' Let x' be the gas velocity component along the x' axis, distributed across the cross-section of the airway, μ be the gas dynamic viscosity coefficient, and p be the gas pressure.

[0249] For velocity u x' Integrating over the cross-section of the DC channel and combining this with the perfect gas law pM = ρRT, the molar flow rate q is obtained. n The relationship between (x') and the pressure gradient is as follows:

[0250]

[0251] Where T is the gas temperature, Γ is an intermediate variable, and q n (x) represents the gas molar flow rate;

[0252] The molar flow rate q in the airway of network unit (i,j) n The first transformation is obtained by integrating the relationship between (x') and the pressure gradient:

[0253]

[0254] Where, p i,j p i,j+1 These represent the pressure at the center of network cell (i,j) and the center of network cell (i,j+1), respectively. Let A be the gas temperature at the center of network cell (i,j) and the center of network cell (i,j+1), respectively. gc,z This refers to the contact area between the airway and the membrane electrode in a single network unit. denoted as the total mass flux density of the mixed gas passing through the gas channel-membrane electrode interface at the center of network unit (i,j) and the center of (i,j+1), respectively. Let be the gas molar flow rate at the inlet of network cell (i,j);

[0255] The molar flow rate q in the airway of network unit (i,j-1) n The second transformation is obtained by integrating the equation relating (x') to the pressure gradient:

[0256]

[0257] Where, p i,j-1 The pressure at the center of network cell (i,j-1) is... Let (i,j-1) be the gas temperature at the center of the network cell. Let be the total mass flux density of the mixed gas at the center of network cell (i,j) through the gas channel-membrane electrode interface. Let $\frac{i}{j-1}$ be the gas molar flow rate at the inlet of network cell (i,j-1).

[0258] Subtract the first transformation from the second transformation and substitute the result into the airway flow control equation. The governing equations for airway pressure in network cell (i,j) are obtained as follows:

[0259]

[0260] Where, α i,j α i,j-1 These are the pressure loss factors in network units (i,j) and (i,j-1), respectively, β i,j β i,j-1 These are the temperature factors at the center of network cell (i,j) and the center of network cell (i,j-1), respectively. The total mass flux density of the mixed gas at the center of network cell (i,j+1) through the gas channel-membrane electrode interface.

[0261] Furthermore, in one embodiment of the present invention, the one-dimensional differential model is:

[0262]

[0263] Where, σ p φ is the conductivity of the proton. p Let z be the proton potential, z be the coordinate along the thickness of the film electrode, and s be the coordinate along the thickness of the film electrode. p s e These are the source terms of the anode and cathode catalyst layers, respectively, σ e φ is the conductivity of electrons. e For the electron potential, j an and j ca These represent the electrochemical reaction rates at the anode and cathode, respectively, in seconds. l This represents the volume fraction of liquid water. These represent the anode and cathode exchange current densities, respectively; F is the Faraday constant; R is the gas constant; T is the temperature; and η is the estrus. an and η ca These are the overpotentials at the anode and cathode, respectively. and These are the four transfer coefficients in the Butler-Volmer model, V ocv ω is the open-circuit voltage of the proton exchange membrane fuel cell, ω is the volume fraction of the ion-polymer, and D is the open-circuit voltage. λ Let ρ be the diffusion coefficient of film water. memLet be the density of the ionomer, EW be the equivalent molar mass of the ionomer, λ be the membrane water content, and i be the density of the ionomer. p For proton current density, s ad D is the source term for water adsorption / desorption by ionomers. s Let ρ be the diffusion coefficient of liquid water. l M is the density of liquid water. H2O Let be the molar mass of water, and s be the volume fraction of liquid water. ec For the liquid water source term, c represents the concentration of the mixed gas. Let x be the effective diffusion coefficient of gas component r. r Let s be the volume fraction of gaseous component r, where r represents hydrogen and water vapor at the anode and oxygen and water vapor at the cathode. r For gaseous component source terms, For the equivalent thermal conductivity, s T This is a heat absorption / heat release source term.

[0264] Furthermore, in one embodiment of the present invention, the coupling module 103 is specifically used for:

[0265] The two-dimensional network model is solved based on the preset physical boundary conditions and the second preset boundary conditions and the temperature field to update the gas pressure, humidity, gas components and liquid water volume fraction in the two-dimensional network model, and serve as the first preset boundary conditions required to solve the one-dimensional differential model.

[0266] Based on the control equation and temperature field of the membrane electrode boundary temperature in the first preset boundary condition, a one-dimensional differential model is solved to update the distribution of electric potential, current density, material composition and endothermic / exothermic source terms in the membrane electrode, and obtain the material flux density of each material through the membrane electrode and the gas channel interface, which is used as the second preset boundary condition required to solve the two-dimensional network model.

[0267] The two-dimensional network model is solved based on physical boundary conditions and second preset boundary conditions to update the temperature distribution of the two-dimensional network model and obtain the temperature of the membrane electrode boundary, which serves as the first preset boundary condition required for solving the one-dimensional differential model.

[0268] Based on the control equation and energy equation of the membrane electrode boundary temperature in the first preset boundary condition, the one-dimensional differential model is solved to update the temperature distribution inside the membrane electrode and obtain the heat flux density at the gas channel-membrane electrode interface, which serves as the second preset boundary condition for the two-dimensional network model.

[0269] The above solution process is iterated, continuously updating the current density, electric potential, temperature, gas pressure, humidity, volume fraction of each gas component and liquid water, mass flux density and heat flux density at the gas channel-membrane electrode interface, until the iteration converges, resulting in a two-dimensional + one-dimensional mathematical model of multi-physics coupling of the proton exchange membrane fuel cell.

[0270] It should be noted that the explanation of the above-described embodiment of the performance prediction method for proton exchange membrane fuel cells also applies to the performance prediction device for proton exchange membrane fuel cells in this embodiment, and will not be repeated here.

[0271] The performance prediction device for proton exchange membrane fuel cells proposed in this invention can analyze the distribution of various physical quantities in the internal components of the membrane electrode assembly (MEA) of a large-size proton exchange membrane fuel cell, reveal the coupling mechanism of heat, electricity, and flow, and provide design ideas for improving the performance of the MEA. Compared with existing technologies, it has higher computational efficiency while maintaining the resolution of the multiphysics field inside the MEA. Furthermore, combined with experimental verification, it can deeply analyze the factors restricting the performance of fuel cells during the development and design of the MEA, quickly verify the effect of parameter optimization design, save design and development costs, and accelerate the iterative design process.

[0272] Figure 12 This is a schematic diagram of an electronic device provided in an embodiment of the present invention. The electronic device may include:

[0273] The memory 1201, the processor 1202, and the computer program stored on the memory 1201 and executable on the processor 1202.

[0274] When the processor 1202 executes the program, it implements the proton exchange membrane fuel cell performance prediction method provided in the above embodiments.

[0275] Furthermore, electronic devices also include:

[0276] Communication interface 1203 is used for communication between memory 1201 and processor 1202.

[0277] The memory 1201 is used to store computer programs that can run on the processor 1202.

[0278] The memory 1201 may include high-speed RAM memory, and may also include non-volatile memory, such as at least one disk storage device.

[0279] If the memory 1201, processor 1202, and communication interface 1203 are implemented independently, then the communication interface 1203, memory 1201, and processor 1202 can be interconnected via a bus to complete communication between them. The bus can be an Industry Standard Architecture (ISA) bus, a Peripheral Component Interconnect (PCI) bus, or an Extended Industry Standard Architecture (EISA) bus, etc. Buses can be divided into address buses, data buses, control buses, etc. For ease of representation, Figure 12 The bus is represented by a single thick line, but this does not mean that there is only one bus or one type of bus.

[0280] Optionally, in a specific implementation, if the memory 1201, processor 1202, and communication interface 1203 are integrated on a single chip, then the memory 1201, processor 1202, and communication interface 1203 can communicate with each other through an internal interface.

[0281] The processor 1202 may be a central processing unit (CPU), an application specific integrated circuit (ASIC), or one or more integrated circuits configured to implement embodiments of the present invention.

[0282] This invention also provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the above-described method for predicting the performance of a proton exchange membrane fuel cell.

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

[0284] Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of indicated technical features. Thus, a feature defined as "first" or "second" may explicitly or implicitly include at least one of that feature. In the description of this invention, "N" means at least two, such as two, three, etc., unless otherwise explicitly specified.

[0285] Any process or method description in the flowchart or otherwise herein can be understood as representing a module, segment, or portion of code comprising one or N executable instructions for implementing custom logic functions or processes, and the scope of preferred embodiments of the invention includes additional implementations in which functions may be performed not in the order shown or discussed, including substantially simultaneously or in reverse order depending on the functions involved, as will be understood by those skilled in the art to which embodiments of the invention pertain.

[0286] The logic and / or steps represented in the flowchart or otherwise described herein, for example, can be considered as a sequenced list of executable instructions for implementing logical functions, and can be embodied in any computer-readable medium for use by, or in conjunction with, an instruction execution system, apparatus, or device (such as a computer-based system, a processor-included system, or other system that can fetch and execute instructions from, an instruction execution system, apparatus, or device). For the purposes of this specification, "computer-readable medium" can be any means that can contain, store, communicate, propagate, or transmit programs for use by, or in conjunction with, an instruction execution system, apparatus, or device. More specific examples (a non-exhaustive list) of computer-readable media include: an electrical connection having one or more wires (electronic device), a portable computer disk drive (magnetic device), random access memory (RAM), read-only memory (ROM), erasable and editable read-only memory (EPROM or flash memory), fiber optic devices, and portable optical disc read-only memory (CDROM). Alternatively, the computer-readable medium may be paper or other suitable media on which the program can be printed, since the program can be obtained electronically by optically scanning the paper or other medium, followed by editing, interpreting, or otherwise processing as necessary, and then stored in a computer memory.

[0287] It should be understood that various parts of the present invention can be implemented in hardware, software, firmware, or a combination thereof. In the above embodiments, the N steps or methods can be implemented in software or firmware stored in memory and executed by a suitable instruction execution system. If implemented in hardware, as in another embodiment, it can be implemented using any one or a combination of the following techniques known in the art: discrete logic circuits having logic gates for implementing logical functions on data signals, application-specific integrated circuits (ASICs) having suitable combinational logic gates, programmable gate arrays (PGAs), field-programmable gate arrays (FPGAs), etc.

[0288] Those skilled in the art will understand that all or part of the steps of the methods in the above embodiments can be implemented by a program instructing related hardware. The program can be stored in a computer-readable storage medium, and when executed, the program includes one or a combination of the steps of the method embodiments.

[0289] Furthermore, the functional units in the various embodiments of the present invention can be integrated into a processing module, or each unit can exist physically separately, or two or more units can be integrated into a module. The integrated module can be implemented in hardware or as a software functional module. If the integrated module is implemented as a software functional module and sold or used as an independent product, it can also be stored in a computer-readable storage medium.

[0290] The storage medium mentioned above can be a read-only memory, a disk, or an optical disk, etc. Although embodiments of the present invention have been shown and described above, it is to be understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those skilled in the art can make changes, modifications, substitutions, and variations to the above embodiments within the scope of the present invention.

Claims

1. A method for predicting the performance of a proton exchange membrane fuel cell, characterized in that, Includes the following steps: The target monolithic fuel cell is divided into multiple network units. Each network unit includes a bipolar plate, a gas channel, a membrane electrode assembly, and coolant on the upper and lower sides for heat dissipation. The division of the target monolithic fuel cell into multiple network units includes: Let the two-dimensional plane containing the target monolithic fuel cell be... flat, The positive direction of the axis is the direction of anode gas flow. The axis is the thickness direction of the target monolithic fuel cell; The target monolithic fuel cell is in Division along the axis Each network unit has a size of [number]. ,exist Uniformly divided along the axial direction Each network unit has a size of [number]. The plurality of network units are obtained, and the unit numbers are as follows: ; Establish the relationship between the multiple network units, construct a two-dimensional network model in the two-dimensional plane direction in combination with the preset physical boundary conditions, and establish a one-dimensional differential model in the membrane electrode of each network unit. The two-dimensional network model and the one-dimensional differential model are coupled to establish a two-dimensional + one-dimensional mathematical model of multi-physics coupling for a proton exchange membrane fuel cell, specifically including: The two-dimensional network model is solved based on the preset physical boundary conditions and the second preset boundary conditions and the temperature field to update the gas pressure, humidity, gas components and liquid water volume fraction in the two-dimensional network model, and to serve as the first preset boundary conditions required to solve the one-dimensional differential model. Based on the control equation and temperature field of the membrane electrode boundary temperature in the first preset boundary condition, the one-dimensional differential model is solved to update the distribution of electric potential, current density, each material component and endothermic / exothermic source terms in the membrane electrode, and to obtain the material flux density of each substance passing through the membrane electrode and the airway interface, which is used as the second preset boundary condition required to solve the two-dimensional network model. The two-dimensional network model is solved based on the physical boundary conditions and the second preset boundary conditions to update the temperature distribution of the two-dimensional network model and obtain the temperature of the membrane electrode boundary, which serves as the first preset boundary condition required for solving the one-dimensional differential model. Based on the control equation and energy equation of the membrane electrode boundary temperature in the first preset boundary condition, the heat absorption / heat release source term of the energy equation is used to solve the one-dimensional differential model to update the temperature distribution in the membrane electrode and obtain the heat flux density at the airway-membrane electrode interface, which serves as the second preset boundary condition for the two-dimensional network model. The above solution process is iterated, continuously updating the current density, the electric potential, the temperature, the gas pressure, the humidity, the volume fraction of each gas component and liquid water, the mass flux density and heat flux density of the gas channel-membrane electrode interface, until the iteration converges, and a two-dimensional + one-dimensional mathematical model of the multi-physics coupling of the proton exchange membrane fuel cell is obtained. The distributions of current density, temperature, gas pressure, flow rate, humidity, and concentration are solved using a two-dimensional + one-dimensional mathematical model of multi-physics coupling in the proton exchange membrane fuel cell to predict the performance of the target monolithic fuel cell.

2. The method for predicting the performance of a proton exchange membrane fuel cell according to claim 1, characterized in that, The method of dividing the target monolithic fuel cell into multiple network units further includes: A layer of boundary virtual point mesh elements with the same size as the network elements is set on the outer wall of the target single fuel cell after partitioning, so as to set the equations corresponding to the boundary conditions.

3. The method for predicting the performance of a proton exchange membrane fuel cell according to claim 1, characterized in that, The two-dimensional network model is as follows: in, , Network units central area Components in the central gas mixture Traffic, The anode represents hydrogen and water, and the cathode represents oxygen. Network units , Gas components passing through the gas passage-membrane electrode interface The mass flux density, This refers to the contact area between the airway and the membrane electrode in a single network unit. , Network units , Pressure loss factor in , , Network units central area central area Pressure at the center , Network units central area Temperature factor at the center , , Network units central area central area The total mass flux density of the mixed gas at the center through the gas channel-membrane electrode interface. This is the convective heat transfer coefficient between the airflow and the bipolar plates in the air passage. This refers to the contact area between the airway and the bipolar plate in a single network unit. This is the convective heat transfer coefficient between the airflow in the gas channel and the boundary of the gas diffusion layer of the membrane electrode. gaseous components k Specific heat capacity, The anode represents only hydrogen gas, and the cathode represents only oxygen gas. , Network units Network Unit Middle gas components k Mass flux density through the airway-membrane electrode interface gaseous components Specific heat capacity, , Network units central area The gas temperature at the center The specific heat capacity of water vapor. , Network units , The mass flux density of water vapor through the gas channel-membrane electrode interface. The specific heat capacity of liquid water, , Network units , The mass flux density of liquid water through the gas channel-membrane electrode interface. , Network units middle, The boundary temperature of the gas diffusion layer in the middle film electrode. , , Network units , , The temperature of the bipolar plate For thermal diffusivity, The thermal conductivity of the bipolar plate is... This refers to the contact area between the bipolar plate and the membrane electrode in a single network unit. For bipolar plates in Thickness in the axial direction, The convective heat transfer coefficient between the coolant and the bipolar plate. This represents the contact area between the bipolar plate and the coolant. For bipolar plates in Cross-sectional area along the axial direction, For each network unit Dimensions in the axial direction For bipolar plates in Cross-sectional area along the axial direction, For each network unit Dimensions in the axial direction For network units The boundary temperature of the bipolar plate at the contact surface with the membrane electrode. , , Network units middle, middle, The temperature of the coolant This refers to the coolant flow rate. This refers to the specific heat capacity of the coolant.

4. The method for predicting the performance of a proton exchange membrane fuel cell according to claim 3, characterized in that, The process of solving for the gas pressure in the two-dimensional network model is as follows: The expression for constructing the flow velocity in the airway is: Among them, coordinate axes , , The shafts are respectively with x , y , z Axis parallel, coordinate The range of values ​​for is ,coordinate The range of values ​​for is ,coordinate The range of values ​​for is , for The axial gas velocity component is distributed across the cross-section of the airway. The gas dynamic viscosity coefficient, p This refers to gas pressure. speed Integrate over the cross-section of the DC channel and combine with the complete gas law. To obtain the molar flow rate The relationship between pressure gradient and pressure gradient: in, For gas temperature, As an intermediate variable, This represents the gas molar flow rate; In network unit Molar flow rate in the airway The first transformation is obtained by integrating the equation relating the pressure gradient to the pressure gradient: in, , Network units central area Pressure at the center , Network units central area The gas temperature at the center This refers to the contact area between the airway and the membrane electrode in a single network unit. , Network units central area The total mass flux density of the mixed gas at the center through the gas channel-membrane electrode interface. For network units Gas molar flow rate at the inlet; In network unit Molar flow rate in the airway The second transformation is obtained by integrating the equation relating the pressure gradient to the pressure gradient: in, For network units Pressure at the center For network units The gas temperature at the center For network units The total mass flux density of the mixed gas at the center through the gas channel-membrane electrode interface. For network units Gas molar flow rate at the inlet; Subtract the first transformation from the second transformation and substitute the result into the airway flow control equation. To obtain network units The governing equation for middle airway pressure: in, , Network units , Pressure loss factor in , Network units central area Temperature factor at the center.

5. The method for predicting the performance of a proton exchange membrane fuel cell according to claim 1, characterized in that, The one-dimensional differential model is as follows: in, The conductivity of the proton. For proton potential, The coordinates are along the thickness direction of the membrane electrode. , These are the source terms of the anode and cathode catalyst layers, respectively. The conductivity of electrons, For electron potential, and These represent the electrochemical reaction rates at the anode and cathode, respectively. This represents the volume fraction of liquid water. , These represent the exchange current densities at the anode and cathode, respectively. F It is Faraday's constant. R The gas constant is... T For temperature, and These are the overpotentials at the anode and cathode, respectively. , , and These are the four transfer coefficients in the Butler-Volmer model. This refers to the open-circuit voltage of a proton exchange membrane fuel cell. This represents the volume fraction of the ionic polymer. Let be the diffusion coefficient of film water. The density of the ionomer. This represents the equivalent molar mass of the ionomer. The content of membrane water, The proton current density, This is the source term for water adsorption / desorption by ion polymers. Let be the diffusion coefficient of liquid water. The density of liquid water, The molar mass of water, This represents the volume fraction of liquid water. For liquid water source items, The concentration of the mixed gas. gaseous components The effective diffusion coefficient, gaseous components volume fraction, The components are gaseous, representing hydrogen and water vapor at the anode and oxygen and water vapor at the cathode. For gaseous component source terms, The equivalent thermal conductivity is This is a heat absorption / heat release source term.

6. A device for predicting the performance of a proton exchange membrane fuel cell, characterized in that, include: The partitioning module is used to divide the target monolithic fuel cell into multiple network units. Each network unit includes bipolar plates, gas channels, membrane electrode assemblies, and coolant on the upper and lower sides for heat dissipation. Specifically, it includes: Let the two-dimensional plane containing the target monolithic fuel cell be... flat, The positive direction of the axis is the direction of anode gas flow. The axis represents the thickness direction of the target monolithic fuel cell; The target monolithic fuel cell in Division along the axis Each network unit has a size of [number]. ,exist Uniformly divided along the axial direction Each network unit has a size of [number]. Multiple network elements were obtained, with element numbers as follows: ; The construction module is used to establish the relationship between the multiple network units, construct a two-dimensional network model in the two-dimensional plane direction in combination with preset physical boundary conditions, and establish a one-dimensional differential model in the membrane electrode of each network unit. The coupling module, used to couple the two-dimensional network model and the one-dimensional differential model, establishes a two-dimensional + one-dimensional mathematical model of multi-physics coupling in a proton exchange membrane fuel cell, specifically including: The two-dimensional network model is solved based on the preset physical boundary conditions and the second preset boundary conditions and the temperature field to update the gas pressure, humidity, gas components and liquid water volume fraction in the two-dimensional network model, and serve as the first preset boundary conditions required to solve the one-dimensional differential model. Based on the control equation and temperature field of the membrane electrode boundary temperature in the first preset boundary condition, a one-dimensional differential model is solved to update the distribution of electric potential, current density, material composition and endothermic / exothermic source terms in the membrane electrode, and obtain the material flux density of each material through the membrane electrode and the gas channel interface, which is used as the second preset boundary condition required to solve the two-dimensional network model. The two-dimensional network model is solved based on physical boundary conditions and second preset boundary conditions to update the temperature distribution of the two-dimensional network model and obtain the temperature of the membrane electrode boundary, which serves as the first preset boundary condition required for solving the one-dimensional differential model. Based on the control equation and energy equation of the membrane electrode boundary temperature in the first preset boundary condition, the one-dimensional differential model is solved to update the temperature distribution inside the membrane electrode and obtain the heat flux density at the gas channel-membrane electrode interface, which serves as the second preset boundary condition for the two-dimensional network model. The above solution process is iterated, continuously updating the current density, electric potential, temperature, gas pressure, humidity, volume fraction of each gas component and liquid water, mass flux density and heat flux density at the gas channel-membrane electrode interface, until the iteration converges, and a two-dimensional + one-dimensional mathematical model of multi-physics coupling of proton exchange membrane fuel cell is obtained. The prediction module is used to solve the distribution of current density, pressure, flow rate, temperature, humidity and gas concentration using a two-dimensional + one-dimensional mathematical model of multi-physics coupling of the proton exchange membrane fuel cell, so as to predict the performance of the target monolithic fuel cell.

7. An electronic device, characterized in that, include: The memory, the processor, and the computer program stored in the memory and executable on the processor, the processor executing the program to implement the proton exchange membrane fuel cell performance prediction method as described in any one of claims 1-5.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that, The program is executed by a processor to implement the proton exchange membrane fuel cell performance prediction method as described in any one of claims 1-5.

Images