Proton exchange membrane electrolytic cell electric-water-hot gas multi-physical field performance prediction software architecture and construction method

Abstract

The invention discloses a proton exchange membrane electrolytic cell electric-water-hot gas multi-physical field performance prediction software framework and a construction method. The software framework comprises a user-defined layer opened to a user and a bottom layer framework invisible to the user; the function of the user-defined layer is realized through a user-defined function; the bottom-layer architecture comprises a storage-oriented data structure, a CFD universal solver input and a solving-oriented data structure; the storage-oriented data structure is an underlying data packet and is used for storing various data; the CFD universal solver input is used for receiving data packet data in the storage-oriented data structure and writing an xml file configuration to the solving-oriented data structure; the solution-oriented data structure is used for solving xml file configuration to obtain different physical fields in the proton exchange membrane electrolytic tank. The method has the characteristics that the degree of freedom of user operation is high, the calculation stability is higher, and the solved physical field better conforms to the reality.

Description

Technical Field

[0001] The present invention relates to the technical field of proton exchange membrane electrolyzers, and particularly relates to a software architecture and construction method for predicting the multi-physical field performance of proton exchange membrane electrolyzers in terms of electricity, water, heat, and gas. Background Technique

[0002] Proton exchange membrane electrolyzers have the advantages of being structurally compact, having high hydrogen production purity, and high electrolysis efficiency. Their main components include: a catalytic layer, a microporous layer (which may not be included), a gas diffusion layer, and the anode and cathode plates. When the electrolytic cell operates, the water introduced from the anode inlet is evenly distributed under the guidance of the flow channel, and then reaches the surface of the catalytic layer through the gas diffusion layer and the microporous layer, where an oxidation reaction occurs to generate oxygen. At the same time, protons are transported across the membrane to the surface of the cathode catalytic layer, where a reduction reaction occurs to generate hydrogen. The equations of the above reactions are shown as follows:

[0003] Anode:

[0004] Cathode: 2H + +2e - →H2

[0005] Total reaction:

[0006] During the operation of a proton exchange membrane electrolyzer, it involves multiphase, multicomponent, and complex multi-physical field coupled transport and electrochemical reactions. Accurately predicting the multi-physical fields inside the electrolyzer is of great significance for guiding its design and performance improvement. The single cell inside a proton exchange membrane electrolyzer involves multi-physical fields, multi-components, and multiphase problems, and the coupling relationships between different physical fields are complex. The gas-liquid two-phase flow also increases the difficulty of solving. With the continuous advancement of the commercialization process of proton exchange membrane electrolyzers, there is an urgent need for a software for predicting the multi-physical field performance of electrolyzers with strong generality and stable calculation.

[0007] Although some existing commercial software has a simulation module for proton exchange membrane electrolyzers, there are still many problems, such as:

[0008] (1) The simulation module of the proton exchange membrane electrolyzer in the existing commercial software is highly encapsulated, and users cannot view or modify it from the outside;

[0009] (2) Using a single-component system results in over-solving and over-determination problems;

[0010] (3) Solving all control equations using a full-field solution instead of solving in the control region where the equation truly corresponds, resulting in abnormal physical fields in some cases;

[0011] (4) The electrolyzer module of some commercial software adopts a single-phase assumption, which has a large deviation from the actual situation;

[0012] (5) Although some commercial software overcomes the single-phase assumption, only one main phase and one secondary phase can be set in the multiphase flow simulation. However, in actual situations, the volume fraction of the anode water in the electrolytic cell is dominant, and the volume fraction of the cathode gas is dominant. The secondary phase usually represents a phase with a smaller volume fraction. It is unreasonable to set the same main and secondary phases at both electrodes, which is not conducive to calculation stability.

[0013] (6) The multiphase model of existing commercial software can only be set in the entire fluid domain. However, under some working conditions, there is almost no liquid water at the cathode of the electrolytic cell. In this case, it is unreasonable to still use a multiphase model for the cathode.

[0014] In the study by Zhuang et al. (Document No.: INTERNATIONAL JOURNAL OF HYDROGEN EN ERGY 49 (2024) 337-352), a model of a single direct current channel of a proton exchange membrane electrolyzer was established and solved using COMSOL Multiphysics software, and the established model adopted a single-phase assumption. It can also be learned from the Fuel Cell & Electrolyzer Module User's Guide of COMSOL Multiphysics 6.2 that the proton exchange membrane electrolyzer module that comes with COMSOL adopts a single-phase model. However, in the actual operation of the electrolyzer, the anode consumes liquid water to produce oxygen, which is a two-phase flow. The single-phase model cannot reflect the actual situation inside the electrolyzer.

[0015] In summary, the proton exchange membrane electrolyzer module of existing software adopts a single-phase assumption. Although some software has overcome the single-phase assumption, the multiphase model can only be set in the entire fluid domain, and there can only be one main phase and secondary phase setting in the multiphase flow simulation, which is obviously inconsistent with the actual physical process. Summary of the invention

[0016] In order to overcome the defects of the above-mentioned prior art, the present invention provides a software framework and construction method for predicting the multi-physical field performance of electricity, water, heat and gas in a proton exchange membrane electrolyzer. The user has more flexible operation space, can set different main phases and secondary phases at the cathode and anode, and overcomes the problems of over-solving and abnormalities of some physical fields. The method has the characteristics of high user operation freedom, higher calculation stability, and the physical field obtained by the solution is more consistent with the actual situation.

[0017] In order to achieve the above object, the technical solution adopted by the present invention is:

[0018] A software framework for predicting the multi-physical performance of electricity, water, heat and gas in a proton exchange membrane electrolyzer, including a user-friendly custom layer and an underlying architecture that is invisible to the user;

[0019] The function of the custom layer is implemented by a custom function;

[0020] The underlying architecture includes a storage-oriented data structure, a CFD general solver input, and a solution-oriented data structure;

[0021] The storage-oriented data structure is the underlying data packet for storing data generated during various solution processes;

[0022] The CFD general solver input is used to receive data from the storage-oriented data structure and write an xml file for configuration to the solution-oriented data structure;

[0023] The solution-oriented data structure is used to solve the xml file configuration to obtain different physical fields inside the proton exchange membrane electrolyzer.

[0024] The custom function includes five types: custom source term, custom diffusion coefficient, custom profile function, solution strategy interfaces before and after iteration, and custom initialization function. Among them, the solution strategy interfaces before and after iteration implement the information transfer between the custom layer and the underlying architecture; the custom profile function is located in the user-defined layer. Inside the custom profile function, the units on the boundary (including but not limited to the flow channel inlet, outlet, end faces of bipolar plates, etc.) can be traversed, and the face center values of the units can be assigned to set the boundary conditions. Then, the defined boundary conditions are configured into the solver by registering the custom profile function; the solver is located in the underlying architecture, including a custom scalar / vector equation solver and also a CFD general solver for solving the temperature field and flow field.

[0025] The custom initialization function initializes the physical field before solving the control equation.

[0026] The solution-oriented data structure includes a flow and heat transfer module, a cathode component module, an anode component module, a cathode electron potential module, an anode electron potential module, a proton potential module, a membrane water content module, a liquid water pressure module, and a multiphase module; the above modules all correspond to different control equations, and there is coupling between different modules; the solution results of the equations in the nine modules correspond to different physical fields inside the proton exchange membrane electrolyzer.

[0027] The solution of multiple physical fields inside the proton exchange membrane electrolyzer includes the following steps:

[0028] Step (1), first, initialize. Assign initial values to some custom scalar / vector fields through the custom function in the custom layer, and also assign initial values to the flow field temperature field solver and related variables;

[0029] Step (2), update the physical properties, source terms, and boundaries in the pre-iteration solution strategy interface of the custom layer;

[0030] Step (3): Solve the continuity equation and momentum equation. There is a source term related to the slip velocity in the momentum equation. First, input the average velocity of the mixed phase, and record the relative velocity at the previous iteration step within the slip velocity model. Use a series of coupling relationships to calculate the new relative velocity. When the iteration ends, calculate the source term related to the slip velocity and add it to the momentum source term;

[0031] Step (4):

[0032] Solve the energy equation;

[0033] Solve the component equations, with hydrogen being solved at the cathode and oxygen being solved at the anode;

[0034] Solve six custom scalar equations for the anode potential, cathode potential, proton potential, membrane water content, liquid water pressure, and secondary phase volume fraction. The membrane water content equation may not need to be solved;

[0035] Step (5): In the iteration post - solution strategy interface of the custom layer, calculate and update the physical quantities to be output, and calculate the residuals;

[0036] Step (6): If convergence is achieved, the calculation ends and the results are output. If not, return to Step 2.

[0037] The specific method steps for the initial value assignment in Step (1) are as follows: In the defined initialization function in the custom layer, traverse the areas that need to be initialized. In each area that needs to be initialized, traverse the cells in the area to assign the specified initial values.

[0038] In Step (2), the specific method steps for updating physical properties, source terms, and boundaries are as follows: In the iteration pre - solution strategy interface function of the custom layer, traverse the areas where unit information such as physical properties and source terms need to be updated. In each area that needs to be updated, traverse the cells in the area. Since the information to be updated is stored in the form of custom scalar / vector fields, updating the corresponding custom scalar / vector fields of the cells can achieve the update of unit information such as physical properties and source terms. For the update of boundaries, also in the iteration pre - solution strategy interface function of the custom layer, traverse the boundaries that need to be updated, traverse the faces on each boundary that needs to be updated, and assign the values to be updated to the custom scalar / vector fields of the unit faces.

[0039] In Step (3), the continuity equation is as follows:

[0040]

[0041] Where: ρ m is the density of the mixed phase; S m is the mass source term.

[0042] The density and velocity of the mixture phase can be expressed as follows:

[0043]

[0044] The momentum equation is as follows:

[0045]

[0046] Where: S u is the momentum source term.

[0047] The coupling relationships of the slip velocity, relative velocity, and mixture phase velocity are as follows:

[0048]

[0049]

[0050] Where: U sec,p is the relative velocity of the secondary phase and the primary phase (U p ), and is defined as follows:

[0051] U sec,p = U sec - U p

[0052] τ sec is the relaxation time; d sec is the particle size of the secondary phase bubbles / droplets; μ p is the dynamic viscosity of the primary phase; f drag is the drag function; Re is the Reynolds number; is the acceleration; is the gravitational acceleration.

[0053] In the step (4), the solution methods for the energy equation and a series of scalar and vector equations are similar. The above equations can be expressed in the following general form. From left to right, the equations are the transient term, convection term, diffusion term, and source term:

[0054]

[0055] When solving the diffusion equation in the above general form, each term needs to be discretized. The discretized form of the transient term is as follows:

[0056]

[0057] It should be noted that the above discretization format is the first-order forward difference;

[0058] The discretized form of the convection term is as follows:

[0059]

[0060] The explicit discrete form of the diffusion term is as follows:

[0061]

[0062] The discrete form of the source term is as follows:

[0063] s φ (φ) = S p φ + S c

[0064]

[0065] Rearranging the above discrete equations gives the following algebraic equation system:

[0066] (A T + A C - A D - A S )·[Φ] = (b T + b C - b D - b S )

[0067] where: the subscripts T, C, D, and S represent the unsteady term, the convective term, the diffusion term, and the source term respectively;

[0068] In step (5), the specific method steps for calculating and updating the physical quantities to be output are as follows. In the iterative post-solving strategy interface in the custom layer, for the cell center values that need to be calculated and updated, traverse the areas that need to be updated or calculated. In each area that needs to be updated or calculated, traverse the cells in the area. Since the information that needs to be updated or calculated is stored in the form of a custom scalar / vector field, calculating or updating the custom scalar / vector field corresponding to the cell can achieve the calculation and update of the physical quantities to be output. For the cell face values that need to be calculated and updated, traverse the internal or external boundary faces that need to be updated or calculated. On each boundary face that needs to be updated or calculated, traverse the cell faces on the boundary. Since the information that needs to be updated or calculated is stored in the form of a custom scalar / vector field, calculating or updating the custom scalar / vector field corresponding to the cell face can achieve the calculation and update of the physical quantities to be output. There is a specified function in the underlying architecture to calculate the residual after each iteration.

[0069] In step (6), the method for determining the end of the calculation is that if the number of iterations reaches the maximum value set by the user, or the residual is less than the set value for terminating the calculation, then the calculation ends.

[0070] The cathode component module, anode component module, anode electron potential module, cathode electron potential module, proton potential module, membrane water content module, liquid water pressure module, and multiphase module are all solved by a custom scalar / vector equation solver;

[0071] All custom scalar / vector equations and their solution regions are stored in the form of a two-dimensional dynamic array. One dimension of the array stores the labels of the equations to be solved, and the other dimension stores the labels of the solution regions corresponding to these equations. By adding the label of the region where a certain equation needs to be solved into the two-dimensional dynamic array, the solution of this governing equation in this region is achieved, thereby enabling the setting of corresponding solution regions for each custom scalar / vector equation.

[0072] The general form of the custom scalar / vector equation is as follows. From left to right, the equation consists of a transient term, a convective term, a diffusion term, and a source term:

[0073]

[0074] Among them, for the cathode component module and the anode component module, the software also defines the component data packets for the cathode and anode respectively. The component data packet consists of two custom scalar / vector data packets, corresponding to the component systems of the cathode and anode respectively.

[0075] The software is used to achieve the adaptive primary phase and secondary phase according to the volume fractions of the gas phase and liquid phase in a connected domain

[0076] The secondary phase volume fraction equation is as follows. From left to right, it consists of a transient term, a convective term, a source term related to the slip velocity, and a gas-liquid phase change source term:

[0077]

[0078] Among them: α sec is the secondary phase volume fraction; ρ sec is the secondary phase density; U m is the mixed phase velocity; S v-l is the gas-liquid phase change source term; U dr,sec is the secondary phase slip velocity, that is, the velocity of the secondary phase relative to the mixed phase, which can be expressed by the following formula:

[0079] U dr,sec = U sec - U m

[0080] When the phase fraction crosses 0.5, the primary phase and secondary phase are reversed, and the governing equation is automatically modified.

[0081] The specific method is as follows:

[0082] First, determine the high and low of the gas-phase and liquid-phase volume fractions in the unit. If the liquid-phase volume fraction is greater than the gas phase, the gas phase should be regarded as the secondary phase. When setting the secondary-phase volume fraction equation in the user-defined scalar / vector equation solver, the convective term uses the gas density, and the source term related to the slip velocity in the secondary-phase volume fraction equation also uses the gas slip velocity. Then, solve the gas-phase volume fraction equation. Otherwise, solve the liquid-phase volume fraction equation;

[0083] The above coupling relationship is as follows:

[0084]

[0085] where: U sec,p is the relative velocity of the secondary phase and the primary phase (U p ), and is defined as follows:

[0086] U sec,p = U sec - U p τ sec is the relaxation time; d sec is the particle size of the secondary-phase bubbles / droplets; μ p is the dynamic viscosity of the primary phase; f drag is the drag function; Re is the Reynolds number; is the acceleration; is the gravitational acceleration.

[0087] A construction method for a software for electro-hydrothermal management analysis and design of a proton exchange membrane electrolytic cell, comprising the following steps;

[0088] Step 1. Construct a fluid flow and heat transfer module:

[0089] Set the solvers for the continuity equation, momentum equation, and energy equation, and add mass source terms, momentum source terms, and energy source terms;

[0090] Step 2. Construct cathode and anode component modules:

[0091] Set the cathode and anode component solvers. The cathode components include hydrogen and water vapor. Add the source terms and diffusion coefficients of the cathode and anode component equations. The anode components include oxygen, and whether to include water vapor can be selected;

[0092] Step 3. Construct a multiphase module:

[0093] Set the user-defined scalar solver for the secondary-phase volume fraction equation, and add the source term of the secondary-phase volume fraction equation;

[0094] Step 4. Construct cathode and anode electron potential modules:

[0095] Set the user-defined scalar solvers for the cathode and anode potential equations, and add the source terms and diffusion coefficients of the cathode and anode potential equations;

[0096] Step 5. Construct the proton potential module:

[0097] Set the custom scalar solver for the proton potential equation, and add the source term and diffusion coefficient of the proton potential equation;

[0098] Step 6. Construct the liquid water pressure module:

[0099] Set the custom scalar solver for the liquid water pressure equation, and add the source term and diffusion coefficient of the liquid water pressure equation;

[0100] Step 7. Construct the membrane water content module:

[0101] Select whether to solve the membrane water content equation according to the working conditions. If not, specify the membrane water content. If so, set the custom scalar solver for the membrane water content equation, and add the source term and diffusion coefficient of the membrane water content equation;

[0102] Step 8. Add coupling relationships:

[0103] Add the custom functions involved in the above modules, including but not limited to a series of custom functions such as the electrochemical reaction rates of the cathode and anode, the effective permeability of the porous medium, the water phase change source term, the open circuit voltage, and the membrane conductivity.

[0104] Advantages of the present invention:

[0105] The present invention fills the blank in the development of proton exchange membrane electrolyzer simulation software in China and provides a basic framework for the development of related software.

[0106] The present invention develops a user-defined layer for the characteristics of the electrolyzer. Since the software architecture proposed by the present invention includes the interface between the user-defined layer and the underlying layer, and the connection between the user-defined layer and the underlying layer has been configured and improved in advance through the function interface, all the contents related to the computational fluid dynamics algorithm are integrated in the underlying layer and encapsulated for the user, and there is no need to change under different working conditions. The user only needs to modify the source term, boundary conditions, diffusion coefficient, coupling relationship, initialization, and the content required before and after iteration in the user-defined layer according to the research needs. Therefore, the user can perform personalized customization according to the research needs without understanding the software underlying architecture and complex computational fluid dynamics algorithms, which greatly reduces the user's learning cost.

[0107] In step 2 of the construction method of the present invention, by defining different component data packets for the cathode and anode respectively, different component systems are implemented for the cathode and anode, thus overcoming the over-solving and over-determination problems caused by the single-component system of the proton exchange membrane electrolyzer module of foreign commercial software.

[0108] Existing foreign commercial software solves all control equations using a full-field solution instead of solving them in the control region where the equations truly apply, resulting in anomalies in some physical fields. The present invention has developed a custom scalar / vector equation module that can achieve solution only in the region corresponding to the equations, thus avoiding this problem.

[0109] The present invention has made significant improvements in dealing with multiphase problems in electrolytic cells compared to existing foreign commercial software. In the multiphase module of the present invention, the secondary phase volume fraction equation is solved by a custom scalar solver, and any region can be set as the solution region of the custom scalar / vector solver. Therefore, the present invention can apply the multiphase model only in the specified fluid domain, and different primary and secondary phases can be set at the cathode and anode. In a connected domain, the primary and secondary phases can be adaptively adjusted according to the volume fractions of the gas phase and the liquid phase, improving the stability of the calculation. BRIEF DESCRIPTION OF THE DRAWINGS

[0110] Figure 1 It is a schematic diagram of the software framework structure proposed by the present invention.

[0111] Figure 2 It is a flow chart for solving the adaptive primary and secondary phases

[0112] Figure 3 It is a software solution flow chart proposed by the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0113] The present invention will be further described in detail below with reference to the accompanying drawings.

[0114] The present invention provides a proton exchange membrane electrolytic cell electro-hydro-thermal multi-physical field performance prediction software framework and construction method with clear structure and convenient use, filling the gap in the development of domestic related software. Compared with foreign corresponding software modules, users have more flexible operation space, can set different primary and secondary phases at the cathode and anode, and overcome the problems of over-solving and anomalies in some physical fields.

[0115] Figure 1 It is a schematic diagram of the software framework structure proposed by the present invention. As shown in the figure, the software mainly includes nine modules, namely: flow and heat transfer module, cathode component module, anode component module, anode electron potential module, cathode electron potential module, proton potential module, membrane water content module, liquid water pressure module, and multiphase module.

[0116] The above modules correspond to different control equations and are coupled with each other: when solving the fluid flow and heat transfer module to obtain the temperature field, velocity field, and pressure field, the calculation of the mass source term and momentum source term in solving this module requires knowledge of the mass of the generated hydrogen, oxygen, and consumed liquid water. These physical fields need to be calculated through the reaction rates in the anode electron potential module and cathode electron potential module, and the phase fraction in the multiphase module; the calculation of the energy source term requires knowledge of the anode and cathode electron potentials, reaction rates, proton potential, anode and cathode overpotentials, and the mass of water undergoing gas-liquid phase change. These physical fields are obtained through the anode electron potential module, cathode electron potential module, proton potential module, and multiphase module. Solving the cathode component module and anode component module to obtain the mass fractions of each component. The convective term in the component equation contains the gas phase velocity, which needs to be obtained by solving the fluid flow and heat transfer module, and the source term of the component equation needs to be calculated through the reaction rates in the anode electron potential module and cathode electron potential module. Solving the anode electron potential module and cathode electron potential module to obtain the electron potential distributions of the cathode and anode, and solving the proton potential module to obtain the proton potential distributions of the anode and cathode catalyst layers and the membrane. Solving the liquid water pressure module to obtain the liquid water pressure. Solving the multiphase module to obtain the phase fraction. The convective term in the secondary phase volume fraction equation contains the mixture velocity, which needs to be obtained through the fluid flow and heat transfer module.

[0117] As Figure 1 shown, the overall software architecture can be divided into a custom layer open to users and a bottom-layer architecture invisible to users. The functions of the custom layer are implemented through custom functions. Therefore, custom functions include five types: custom source terms, custom diffusion coefficients, custom profile functions, solution strategy interfaces before and after iteration, and custom initialization. Among them, the solution strategy interfaces before and after iteration implement the information transfer between the custom layer and the bottom layer.

[0118] The bottom-layer architecture contains two data structures, namely the storage-oriented data structure and the solution-oriented data structure. The storage-oriented data structure is the bottom-layer data packet that stores the following data:

[0119] (1) Temperature field, velocity field, pressure field, physical properties (including but not limited to density, viscosity, specific heat, etc.). The temperature field, velocity field, and pressure field are obtained by solving the fluid flow and heat transfer module, and the physical properties are essential parameters for solving the fluid flow and heat transfer module;

[0120] (2) Component data packet, which is composed of custom scalar-vector equation data packets. The data of the anode component module and cathode component module are both stored in the component data packet;

[0121] (3) Custom scalar / vector equation data packet. The data of the cathode electron potential module, anode electron potential module, proton potential module, membrane water content module, liquid water pressure module, and multiphase module are all stored in the custom scalar / vector equation data packet;

[0122] (4) Custom scalar / vector field, which is the information transmitted between the custom layer and the underlying layer;

[0123] (5) Custom profile line, which sets boundary conditions for the fluid flow and heat transfer module, anode component module, cathode component module, cathode electron potential module, anode electron potential module, liquid water pressure module, and multiphase module. In the user-defined layer, the boundary conditions are set through the custom profile line function, and the custom profile line is registered with the underlying architecture, and the data is stored in the custom profile line data packet.

[0124] The solution-oriented data structure includes the flow field solver and temperature solver for the fluid flow and heat transfer module, as well as the custom scalar / vector equation solver. The general solver input for CFD is configured for the flow field solver and temperature solver by writing an xml file. All other modules except the fluid flow and heat transfer module are solved by the custom scalar / vector equation solver.

[0125] The functions of each module are as follows:

[0126] (1) Fluid flow and heat transfer module: Solve the continuity equation, momentum equation, and energy equation. The continuity equation and momentum equation are solved in all flow regions (i.e., fluid regions), and the energy equation is solved in all computational domains (including fluid regions and solid regions). Among them, the pressure-velocity coupling algorithms used to solve the continuity equation and momentum equation include but are not limited to SIMPLE, IDEAL, etc.

[0127] (2) Cathode component module: The cathode component system includes hydrogen and water vapor. The component equation of cathode hydrogen is as follows, and the equation is solved in the cathode flow region.

[0128]

[0129] Where: is the mass fraction of hydrogen; U g is the gas-phase velocity; is the effective diffusion coefficient of the component in the porous medium; is the source term of hydrogen.

[0130] (3) Anode component module: The anode component system includes oxygen and water vapor (whether to include water vapor can be selected according to the working conditions). The component equation of anode oxygen is as follows, and the equation is solved in the anode flow region.

[0131]

[0132] Where: is the mass fraction of oxygen; U g is the gas-phase velocity; is the effective diffusion coefficient of the component in the porous medium; is the source term of oxygen.

[0133] (4) Anode electron potential module: The anode electron potential equation is as follows and is solved in the anode region except for the flowing area.

[0134]

[0135] where: φ ele is the electron potential energy; is the equivalent conductivity in the porous medium material; S ele is the charge generation rate.

[0136] (5) Cathode electron potential module: The cathode electron potential equation is as follows and is solved in the cathode region except for the flowing area.

[0137]

[0138] where: φ ele is the electron potential energy; is the equivalent conductivity in the porous medium material; S ele is the charge generation rate.

[0139] (6) Proton potential module: Solve the proton potential equation in the membrane, cathode, and anode catalyst layers.

[0140]

[0141] where: φ ele is the electron potential energy; is the equivalent conductivity in the porous medium material; S ele is the charge generation rate.

[0142] (7) Membrane water content module: It is possible to choose whether to solve the membrane water content equation. If solved, it is solved in the membrane, cathode, and anode catalyst layers. Given the working characteristics of the electrolytic cell, it can also be considered that the membrane is fully hydrated, that is, the equation is not solved, and a fixed value is taken for the membrane water content in the membrane, cathode, and anode catalyst layers.

[0143] (8) Liquid water pressure module: Solve the liquid water pressure equation in the cathode catalyst layer, microporous layer, and gas diffusion layer.

[0144] (9) Multiphase module: Solve the secondary phase volume fraction equation of the mixture model and a series of algebraic expressions related to the slip velocity. The secondary phase volume fraction equation is as follows and is solved in the anode catalyst layer, microporous layer, gas diffusion layer, and flow channel.

[0145]

[0146] where: α sec is the secondary phase volume fraction; ρsec is the secondary phase density; U m is the mixed phase velocity; U dr,sec is the secondary phase slip velocity; S v-l is the gas-liquid phase change source term.

[0147] Among the above nine modules, the cathode component module, anode component module, anode electron potential module, cathode electron potential module, proton potential module, membrane water content module, liquid water pressure module and multiphase module are all solved by a user-defined scalar / vector equation solver. The main function of the user-defined scalar / vector equation solver is to define a scalar / vector equation according to requirements and solve it. Each user-defined scalar / vector equation can set the corresponding solution region. The user-defined scalar / vector equation solver essentially solves the convection-diffusion equation, and the specific solution method is similar to that of the flow and heat transfer module. All user-defined scalar / vector equations and their solution regions are stored in the form of a two-dimensional dynamic array. One dimension of the array stores the labels of the equations to be solved, and the other dimension stores the labels of the solution regions corresponding to these equations. By adding the label of the region where a certain equation needs to be solved into the two-dimensional dynamic array, the solution of this control equation in this region is realized, so as to set the corresponding solution region for each user-defined scalar / vector equation.

[0148] The general form of the user-defined scalar / vector equation is as follows. From left to right, the equation is the transient term, convection term, diffusion term and source term:

[0149]

[0150] Since the user-defined scalar / vector equation involves a series of data such as the gradients of user-defined scalars / vectors, diffusion coefficients (Γ φ in the above formula), source terms, etc., a user-defined scalar / vector equation data packet is also defined and added to the data structure used for storage at the software bottom layer. In addition, relevant functions are defined to access or modify the diffusion coefficients, source terms, etc. of the user-defined scalar / vector equation. It should be noted that physical fields such as density field, velocity field, temperature field, pressure field, viscosity, specific heat, etc. are directly included in the software bottom layer data packet and do not belong to the user-defined scalar / vector equation data packet.

[0151] Among them, for the cathode component module and anode component module, the software also defines the component data packets for the cathode and anode respectively. The component data packet consists of two user-defined scalar / vector data packets, corresponding to the component systems of the cathode and anode respectively.

[0152] During the solution Figure 1When solving the nine modules shown, a series of intermediate physical fields are involved. Therefore, a custom scalar / vector field module is also added to the custom layer of the software. The custom scalar / vector field is similar to the custom scalar / vector equation, but the difference is that the custom scalar / vector field is only used as a physical field and not the physical quantity to be solved by the equation. Similarly, relevant functions are also defined to access or modify the values, gradients, etc. of the custom scalar / vector field.

[0153] When solving Figure 1 the nine modules shown, boundary conditions need to be set for each control equation. Therefore, a custom profile is also added to the custom layer of the software. The complex boundary conditions are set through the custom profile function. The custom profile function is located in the user-defined layer. Inside the custom function, the elements on the boundary can be traversed and the face-centered values of the elements can be assigned, so as to achieve the setting of boundary conditions. Then, the defined boundary conditions are configured into the solver by registering the custom profile function.

[0154] There are a large number of control equations in the electrolytic cell, and there are complex coupling relationships between physical quantities. Specific functions are required to handle these coupling relationships and access or modify the values of the above-mentioned custom scalar / vector fields. Therefore, a series of custom functions are added to the custom layer of the software to access or modify the body-centered values, face-centered values of each physical quantity, traverse regions, elements, and faces, etc. Since the values of some custom scalar / vector fields need to be updated during the iteration process, relevant custom functions are added to the custom layer of the software as the solution strategy interface before and after the iteration. Before solving the control equation, the physical field needs to be initialized. Therefore, a custom initialization function is added to the custom layer of the software.

[0155] The software can achieve adaptive primary and secondary phases according to the volume fractions of the gas phase and the liquid phase in a connected domain. The process is as Figure 2 shown. The secondary phase volume fraction equation is as follows. From left to right, they are the transient term, the convection term, the source term related to the slip velocity, and the gas-liquid phase change source term:

[0156]

[0157] where: α sec is the secondary phase volume fraction; ρ sec is the secondary phase density; U m is the mixture phase velocity; S v-l is the gas-liquid phase change source term; U dr,sec is the secondary phase slip velocity, that is, the velocity of the secondary phase relative to the mixture phase, which can be expressed by the following formula:

[0158] U dr,sec = U sec - U m

[0159] When the phase fraction crosses 0.5, the primary phase and the secondary phase are reversed, and the control equations are automatically modified. The specific method is as follows:

[0160] First, determine the levels of the gas-phase and liquid-phase volume fractions within the unit. If the liquid-phase volume fraction is greater than the gas-phase, the gas phase should be taken as the secondary phase. When setting the secondary-phase volume fraction equation in the user-defined scalar / vector equation solver, the convective term uses the gas density, and the source term related to the slip velocity in the secondary-phase volume fraction equation also uses the gas slip velocity. Then, solve the gas-phase volume fraction equation. Conversely, solve the liquid-phase volume fraction equation. It should be noted that after knowing the slip velocity and phase fraction of one phase, the slip velocity and phase fraction of the other phase as the secondary phase can be obtained through a series of coupling relationships. Even if the adjacent unit solves the volume fraction of the other phase, information transfer can be carried out. The above coupling relationships are as follows:

[0161]

[0162] Where: U sec,p is the relative velocity between the secondary phase and the primary phase (U p ), and is defined as follows:

[0163] U sec,p = U sec - U p τ sec is the relaxation time; d sec is the particle size of the secondary-phase bubbles / droplets; μ p is the dynamic viscosity of the primary phase; f drag is the drag function; Re is the Reynolds number; is the acceleration; is the gravitational acceleration. As Figure 1 shown, the user-defined layer of the software is open to users. Users can make relevant settings in the solution strategy interface before the program runs in the software development according to the specific parameters of the electrolytic cell, define the required physical quantities and functions. The user-defined operations that can be performed include user-defined initialization, user-defined profile function, user-defined diffusion coefficient, user-defined source term, and perform the required calculations in the solution strategy interface before and after iteration.

[0164] As Figure 1 shown, the solution results of the equations in the nine modules within the software framework correspond to different physical fields such as the temperature field, pressure field, velocity field, and electric potential field within the proton exchange membrane electrolytic cell. As Figure 3 shown, the solution steps of the multi-physical fields within the proton exchange membrane electrolytic cell of the present invention are as follows:

[0165] Step 1, first initialize. Assign initial values to some scalar / vector fields through the user-defined function in the user-defined layer, and also assign initial values to the flow field temperature field solver and related variables.

[0166] Step 2: Update physical properties, source terms, boundaries, etc. in the solution strategy interface before the iteration of the custom layer.

[0167] Step 3: Solve the continuity equation and the momentum equation. There is a source term related to the slip velocity in the momentum equation. Therefore, the coupling with the multiphase module is involved in solving the momentum equation. Since the calculation of the drag function is involved in solving the coupling relationship of the relative velocity, and the drag function is related to the relative velocity, the relative velocity needs to be solved iteratively. Specifically: First, input the average velocity of the mixed phase and record the relative velocity of the previous iteration step in the slip velocity model. Use a series of coupling relationships to calculate the new relative velocity. When the iteration ends, calculate the source term related to the slip velocity and add it to the momentum source term. In the above description, the relative velocity refers to the velocity of the secondary phase relative to the primary phase, and the slip velocity refers to the velocity of the secondary phase relative to the mixed phase.

[0168] Step 4: Solve the energy equation.

[0169] Step 5: Solve the component equations, with hydrogen being solved at the cathode and oxygen being solved at the anode.

[0170] Step 6: Solve six custom scalar equations for the anode potential, cathode potential, proton potential, membrane water content, liquid water pressure, and secondary phase volume fraction. The membrane water content equation may not be solved.

[0171] Step 7: Calculate and update the physical quantities to be output in the solution strategy interface after the iteration of the custom layer, and calculate the residuals.

[0172] Step 8: If convergence is achieved, the calculation ends and the results are output. If not, return to Step 2.

[0173] The specific method steps for initializing values in Step 1 are as follows: In the defined initialization function in the custom layer, traverse the areas that need to be initialized. In each area that needs to be initialized, traverse the cells in the area to assign the specified initial values.

[0174] In Step 2, the specific method steps for updating physical properties, source terms, and boundaries are as follows: In the solution strategy interface function before the iteration of the custom layer, traverse the areas that need to update the unit information such as physical properties and source terms. In each area that needs to be updated, traverse the cells in the area. Since the information to be updated is stored in the form of custom scalar / vector fields, updating the corresponding custom scalar / vector fields of the cells can achieve the update of the unit information such as physical properties and source terms. For the update of boundaries, also traverse the boundaries that need to be updated in the solution strategy interface function before the iteration of the custom layer, traverse the faces on each boundary that needs to be updated, and assign the values to be updated to the custom scalar / vector fields of the cell faces.

[0175] In step 3, the continuity equation is as follows:

[0176]

[0177] Where: ρ m is the density of the mixed phase; S m is the mass source term.

[0178] The density and velocity of the mixed phase can be expressed as follows:

[0179]

[0180] In step 3, the momentum equation is as follows:

[0181]

[0182] Where: S u is the momentum source term.

[0183] In step 3, the coupling relationships of the slip velocity, relative velocity, and mixed-phase velocity are as follows:

[0184]

[0185] Where: U sec,p is the relative velocity of the secondary phase and the primary phase (U p ), and is defined as follows:

[0186] U sec,p = U sec - U p

[0187] τ sec is the relaxation time; d sec is the particle size of the secondary-phase bubbles / droplets; μ p is the dynamic viscosity of the primary phase; f drag is the drag function; Re is the Reynolds number; is the acceleration; is the gravitational acceleration.

[0188] In steps 4, 5, and 6, the solution methods for the energy equation and a series of scalar-vector equations are similar. The above equations can be expressed in the following general form. From left to right, the equation is the transient term, convection term, diffusion term, and source term:

[0189]

[0190] When solving the diffusion equation in the above general form, each term needs to be discretized. The discrete form of the transient term is as follows:

[0191]

[0192] It should be noted that the above discrete format is the first-order forward difference. The discrete formats that can be implemented in the present invention include but are not limited to the first-order forward difference, the first-order backward difference, and the second-order format, etc.

[0193] The discrete form of the convective term is as follows:

[0194]

[0195] The explicit discrete form of the diffusive term is as follows:

[0196]

[0197] It should be noted that the above discrete format is the explicit discrete format of the diffusive term. An implicit discrete format can also be implemented in the present invention.

[0198] The discrete form of the source term is as follows:

[0199] s φ (φ) = S p φ + S c

[0200]

[0201] It should be noted that the above discrete form of the source term is only one of them. The discrete forms of the source term that can be implemented in the present invention include but are not limited to the above form.

[0202] Rearranging the above discrete equations can obtain the following algebraic equation set:

[0203] (A T + A C - A D - A S )·[Φ] = (b T + b C - b D - b S )

[0204] Where: the subscripts T, C, D, and S represent the unsteady term, the convective term, the diffusive term, and the source term, respectively.

[0205] The rearranged algebraic equation set can be solved by numerical methods including but not limited to the Jacobi iterative method, the GS method, the over-relaxation SOR method, the GMRES algorithm, the BiCGStab algorithm, etc.

[0206] In step 7, the specific method steps for calculating and updating the physical quantities to be output are as follows. In the iterative post-solving strategy interface in the custom layer, for the centroid values of the elements that need to be calculated and updated, traverse the areas that need to be updated or calculated. In each area that needs to be updated or calculated, traverse the elements in the area. Since the information that needs to be updated or calculated is stored in the form of a custom scalar / vector field, calculating or updating the custom scalar / vector field corresponding to the element can achieve the calculation and update of the physical quantities to be output. For the face centroid values of the elements that need to be calculated and updated, traverse the internal or external boundary faces that need to be updated or calculated. In each boundary face that needs to be updated or calculated, traverse the element faces on the boundary. Since the information that needs to be updated or calculated is stored in the form of a custom scalar / vector field, calculating or updating the custom scalar / vector field corresponding to the element face can achieve the calculation and update of the physical quantities to be output. There is a specified function in the underlying architecture to calculate the residual after each iteration.

[0207] In step 8, the method for determining the end of the calculation is that the calculation ends if the number of iterations reaches the maximum value set by the user or the residual is less than the set value for terminating the calculation.

Claims

1. A software architecture for predicting the multi-physical field performance of a proton exchange membrane electrolytic cell electro-hydrothermal gas, characterized in that, It includes a custom layer open to users and an underlying architecture invisible to users; The functions of the custom layer are implemented through custom functions; The underlying architecture includes a storage-oriented data structure, a CFD general solver input, and a solution-oriented data structure; The storage-oriented data structure is an underlying data packet for storing data generated during various solution processes; The CFD general solver input is used to receive data from the storage-oriented data structure and write an xml file for configuration to the solution-oriented data structure; The solution-oriented data structure is used to solve the xml file configuration to obtain different physical fields in the proton exchange membrane electrolyzer; 2. The software architecture for proton exchange membrane electrolyzer electro-hydrothermal management analysis and design according to claim 1, characterized in that, The custom functions include five types: custom source term, custom diffusion coefficient, custom profile function, solution strategy interfaces before and after iteration, and custom initialization function. Among them, the solution strategy interfaces before and after iteration implement information transfer between the custom layer and the underlying architecture; Inside the custom profile function, the cells on the boundary can be traversed and the face-centered values of the cells can be assigned to set the boundary conditions. Then, the defined boundary conditions are configured into the solver by registering the custom profile function. The solver is located in the underlying architecture and includes a custom scalar / vector equation solver and also a CFD general solver for solving the temperature field and flow field; The custom initialization function initializes the physical field before solving the control equation; 3. The software architecture for proton exchange membrane electrolyzer electro-hydrothermal management analysis and design according to claim 1, characterized in that, The solution-oriented data structure includes a flow and heat transfer module, a cathode component module, an anode component module, a cathode electron potential module, an anode electron potential module, a proton potential module, a membrane water content module, a liquid water pressure module, and a multiphase module. The above modules all correspond to different control equations and there are couplings between different modules. The solution results of the equations in the modules correspond to different physical fields in the proton exchange membrane electrolyzer; 4. A software architecture for proton exchange membrane electrolyzer electro-hydrothermal management analysis and design according to claim 3, characterized in that, The solution of multiple physical fields in the proton exchange membrane electrolyzer includes the following steps: Step (1), assign initial values to some custom scalar / vector fields through the custom functions in the custom layer, and also assign initial values to the flow field temperature field solver and related variables; Step (2), update the physical properties, source terms, and boundaries in the pre-iteration solution strategy interface of the custom layer; Step (3), solve the continuity equation and momentum equation. There is a source term related to the slip velocity in the momentum equation. First, input the average velocity of the mixed phase and record the relative velocity of the previous iteration step in the slip velocity model. Use a series of coupling relationships to find the new relative velocity. When the iteration ends, calculate the source term related to the slip velocity and add it to the momentum source term; Step (4): Solve the energy equation; Solve the component equations, solve hydrogen at the cathode and oxygen at the anode; Solve six custom scalar equations of anode potential, cathode potential, proton potential, membrane water content, liquid water pressure, and secondary phase volume fraction. The membrane water content equation may not be solved; Step (5), calculate and update the physical quantities to be output in the post-iteration solution strategy interface of the custom layer and calculate the residuals; Step (6), if it converges, the calculation ends and the results are output. If it does not converge, return to step 2.

5. A software architecture for proton exchange membrane electrolyzer electro-hydrothermal management analysis and design according to claim 4, characterized in that, The specific method steps for initial value assignment in step (1) are as follows: in the custom initialization function in the custom layer, traverse the areas that need to be initialized. In each area that needs to be initialized, traverse the cells in the area to assign the specified initial values. In step (2), the specific method steps for updating physical properties, source terms, and boundaries are as follows: In the iteration pre-solving strategy interface function in the custom layer, traverse the areas of cell information where physical properties and source terms need to be updated. In each area that needs to be updated, traverse the cells in the area to update the cell information of physical properties and source terms by updating the custom scalar / vector field corresponding to the cell; for boundary updates, also in the iteration pre-solving strategy interface function in the custom layer, traverse the boundaries that need to be updated, traverse the faces on each boundary that needs to be updated, and assign the values that need to be updated to the custom scalar / vector field of the cell face.

6. A software architecture for proton exchange membrane electrolyzer electro-hydrothermal management analysis and design according to claim 4, characterized in that In step (3), the continuity equation is as follows: where: ρ m is the density of the mixed phase; S m is the mass source term; The density and velocity of the mixed phase are expressed as follows: The momentum equation is as follows: Where: S u is the momentum source term; The coupling relationships of slip velocity, relative velocity, and mixed-phase velocity are as follows: Where: U sec,p is the relative velocity between the secondary phase and the primary phase (U p ), and is defined as follows: U sec,p = U sec - U p τ sec is the relaxation time; d sec is the particle size of the secondary phase bubbles / droplets; μ p is the dynamic viscosity of the primary phase; f drag is the drag function; Re is the Reynolds number; is the acceleration; is the gravitational acceleration; In step (4), the solution methods for the energy equation and a series of scalar-vector equations are similar. The above equations are expressed in the following general form. From left to right, the equation is the transient term, convection term, diffusion term, and source term: When solving the diffusion equation in the above general form, each term needs to be discretized. The discretized form of the transient term is as follows: It should be noted that the above discretization format is first-order forward difference; The discretized form of the convection term is as follows: The explicit discretized form of the diffusion term is as follows: The discretized form of the source term is as follows: s φ (φ) = S p φ + S c Rearranging the above discretized equations gives the following algebraic equations: (A T +A C -A D -A S )·[Φ]=(b T +b C -b D -b S ) Where: the subscripts T, C, D, and S represent the unsteady term, convection term, diffusion term, and source term respectively.

7. A software architecture for proton exchange membrane electrolyzer electro-hydrothermal gas management analysis and design according to claim 4, characterized in that In step (5), the specific method steps for calculating and updating the physical quantities to be output are as follows: in the iteration post-solving strategy interface in the custom layer, for the cell center values that need to be calculated and updated, traverse the areas that need to be updated or calculated. In each area that needs to be updated or calculated, traverse the cells in the area, and calculate or update the custom scalar / vector field corresponding to the cell to achieve the calculation and update of the physical quantities to be output; For the cell face center values that need to be calculated and updated, traverse the internal or external boundary faces that need to be updated or calculated. In each boundary face that needs to be updated or calculated, traverse the cell faces on the boundary, and calculate or update the custom scalar / vector field corresponding to the cell face to achieve the calculation and update of the physical quantities to be output. There is a specified function in the underlying architecture to calculate the residual after each iteration; In step (6), the method for determining the end of the calculation is that if the number of iterations reaches the maximum value set by the user, or the residual is less than the set value for terminating the calculation, the calculation ends.

8. A software architecture for proton exchange membrane electrolyzer electro-hydrothermal management analysis and design according to claim 4, characterized in that, The cathode component module, anode component module, anode electronic potential module, cathode electronic potential module, proton potential module, membrane water content module, liquid water pressure module, and multiphase module are all solved by a custom scalar / vector equation solver; All custom scalar / vector equations and their solution regions are stored in the form of a two-dimensional dynamic array. One dimension of the array stores the labels of the equations to be solved, and the other dimension stores the labels of the solution regions corresponding to the equations. By adding the label of the region where a certain equation needs to be solved to the two-dimensional dynamic array, the solution of this governing equation in this region is realized, thereby achieving the setting of the corresponding solution region for each custom scalar / vector equation; The general form of the custom scalar / vector equation is as follows. From left to right, the equation consists of a transient term, a convective term, a diffusion term, and a source term: Among them, for the cathode component module and the anode component module, the software also defines the component data packets for the cathode and anode respectively. The component data packet consists of two custom scalar / vector data packets, corresponding to the component systems of the cathode and anode respectively.

9. Application of a software architecture for proton exchange membrane electrolyzer electro-hydrothermal management analysis and design, characterized in that, The software is used to realize the adaptive primary and secondary phases according to the volume fractions of the gas phase and the liquid phase in a connected domain; The secondary phase volume fraction equation is as follows. From left to right, it consists of a transient term, a convective term, a source term related to the slip velocity, and a gas-liquid phase change source term: Where: α sec is the volume fraction of the secondary phase; ρ sec is the density of the secondary phase; U m is the velocity of the mixed phase; S v-l is the source term of gas-liquid phase change; U dr,sec is the slip velocity of the secondary phase, i.e., the velocity of the secondary phase relative to the mixed phase, which can be expressed by the following formula: U dr,sec = U sec - U m The specific method is as follows: First, judge the high and low of the gas phase and liquid phase volume fractions in the unit. If the liquid phase volume fraction is greater than the gas phase, the gas phase should be used as the secondary phase. When setting the secondary phase volume fraction equation in the custom scalar / vector equation solver, the convective term uses the gas density, and the source term related to the slip velocity in the secondary phase volume fraction equation also uses the gas slip velocity. Then, solve the gas phase volume fraction equation. Otherwise, solve the liquid phase volume fraction equation; The above coupling relationship is as follows: Where: U sec,p is the relative velocity between the secondary phase and the primary phase (U p ), and is defined as follows: U sec,p = U sec - U p τ sec is the relaxation time; d sec is the particle size of the secondary-phase bubbles / droplets; μ p is the dynamic viscosity of the primary phase; f drag is the drag function; Re is the Reynolds number; is the acceleration; is the gravitational acceleration.

10. A construction method of a software for proton exchange membrane electrolyzer electro-hydrothermal management analysis and design, characterized in that, It includes the following steps; Step 1. Build the fluid flow and heat transfer module: Set the solvers for the continuity equation, momentum equation, and energy equation, and add the mass source term, momentum source term, and energy source term; Step 2. Build the cathode and anode component modules: Set the cathode and anode component solvers. The cathode components include hydrogen and water vapor. Add the source terms and diffusion coefficients of the cathode and anode component equations. The anode components include oxygen, and it can be selected whether to include water vapor; Step 3. Build the multiphase module: Set the custom scalar solver for the secondary phase volume fraction equation, and add the source term of the secondary phase volume fraction equation; Step 4. Build the cathode and anode electronic potential modules: Set the custom scalar solvers for the cathode and anode potential equations, and add the source terms and diffusion coefficients of the cathode and anode potential equations; Step 5. Build the proton potential module: Set the custom scalar solver for the proton potential equation, and add the source terms and diffusion coefficients of the proton potential equation; Step 6. Build the liquid water pressure module: Set the custom scalar solver for the liquid water pressure equation, and add the source terms and diffusion coefficients of the liquid water pressure equation; Step 7. Build the membrane water content module: Select whether to solve the film water content equation according to the working conditions. If not, specify the film water content. If so, set the custom scalar solver for the film water content equation, and add the source term and diffusion coefficient of the film water content equation; Step 8. Add coupling relationships: Add the custom functions involved in the above modules, including but not limited to a series of custom functions such as the anodic and cathodic electrochemical reaction rates, the effective permeability of the porous medium, the phase change source term of water, the open circuit voltage, and the conductivity of the membrane.

Images