Consequence simulation method based on water electrolysis hydrogen production

Abstract

The invention discloses a consequence simulation method based on water electrolysis hydrogen production, which comprises the following steps: constructing a multi-physical field coupling simulation model comprising an electric field, a thermal field, a flow field and a stress field, and setting and calibrating a geometric structure, physical property parameters, boundary conditions and a numerical strategy to obtain base field distribution of an electrode system. A bubble dynamics whole-process model is introduced, nucleation, diffusion growth and separation of bubbles are simulated, gas-liquid two-phase dynamic feedback is achieved through a volume fraction method, bubble interface force and a thermal expansion effect are coupled, and micro-crack formation of the electrode is judged by utilizing a stress criterion and a crack initiation threshold. And carrying out curing treatment on crack evolution by adopting a virtual crack unit or a phase field method. A crack propagation process is simulated in combination with a Miner criterion and a Paris law, an acoustic source item is introduced, an acoustic wave equation is solved, and the effects of full-link prediction, early warning, reliability evaluation and life optimization of electrode damage and hydrogen production performance are achieved in the operation process of the water electrolysis hydrogen production system.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of computer modeling hydrogen production simulation, and particularly relates to a consequence simulation method based on water electrolysis hydrogen production. BACKGROUND

[0002] As an important carrier of clean energy, water electrolysis hydrogen production has become a focus of current research and application due to its clean process, high product purity and other advantages. Existing hydrogen production simulation mostly uses single physical field simulation methods, such as electric field simulation for evaluating reaction rate, thermal field simulation for analyzing heat transfer efficiency, flow field simulation for studying bubble distribution, or stress field simulation for electrode mechanical property evaluation. Such methods are often independent of each other, lacking overall consideration of the multi-physical field coupling between electric field, thermal field, flow field and stress field.

[0003] In actual operation, the nucleation, growth and detachment behavior of bubbles on the electrode surface will not only change the local concentration and temperature distribution, but also cause stress concentration and structural damage on the electrode surface, thereby affecting hydrogen evolution efficiency and electrode life. If the bubble dynamics process is not combined with multi-physical field to describe, it is difficult to accurately simulate the real behavior of water electrolysis hydrogen production system in long-term operation. Therefore, we propose a consequence simulation method based on water electrolysis hydrogen production. SUMMARY

[0004] The purpose of the present application is to solve the technical problems raised in the background art by providing a consequence simulation method based on water electrolysis hydrogen production.

[0005] To achieve the above purpose, the present application provides the following technical scheme: A consequence simulation method based on water electrolysis hydrogen production, comprising the following steps: S1, constructing a multi-physical field coupling simulation model of a water electrolysis hydrogen production system, the simulation model including electric field, thermal field, flow field and stress field, inputting electrode material parameters, electrolyte physical property parameters and operation boundary conditions in the model; S2, introducing a full-process model of bubble nucleation, growth and detachment in the simulation model, coupling calculation of bubble radius evolution with electric field strength, flow field velocity distribution and thermal conduction boundary conditions; S3, further introducing stress concentration response in the simulation model, coupling interfacial tension in the bubble dynamics process with local stress distribution on the electrode surface to obtain a stress-triggered micro-crack formation threshold; S4, associating the micro-crack formation threshold with cycle number, constructing an electrode micro-crack evolution curve, and simulating characteristic acoustic emission signals generated in the crack propagation process by using acoustic wave equation; S5. Based on the acoustic emission signal and the microcrack evolution curve, the dynamic change of the hydrogen evolution rate with the accumulation of electrode damage is predicted, and the consequence simulation results of the water electrolysis hydrogen production system under different operating conditions are formed.

[0006] S1 is specifically as follows: a three-dimensional geometric model is established based on the electrode, flow channel, ion membrane and shell, and an encrypted grid is used in the vicinity of the electrode surface; temperature and concentration-dependent physical parameters are assigned to each subdomain, and the effective medium approximation method is used to treat the electrode porous layer; control equations are set in the electric field, thermal field, flow field and stress field respectively, and coupling is achieved through source terms and physical property corrections; operating boundary conditions and initial conditions are set to determine the current density, temperature, flow rate and stress state; a numerical strategy combining step-by-step weak coupling and iterative strong coupling is adopted in combination with adaptive time step control to ensure convergence and stability; the model is calibrated with experimental data and the potential, temperature, velocity, pressure and stress fields are output as base field data.

[0007] S2 is specifically as follows: based on the local supersaturation, the bubble nucleation conditions are determined and the critical radius is calculated to trigger the bubble nucleation event; after nucleation, the growth process of the bubble radius with time is simulated according to the diffusion mass transfer control law, and the bubble growth rate is coupled with the electric field intensity, flow field velocity distribution and heat conduction conditions; when the bubble reaches a certain scale, the sum of the buoyancy and fluid shear force is compared with the interface adhesion force, and when the former is greater than or equal to the latter, it is determined that the bubble has detached from the electrode surface; and the gas-liquid two-phase distribution is dynamically updated through the volume fraction method to correct the local electrical conductivity, thermal conductivity, density and viscosity, thereby realizing the feedback coupling of the bubble to the electric field, thermal field and flow field.

[0008] S3 specifically maps the bubble interfacial tension to a local load acting on the electrode surface and dynamically updates it with the bubble radius and contact angle; introduces the thermal expansion effect caused by Joule heat and reaction exotherm into the stress calculation, and corrects the thermal stress distribution of the material in combination with the linear expansion coefficient; and adopts The stress criterion is used to calculate the equivalent stress and identify stress concentration areas. The crack initiation threshold is determined based on the material yield strength, notch coefficient, and temperature and corrosion environment correction factors. Microcrack formation is determined when the equivalent stress exceeds the threshold and the duration reaches the set conditions. Crack solidification is performed at the crack initiation location using virtual crack units or phase field methods, and the local mesh is subdivided to ensure the accuracy of crack tip stress calculation. The model is calibrated based on experimental results.

[0009] S4 is specifically: according to the current density, temperature and electrolyte flow to establish a cyclic load spectrum and as a stress environment input; the cumulative damage variable is calculated by using the Miner linear damage criterion to judge the material fatigue state; in the crack propagation stage, Paris law is introduced, the crack growth rate is calculated according to the stress intensity factor amplitude, and the crack length is determined according to the material constant; the acoustic source term is introduced at the crack tip to represent the local strain energy instantaneous release, and it is input into the acoustic wave equation as a driving term to calculate the stress wave propagation process generated by crack propagation; virtual acoustic sensors are set at the back plate or shell of the electrode, and the displacement or acceleration signals are collected, and the peak value, root mean square value, counting rate and spectral centroid acoustic characteristic parameters are extracted; the simulation predicted acoustic emission signal is compared and verified with the actual observation result through experiment, so as to realize the calibration and correction of the model.

[0010] S5 is specifically: based on The equation is used to establish an electrode dynamics model, and the parameters are calibrated through the experimental polarization curve; the dynamics model is corrected in combination with the crack evolution curve and the acoustic emission characteristics, and the effective reaction area factor and impedance parameters are used to reflect the performance attenuation of the electrode; the hydrogen evolution rate is calculated in the corrected model in combination with the distribution of the electric field, the flow field and the temperature field, and the Faraday efficiency is further calculated; a multi-dimensional consequence evaluation index system is constructed, including hydrogen production attenuation rate, life prediction, risk level and triggering working condition; simulation is carried out on typical working conditions such as steady-state operation, high-load operation, low-temperature start, high-temperature continuous operation and insufficient flow, and the hydrogen production performance curve and risk result matrix are output; the Latin hypercube sampling method is used for uncertainty analysis of the key input parameters, and the confidence interval of the hydrogen production rate and the life prediction is output; The simulation results are verified through long-term running experiments to calibrate the correction factor and the impedance parameter, finally the simulation results are output in the form of curve graph, distribution graph, acoustic emission spectrum graph and risk level graph, and a data interface is provided to access a digital twin platform or a host computer system.

[0011] The beneficial effects of the present application are: The present application can simultaneously consider various physical effects such as electric field, thermal field, flow field and stress field through the multi-physical field coupling simulation method, systematically simulate the nucleation, expansion and detachment process of the gas bubble, and further improve the prediction accuracy of the electrode reaction rate and the hydrogen production efficiency. Through the simulation process of coupling multiple physical fields, the synergistic effect of various physical effects in the system is realized, and the overall understanding of the water electrolysis hydrogen production process is further improved.

[0012] The present application introduces a micro-crack evolution model combined with bubble dynamics and stress concentration effect, which can predict the formation and expansion of micro-cracks on the electrode surface and evaluate the service life of the electrode in advance. This prediction method based on micro-crack evolution improves the hysteresis of the electrode fatigue process in the traditional technology, and has the advantage of early warning.

[0013] By combining simulations of operating conditions such as current density, temperature, and flow, this invention can predict the dynamic changes in hydrogen evolution rate and electrode damage. Through multi-condition simulation and uncertainty analysis, evaluation indicators such as hydrogen production decay rate, lifespan prediction, and risk level are derived, providing a reliable basis for performance evaluation and optimization of actual water electrolysis hydrogen production systems. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] Figure 1 This is a schematic diagram of a consequence simulation method for producing hydrogen based on water electrolysis according to the present invention. DETAILED DESCRIPTION

[0015] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0016] Example 1: Figure 1 As shown, this embodiment provides a consequence simulation method based on water electrolysis to produce hydrogen, comprising the following steps: S1. Construct a multi-physics field coupling simulation model of the water electrolysis hydrogen production system, wherein the simulation model includes electric field, thermal field, flow field and stress field, and input electrode material parameters, electrolyte physical properties and operating boundary conditions into the model; S2. Introducing a full-process model of bubble nucleation-growth-detachment into the simulation model, coupling the bubble radius evolution with the electric field intensity, flow field velocity distribution, and heat conduction boundary conditions; S3. Further introducing a stress concentration response into the simulation model, coupling the interfacial tension in the bubble dynamics process with the local stress distribution on the electrode surface, and obtaining a stress-triggered microcrack formation threshold; S4. Correlating the microcrack formation threshold with the number of cycles, constructing an electrode microcrack evolution curve, and using an acoustic wave equation to simulate a characteristic acoustic emission signal generated during crack propagation; S5. Based on the acoustic emission signal and the microcrack evolution curve, the dynamic change of the hydrogen evolution rate with the accumulation of electrode damage is predicted, and the consequence simulation results of the water electrolysis hydrogen production system under different operating conditions are formed.

[0017] S1. Construct a multi-physics field coupling simulation model of the water electrolysis hydrogen production system. The simulation model includes electric field, thermal field, flow field and stress field, and inputs electrode material parameters, electrolyte physical properties and operating boundary conditions into the model. Specifically, the model includes the following sub-steps: S110, geometry modeling and meshing: a three-dimensional geometric model of the electrolytic cell is established, and the geometric structure at least includes the porous layer of the anode and cathode, the flow channel, the ion membrane, and the shell and other main parts. The electrode surface adjacent to the membrane-electrode interface region adopts local encryption mesh, and the typical unit size is 1-10 μm, and the number of overall meshes is about 10 5 -10 7 orders of magnitude to ensure numerical convergence and result accuracy. Five to ten layers of boundary layer meshes are set on the electrode surface to analyze the bubble nucleation and mass transfer region. The flow channel width ranges from 0.5 to 2 mm, and the depth ranges from 0.5 to 2 mm; the electrode porous layer pore size ranges from 1 to 20 μm.

[0018] S120, material parameter setting: each sub-domain is given with physical property parameters varying with temperature and concentration. For example: electrolyte viscosity: at 25 °C, the viscosity η ≈ 1.0 mPa·s; at 60 °C, η ≈ 0.6 mPa·s; wherein η represents the viscosity of the electrolyte; the empirical formula of viscosity with temperature can adopt Arrhenius form: wherein is the reference viscosity, E is the activation energy, R is the gas constant, and T is the temperature (K).

[0019] electrolyte conductivity: the typical range is κ = 5-20 S / m. The change of conductivity κ with temperature T and concentration C can be expressed as: wherein κ0 is the reference conductivity, a and b are fitting coefficients, T is the temperature (K), and C is the electrolyte concentration (mass fraction); thermal conductivity: the typical range is k = 0.5-1.0 W / (m·K), and the approximate relationship with temperature T is: wherein k0 is the reference thermal conductivity (W / (m·K)), and α is the thermal conductivity temperature coefficient (W / (m·K 2 ) ).

[0020] specific heat capacity: the typical range is Cp = 3000-4200 J / (kg·K), and the approximate relationship with temperature T is: wherein is the reference specific heat capacity (J / (kg·K), is the specific heat capacity temperature coefficient (J / (kg·K 2 ) ); electrode porous layer structure parameters: porosity ε = 0.4-0.6 (dimensionless); permeability K = 10 -12 -10 -10 m 2 ; wherein ε represents the pore volume fraction, and K represents the permeability of the porous medium; the relationship between the permeability and the average pore size can be represented as: The equation approximation is: in: is the porosity (dimensionless), ranging from 0.4–0.6; is the average pore size, in m; K: permeability, in m 2 .

[0021] Interfacial tension is the interfacial force between bubbles and liquids (such as electrolytes). It determines whether bubbles can detach from the electrode surface and affects the growth and number of bubbles. The typical range is γ = 0.07–0.09 N / m. The relationship with temperature T can be approximated as: ,in is the reference interfacial tension (N / m), is the temperature coefficient of interfacial tension (N / (m·K)).

[0022] S130, Control equation setting and field coupling: The electric field adopts the steady-state conductivity equation: in, is the conductivity, is the electric potential. Boundary conditions include constant current density boundary, constant potential boundary, and insulating boundary. The electric field distribution calculation results further generate the Joule heat source term: And used as thermal field input.

[0023] The thermal field adopts the transient heat conduction equation and introduces Joule heat and reaction heat: in, is the density, is the specific heat capacity, k is the thermal conductivity, is the heat of reaction, which can be expressed as: ,in is the reaction enthalpy change of the electrode reaction , r is the reaction rate , which can be typically calculated based on The dynamic model is obtained. Boundary conditions can include constant temperature boundaries, convective heat transfer boundaries, and radiation heat transfer boundaries. The calculation results can reflect the temperature field distribution and reversely correct physical properties such as conductivity and viscosity.

[0024] In terms of flow field, incompressible equation: in, is the velocity, p is the pressure, For dynamic viscosity, F is the body force. Inlet boundary can be set with flow rate or volume flux, outlet can be set with constant pressure or zero normal stress, solid wall can be set with no-slip condition. The resulting flow field distribution can be used for mass transfer and bubble dynamics calculation.

[0025] Stress field adopts small deformation linear elasticity model, including thermal expansion term: where, is stress tensor, is strain tensor, is elastic constant tensor, a is thermal expansion coefficient, is temperature rise. Boundary conditions can include fixed constraint or external load, interface conditions need to ensure displacement and stress continuity. The calculation results can be used for crack initiation determination.

[0026] Crack initiation criterion adopts equivalent stress condition: where, is equivalent stress, is critical stress threshold (typical range 200-400 MPa); Coupling between multiple physical fields is achieved through source term and property modification. For example: Joule heat of electric field is transferred to thermal field, temperature field affects electrical conductivity and viscosity, and feedback to electric field and flow field calculation. Coupling calculation can use weak coupling iteration strategy (alternately update each field, residual threshold ) or strong coupling strategy (solve all equation groups at the same time, residual threshold ). Time step is taken as 10 -5 -10 -3 , to ensure convergence and numerical stability.

[0027] S140, run boundary condition setting: electrode boundary can be set as constant current or constant voltage. Under constant current, the constraint is satisfied: where, J is current density (A / m 2 ), I is total current (A). Typical range is current density 0.1-2.0 A / cm 2 . Under constant voltage, electrode voltage U=1.8-2.2 V can be set, covering the common water electrolysis working interval.

[0028] Inlet boundary needs to define electrolyte flow rate, temperature and concentration. Inlet velocity boundary condition can be expressed as: where, is inlet velocity (m / s), is volume flux (m3 / s), A is the inlet cross-sectional area (m 2 ). Typical inlet velocity ranges from 0.01 to 0.1 m / s, inlet temperature ranges from 25 to 80 °C, and electrolyte concentration ranges from 20 to 40 wt%.

[0029] The outer wall boundary can be set to convective heat transfer or radiative heat transfer conditions: convective heat transfer boundary: where h is the convective heat transfer coefficient (W / (m 2 ·K)), T is the wall temperature, and T is the ambient temperature. Typical h ranges from 100 to 500 W / (m 2 ·K).

[0030] Radiative heat transfer boundary: where is the surface emissivity (0.7-0.9), is the Stefan-Boltzmann constant; the outer surface of the solid is subject to a free-frozen mixed constraint to ensure the stability of the overall structure calculation.

[0031] S150, initial condition setting: the initial temperature field is generally set to a uniform distribution of 25 °C; the initial flow field is a stable laminar distribution, and the velocity gradient is gradually established from the inlet to the outlet. The initial stress field can be set to zero or include installation residual stress, and the initial electric potential field can be obtained by rapid preheating under constant working conditions.

[0032] S160, numerical strategy and convergence criterion: a numerical strategy combining stepwise weak coupling and iterative strong coupling is used, and each physical field is iterated alternately within each time step until convergence. The convergence conditions include double criteria of residual threshold and relative error: the residual threshold typically ranges from <10 -6 ; the relative change rate of cross-field variables (temperature, velocity, stress, and electric potential) needs to be less than 10 -3 .

[0033] S170, time step and stability control: the time step typically ranges from 10 -5 to 10 -3 s. An adaptive time step rule is used: when the heat source or bubble volume fraction changes abruptly, the time step is automatically reduced; when running stably, the time step is gradually increased to improve calculation efficiency. Numerical stability can be controlled within the range of Co<1 to ensure the stability of the explicit solver.

[0034] S180, model calibration: the simulation results are calibrated by experimental data. For example, at 25 °C, 0.5 A / cm 2 ​Under operating conditions, the measured terminal voltage was 1.75 V, while the simulated value was 1.72 V, with a deviation of less than 2%. The measured pressure drop was 0.12 bar, while the simulated value was 0.13 bar. This comparison of experimental and simulation data ensured the accuracy of the model's predictions.

[0035] S190, Output and Data Interface Establishment: Output potential distribution, temperature field, velocity field, pressure field, stress field, and equivalent stress field data at each time step, and store acoustic response characteristics (such as crack propagation acoustic emission signals). Establish a data interface to facilitate subsequent coupling with bubble dynamics and stress concentration crack propagation models. Data can be exported to standard formats (such as CSV and VTK) for subsequent visualization and failure prediction modeling.

[0036] S2. Introducing a full-process model of bubble nucleation, growth, and detachment into the simulation model, coupling the bubble radius evolution with the electric field intensity, flow field velocity distribution, and heat conduction boundary conditions; specifically, the following sub-steps are included: S210, bubble nucleation site setting: Set bubble nucleation sites on the electrode surface. The distribution of nucleation sites can be regular or random, depending on the roughness and wettability of the electrode surface. The density range of nucleation sites is 10 4 –10 6 m -2 , to cover the surface of common electrode materials. For Ni-based electrodes, the typical site density is about 10 5 m -2 ; For Pt / C electrodes, the site density is about 10 6 m -2 .

[0037] S220, Nucleation Criteria and Critical Radius: Bubble nucleation is triggered when the local hydrogen supersaturation reaches the nucleation condition. The critical radius can be calculated by the following formula: in, is the interfacial tension, in N / m, with a typical range of 0.07 to 0.09 N·m -1 (Hydrogen-water interface). is the pressure difference between the internal pressure of the bubble and the static pressure of the electrolyte, measured in Pa. When the bubble radius is smaller than the critical radius, the bubble tends to dissolve due to interfacial tension; when the bubble radius is larger than the critical radius, the bubble can exist stably and enter the growth stage.

[0038] In the numerical implementation, the bubble pressure can be approximated by: in is the liquid static pressure, r is the bubble radius. By calculating the relationship between bubble radius and critical radius in real time, it can be determined whether the bubble nucleation event is triggered. For example, under the condition of 25℃ and 1atm, the gas-liquid interfacial tension of hydrogen in the aqueous solution is taken as 0.072 N·m -1 When the pressure difference caused by the supersaturation is 1000 Pa, the corresponding nucleation critical radius is: Therefore, it can be seen that the critical radius decreases with the increase of the pressure difference. The criterion can accurately trigger the bubble nucleation in the simulation model.

[0039] S230, bubble growth process: the change of bubble radius with time is determined by solute diffusion and convective supply: wherein, is the bubble radius, D is the diffusion coefficient of hydrogen in the electrolyte, and the typical value is about 4.5×10 -9  m 2 / s (in the aqueous solution at 25℃). C is the local concentration of hydrogen in the electrolyte, which is determined by the electrode reaction rate and convective mass transfer; is the saturation concentration of hydrogen at the temperature and pressure, which can be calculated by Henry's law; is the density corresponding to the molar volume of hydrogen, which is about 0.089 kg·m -3 When C> , the bubble obtains the driving of supersaturation, and the bubble continues to grow; when C≈ , the bubble tends to be stable; when C< , the bubble radius decreases and may be dissolved.

[0040] For example, under the condition of 25℃, assuming that the diffusion coefficient D of hydrogen in the electrolyte is 4.5×10 -9  m 2 / s, the local concentration difference is 1.0×10 -3 mol / m 3 , the bubble radius r is 50μm, and the hydrogen density is 0.089 kg / m 3 . Substituting these parameters into the bubble radius growth rate formula, the bubble radius growth rate is about 1.0×10 -6 m / s.

[0041] S240, Bubble interaction: When the bubble spacing is less than 1.5 times the radius, the merging process is triggered, and the bubble radius is updated; under the condition of high shear flow field, large bubbles may break into multiple small bubbles. The interaction process satisfies the conservation of mass and volume. In detail, the "merging process" refers to the process of two or more adjacent bubbles merging into a larger bubble after contacting each other.

[0042] S250, Bubble detachment criterion: When the bubble on the electrode surface is subjected to the action of buoyancy , fluid shear force and interfacial adhesion , when the following conditions are met, the bubble detaches from the electrode surface: Where r is the bubble radius, θ is the contact angle (typical range 30°-120°), the interfacial tension is , and the right end represents the size of the interfacial adhesion between the bubble and the electrode surface.

[0043] S260, Multiphase coupling and volume fraction update: The volume fraction (VOF, Volume of Fluid) method is used to update the bubble distribution, and the control equation is: Where represents the volume fraction of the bubble in the calculation unit, the value range is 0-1, α=0 indicates that the unit is completely liquid phase, α=1 indicates that the unit is completely gas phase, and 0<α<1 indicates that the unit is in gas-liquid mixed state; u is the velocity vector. This equation ensures that the volume fraction of the bubble phase remains constant during time advancement.

[0044] In numerical implementation, the VOF method can track the movement and deformation of the bubble interface, and dynamically update the local volume fraction field. The obtained volume fraction distribution not only determines the spatial position of the bubble in the flow field, but also feeds back to the parameters such as electrical conductivity, thermal conductivity, density and viscosity, so as to realize the multi-physical coupling influence of the bubble on the electric field, thermal field and flow field. For example, in a calculation unit, the initial bubble volume fraction α is 0.2, i.e. the unit contains 20% gas phase and 80% liquid phase. When the local velocity field u carries the bubble upward, after a certain time step, the bubble volume fraction in the unit is updated to α=0.35, indicating that the proportion of gas phase increases. At the same time, the α value of the adjacent unit will decrease or increase accordingly to ensure the overall gas volume conservation.

[0045] S270, Numerical stability and grid re-meshing: The time step range is set to 10 -5 -10 -3s. When the bubble grows rapidly or is about to break away, the time step is automatically halved to ensure the stability of the calculation. 3 s -1 In the area with a volume fraction gradient greater than 0.2, local mesh re-division or adaptive mesh refinement is used.

[0046] S280, Experimental comparison and calibration: Calibrate the simulation results through experiments. For example, at a current density of 0.5 A / cm 2 Under the working condition of , the bubble diameter is about 50μm in the experiment, and the simulation result is 48μm, with a deviation of less than 5%; at 1.0 A / cm 2 The experimental result was 80 μm, while the simulation result was 85 μm, with a deviation of less than 10%. Calibration parameters include the diffusion coefficient D, interfacial tension γ, and contact angle θ. This experimental comparison ensures the accuracy of the simulation model under different operating conditions.

[0047] S3. Further introducing stress concentration response into the simulation model, coupling the interfacial tension in the bubble dynamics process with the local stress distribution on the electrode surface, and obtaining the stress-triggered microcrack formation threshold; specifically, the following sub-steps are included: S310, interface force mapping: Map the interface tension between the bubble and the electrode surface as a local load acting on the electrode surface. The interface force calculation formula is: Where r is the bubble radius (10 -6 –10 -3 m), is the interfacial tension (0.07–0.09 N / m), and θ is the contact angle (30°–120°). The interfacial force is dynamically updated with the bubble radius and contact angle, thus forming a distributed load acting on the electrode surface.

[0048] S320, Thermal field effect: The temperature rise effect caused by Joule heat and reaction exotherm is introduced into the stress calculation. The thermal strain expression is: ,in, is the coefficient of thermal expansion (typical range (1–2)×10 -5 K -1 ), is the temperature rise (10–100 K). This thermal strain is added to the total strain tensor: is the total strain (dimensionless), which represents the total deformation of the material during operation. is the elastic strain (dimensionless) that follows Hooke's law and is proportional to the stress, ensuring that the stress distribution correction caused by temperature rise is accurately reflected.

[0049] S330, Stress calculation and criterion: Adopt Stress as equivalent stress criterion, its calculation formula is: Where, is Equivalent stress (Pa), , , The normal stress component (Pa) of the material in the three principal stress directions, respectively. The formula reflects the shear stress of the material by the difference of the three principal stresses, which is essentially equivalent to the distortion energy criterion.

[0050] When Stress exceeds the critical threshold of the material, that is, the material occurs plastic yield or micro crack initiation. For example, if the yield strength of an electrode material is 250 MPa, when the calculated is 270 MPa, it indicates that the local area has exceeded the yield limit, which may trigger crack formation.

[0051] S340, Micro crack threshold calculation: crack initiation threshold Determined by the following formula: Where, The yield strength of the material at operating temperature (typical range 200-400 MPa), The surface notch coefficient (1.5-3.0), Temperature and corrosion environment correction factor (typical value 0.8-0.9).

[0052] S350, initiation condition: when the local Greater than And the duration exceeds the minimum residence threshold (10 -6 -10 -3 s), it is determined that the area has micro crack initiation.

[0053] S360, crack solidification treatment: in the position where crack initiation is determined, the crack is numerically solidified. The crack solidification method can be realized by virtual crack element (Cohesive Zone Model, CZM) or phase field method, wherein: CZM method: the interface force and displacement relationship is introduced on the crack surface of the electrode material, and the damage evolution law is defined. The cohesive force and displacement relationship of the crack surface can adopt linear degradation model: Where, The normal stress (Pa), is the maximum strength of the interface (typical range 100–300 MPa), is the crack opening displacement (m), is the critical separation displacement (10 -8 –10 -6 sm). The relationship between interface energy and parameters is: ,in is the interfacial energy release rate (typical range 10–100 J / m 2 ).when The crack is judged to be completely broken.

[0054] Phase field method: introducing a continuous damage field in the material domain Its value range is 0–1, =0 means complete, =1 indicates complete fracture. The phase field governing equation can be expressed as: in, is the phase field characteristic length (typical range 1–10 μm), : Laplace operator, corresponding to the diffusion of the phase field in space, is the critical energy release rate (10–100 J / m 2 ), is the positive strain energy density (J / m 3 ). Indicates the degree of remaining integrity of the material. The phase field method automatically captures the crack propagation path within the finite element mesh. Meshing and Accuracy Assurance: To ensure the accuracy of the numerical analysis of the stress gradient at the crack tip, localized meshing is performed as necessary, with a typical element size of 1–5 μm.

[0055] S370, experimental comparison and verification: Scanning electron microscopy (SEM), in-situ tensile observation or X-ray tomography are used to detect cracks. Comparison indicators include: crack initiation location (error <15%), crack density (number / area, deviation <10%), and crack length distribution (deviation <10%). For example, at 0.5 A / cm 2 Under these conditions, the crack initiation locations predicted by simulation matched actual observations by over 85%, with the average crack length deviation less than 10%. These comparisons allowed us to calibrate thresholds and model parameters.

[0056] S4, correlating the microcrack formation threshold with the number of cycles, constructing an electrode microcrack evolution curve, and using an acoustic wave equation to simulate a characteristic acoustic emission signal generated during crack propagation; specifically comprising the following sub-steps: S410, Cycle load spectrum definition: According to the operating conditions of the electrolytic cell, a cycle spectrum of current density, temperature and electrolyte flow rate varying with time is established to characterize the dynamic stress environment of the electrode material. The cycle spectrum can be expressed in the form of a time series function: wherein, is the current density (A / cm 2 ), is the temperature (K), is the flow rate (m 3 / s). In this embodiment, the current density ranges from 0.1 to 2.0 A / cm 2 , the temperature ranges from 298 to 353 K, and the flow rate ranges from 1 x 10 -6 s to 1 x 10 -4 m / s. 3

[0057] S420, Cumulative damage calculation: The material life is evaluated using the Miner linear damage criterion, which is expressed as: wherein, D is the cumulative damage variable (dimensionless); is the number of cycles experienced by the material at stress level i; is the limit number of cycles that the material can withstand at the same stress level. When the material is repeatedly loaded at different stress levels, the cumulative damage can be approximately the sum of the partial damages at each stress level. When D≥1, it is determined that the material reaches the fatigue failure limit. For example, in the cycle load spectrum of a certain electrode, the material has experienced 5,000 cycles at a high stress level, and the life limit at this stress level is 20,000 cycles; it has experienced 30,000 cycles at a low stress level, and the life limit is 100,000 cycles. Then the cumulative damage is D=5000 / 20000+30000 / 100000=0.55, and at this time D is 0.55<1, indicating that the material has not failed. When the number of cycles further increases to cause D to reach or exceed 1, it is considered that the material enters the fatigue fracture state.

[0058] S430, Crack propagation law: The Paris law is used to calculate the change of crack length with the number of cycles: wherein, is the crack length, N is the number of cycles, is the stress intensity factor amplitude, and C and m are material constants obtained by experimental calibration.

[0059] The stress intensity factor can be calculated by the following formula: ​ wherein, is the stress amplitude (Pa), a is the crack length (m), Y is the geometry correction factor (usually 1.0-1.2), for example, if the Paris constant of a certain electrode material is C=1×10 -10 , m=3, under the condition of stress amplitude 50 MPa and crack length 1 mm, the stress intensity factor amplitude is calculated to be about 2.8 MPa·m 1 / 2 . Substituting into the Paris law, the crack propagation rate is 2.2×10 −9 m / cycle, that is, the crack propagation is about 2.2 mm per 10 6 cycles.

[0060] S440, acoustic source term construction: the strain energy release in the crack propagation process is converted into an acoustic source term, and the expression is: wherein, is the acoustic source term at time t (unit: J·s -1 ·m -3 ), E is the elastic modulus of the material (Pa), is the local strain at the crack propagation (dimensionless), is the volume element participating in crack propagation (m 3 ), is the Dirac function, which is used to represent the energy release in the instantaneous occurrence of crack propagation. The formula shows that the local strain energy at the moment of crack propagation is converted into acoustic energy in the form of a pulse and is input as a source term into the acoustic wave equation, which is used to simulate the propagation process of acoustic emission signals in the material. For example, for a nickel-based electrode with an elastic modulus of 200 GPa, when the local strain at the crack tip is 0.002 and the volume element participating in crack propagation is 1×10 -15  m 3 , the acoustic source pulse intensity corresponding to a single energy release is calculated to be about 4.0×10 -4 J. This energy is then propagated outward in the form of a stress wave and can be received by an acoustic emission sensor arranged on the electrode back plate.

[0061] S450, acoustic wave solving: solving the acoustic wave equation in the solid domain: wherein is the density (typical value of nickel-based electrode is about 8900 kg / m 3 ), is the displacement field, is the stress tensor. In numerical implementation: time step 10 -6-10 -3 s, with explicit or central difference method. Boundary conditions: free boundary for electrode outer wall; damping absorption boundary for shell boundary to avoid reflection. Grid scale is 10–50 μm to ensure analytical accuracy of sound wave propagation.

[0062] S460, acoustic signal acquisition and feature extraction: virtual acoustic sensors are arranged on the back plate of the electrode or the surface of the shell, and the sampling rate is set to 2 MHz. After signal acquisition, the following features are extracted: peak voltage: 10–200 μV; root mean square value: 5–50 μV; count rate: 100–1000 times / s; spectral centroid: 50–300 kHz. Signal processing includes high-pass filtering (cutoff frequency 20 kHz), FFT transform, envelope analysis, for identifying typical acoustic emission events.

[0063] S470, experimental comparison and model calibration: compare the acoustic emission experiment with the simulation results, for example, under the condition of 0.5 A / cm 2 , the measured acoustic emission peak is 120 μV, and the simulation predicted value is 115 μV, with a deviation of less than 5%. Through the comparison results, the model is calibrated to ensure the reliability of the simulation prediction.

[0064] S5, based on the acoustic emission signal and micro-crack evolution curve, the dynamic change of hydrogen evolution rate with electrode damage accumulation is predicted, and the consequence simulation results of the water electrolysis hydrogen production system under different operating conditions are formed, including the following sub-steps: S510, kinetic model establishment: establish a kinetic model for the hydrogen evolution reaction of the electrode, which can use equation: where η is the overpotential, i is the current density, a and b are fitting parameters; or use equation: where i0 is the exchange current density (typical value 10 -3 ~ 10 -1 A / cm 2 , α is the electron transfer coefficient (0.3~0.7), n is the number of electrons, F is the Faraday constant, R is the gas constant, and T is the temperature. The kinetic parameters are calibrated by the experimental polarization curve.

[0065] S520, model correction mechanism: introduce crack evolution curve and acoustic emission feature correction kinetic model: in the early stage of crack, due to the increase of active area, the effective reaction area factor rises to 1.0–1.2; with crack propagation, the factor drops to 0.5–0.8; impedance rises is 10-30%. The function of ohmic impedance with crack length L can be expressed as: wherein, is the corrected impedance (Ω), is the initial impedance (Ω), L is the real-time crack length, L0 is the initial crack length, k is the experimental fitting constant (0.1-0.3), which is obtained by fitting the experimental results of impedance under different crack lengths by the least square method.

[0066] S530, hydrogen production rate and efficiency calculation: based on the corrected kinetic model, and combined with the distribution of electric field, flow field and temperature field, the evolution curve of hydrogen evolution rate with time is calculated. Further calculation of Faraday efficiency: wherein is the hydrogen production volume, I is the input current, t is the running time, and the formula represents the level of charge utilization in the electrolysis process, i.e. the proportion of unit charge converted into hydrogen molecules. When approaches 1, it means that all the current is used for hydrogen evolution; when is less than 1, it means that there are side reactions or electrode losses.

[0067] S540, consequence evaluation index: a multi-dimensional index system is established, including: hydrogen production decay rate: the percentage decrease of hydrogen production rate per unit time; life prediction: when the hydrogen production rate decreases by more than 10%, it is judged as mild failure, and when it decreases by more than 20%, it is judged as serious failure; risk level: based on acoustic emission signals (such as RMS, count rate) and crack evolution speed for grading determination; trigger working condition: record the current density, temperature and flow conditions that cause performance degradation or risk level improvement.

[0068] S550, typical working condition simulation: simulation is established for various working conditions: steady-state operation: current density 0.3-0.6 A / cm 2 , temperature 25-40 ℃; high load operation: current density 1.0-2.0 A / cm 2 , temperature 60-80 ℃; low temperature start-up: temperature <10 ℃, initial heating rate 0.5-2 K / min; high temperature continuous operation: temperature ≥85 ℃, time ≥100 h; insufficient flow: electrolyte flow <0.01 L / min. Output hydrogen production performance curve, life prediction curve and risk level matrix.

[0069] S560, uncertainty analysis: Latin hypercube sampling is performed for key input parameters (such as crack formation threshold, nucleation site density, contact angle, membrane conductivity) with no less than 1000 times, assuming that the parameters follow normal or uniform distribution, and output the 95% confidence interval of hydrogen production rate and life prediction.

[0070] S570, experimental comparison and model calibration: the simulation results are verified by long-term running experiments. For example, under the condition of 0.5 A / cm 2 , 60 ℃, the measured hydrogen production rate after 100 hours is attenuated by 12%, and the simulation prediction value is 11%, with an error of less than 10%. This experimental comparison is used to calibrate the correction factor and impedance parameter.

[0071] S580, result visualization and interface output: the simulation results are output in the form of curve graph, distribution graph, acoustic emission spectrum graph and risk level radar graph, and data interface is provided to access digital twin platform or host computer, realizing real-time monitoring and online early warning, and the output data can be in CSV / JSON format, and supports OPC UA protocol to interface with industrial control system.

[0072] Embodiment two: the embodiment provides a consequence simulation system based on electrolytic water hydrogen production, comprising: a modeling module: for establishing a three-dimensional geometric model containing electrodes, flow channels, ion membranes and housings; the module supports meshing and physical property assignment for electric field, thermal field, flow field and stress field, and uses encrypted mesh on the local electrode surface; a multi-physical field coupling module: for coupling electric field, thermal field, flow field and stress field; the module introduces source terms and correction terms in the electric field equation, energy equation, equation and mechanical equilibrium equation, realizing the interaction of multi-physical fields; a bubble dynamics simulation module: for simulating the whole process of bubble nucleation, growth and detachment on the electrode surface; including bubble critical radius judgment unit, growth rate calculation unit and bubble detachment criterion unit; using volume fraction (VOF) method to dynamically update the distribution of gas-liquid two-phase, and feeding back to electric field, flow field and thermal field calculation; a crack evolution simulation module: for simulating the stress concentration caused by bubble interfacial force and thermal expansion effect, and predicting the initiation and propagation of electrode micro-cracks; including stress criterion unit, micro-crack threshold calculation unit, crack solidification unit; in the crack propagation stage, Miner rule and Paris law are used to calculate fatigue accumulation and crack growth rate; an acoustic simulation module: for simulating the strain energy release and acoustic emission signal propagation caused by crack propagation; including acoustic source term construction unit and acoustic wave equation solving unit; The peak voltage, root mean square value, count rate, and spectral centroid can be output at the virtual sensor position. Hydrogen production prediction module: based on Equations and The equations establish an electrode kinetics model. The kinetics parameters are corrected by combining crack evolution and acoustic emission characteristics to predict hydrogen evolution rate and electrode performance degradation. The current conversion rate is calculated through the Faraday efficiency formula.

[0073] Consequence simulation and evaluation module: used to output hydrogen production performance and electrode life prediction results under different working conditions. Including multi-working condition simulation unit (steady state, high load, low temperature start-up, etc.), risk level evaluation unit and uncertainty analysis unit (using Latin hypercube sampling). The output results can be visualized as hydrogen production curves, acoustic emission spectra, crack propagation distribution maps and risk matrices.

[0074] Data interface module: used to connect simulation data to digital twin platform or host computer system, and realize docking with actual running system.

[0075] The above formulas are all dimensionless numerical calculations, the formulas are obtained by collecting a large amount of data to simulate the most recent real situation, and the preset parameters and threshold values in the formulas are set by the person skilled in the art according to the actual situation.

[0076] The above embodiments can be realized wholly or partially by software, hardware, firmware or any other combination. When realized by software, the above embodiments can be realized wholly or partially in the form of a computer program product. The computer program product includes one or more computer instructions or computer programs. When the computer instructions or computer programs are loaded or executed on a computer, the processes or functions described in the embodiments of the present application are wholly or partially generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network or other programmable devices. The computer instructions can be stored in a computer readable storage medium or transferred from one computer readable storage medium to another, for example, the computer instructions can be transferred from one website, computer, server or data center to another website, computer, server or data center through wired (such as infrared, wireless, microwave, etc.) mode. The computer readable storage medium can be any available medium that can be accessed by a computer or a data storage device such as a server, data center, etc. containing one or more available medium collections. The available medium can be a magnetic medium (such as a floppy disk, a hard disk, a magnetic tape), an optical medium (such as a DVD), or a semiconductor medium. The semiconductor medium can be a solid state disk.

[0077] Those of skill in the art would understand that the modules and algorithms described in connection with the examples described herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. The

[0078] Those of skill in the art would understand that, for the purposes of description and enabling clarity, the specific process of the system, device and module described above can refer to the corresponding process in the foregoing method embodiments, which will not be described again.

[0079] In several embodiments provided in the present application, it should be understood that the disclosed system, device and method can be implemented in other ways. For example, the device embodiments described above are merely illustrative, for example, the division of the modules is only a logical function division, and actual implementation can have another division manner, for example, multiple modules or components can be combined or integrated into another system, or some features can be ignored or not executed. In addition, the coupling or direct coupling or communication connection between the displayed or discussed modules can be indirect coupling or communication connection through some interfaces, devices or modules, which can be electrical, mechanical or other forms.

[0080] The modules described as separate components can or can not be physically separated, and the components displayed as modules can or can not be physical modules, which can be located in one place or distributed on multiple network modules. Part or all of the modules can be selected according to actual needs to achieve the purpose of the embodiment scheme.

[0081] In addition, the functional modules in each embodiment of the present application can be integrated in one processing module, or each module can exist physically, or two or more modules can be integrated in one module.

[0082] If the functions are implemented in the form of software function modules and sold or used as independent products, they can be stored in a computer readable storage medium. Based on this understanding, the technical solutions of the present application essentially or the parts that contribute to the prior art or parts of the technical solutions can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes a number of instructions for causing a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of the present application. The aforementioned storage medium includes: a U disk, a mobile hard disk, a read-only memory (ROM), a random access memory (RAM), a magnetic disk or an optical disk, and various media that can store program codes.

[0083] The above is only a specific implementation of the present application, but the protection scope of the present application is not limited thereto. Any person skilled in the art can easily think of changes or replacements within the technical scope disclosed in the present application, which should be included in the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claims.

[0084] Finally: the above is only a preferred embodiment of the present application and is not intended to limit the present application. Any modification, equivalent replacement, improvement, etc. made within the spirit and principles of the present application should be included in the protection scope of the present application.

Claims

1. A consequence simulation method based on water electrolysis to produce hydrogen, characterized in that: The steps include: S1. Construct a multi-physics field coupling simulation model of the water electrolysis hydrogen production system, wherein the simulation model includes electric field, thermal field, flow field and stress field, and inputs electrode material parameters, electrolyte physical properties parameters and operating boundary conditions into the model; S2. Introducing a full-process model of bubble nucleation, growth, and detachment into the simulation model, coupling the bubble radius evolution with the electric field intensity, flow field velocity distribution, and heat conduction boundary conditions; S3. Further introducing a stress concentration response into the simulation model, coupling the interfacial tension in the bubble dynamics process with the local stress distribution on the electrode surface, and obtaining a stress-triggered microcrack formation threshold; S4. Correlating the microcrack formation threshold with the number of cycles, constructing an electrode microcrack evolution curve, and using an acoustic wave equation to simulate a characteristic acoustic emission signal generated during crack propagation; S5. Based on the acoustic emission signal and the microcrack evolution curve, the dynamic change of the hydrogen evolution rate with the accumulation of electrode damage is predicted, and the consequence simulation results of the water electrolysis hydrogen production system under different operating conditions are formed.

2. The method for simulating the consequences of hydrogen production based on water electrolysis according to claim 1, characterized in that: S1 is specifically: A three-dimensional geometric model is established based on the electrodes, flow channels, ion membranes and shells, and a refined mesh is used in the vicinity of the electrode surface; Temperature- and concentration-dependent physical parameters are assigned to each subdomain, and the effective medium approximation method is used to treat the electrode porous layer. Set the control equations in the electric field, thermal field, flow field and stress field respectively and realize coupling through source terms and physical property correction; set the operating boundary conditions and initial conditions to determine the current density, temperature, flow rate and stress state; A numerical strategy combining step-by-step weak coupling and iterative strong coupling is adopted, combined with adaptive time step control to ensure convergence and stability; the model is calibrated with experimental data and then the potential, temperature, velocity, pressure and stress fields are output as base field data.

3. The consequence simulation method based on water electrolysis hydrogen production according to claim 1, characterized in that: S2 is specifically: Determine the bubble nucleation conditions based on the local supersaturation and calculate the critical radius to trigger the bubble nucleation event; After nucleation, the growth of bubble radius over time is simulated according to the diffusion mass transfer control law, and the bubble growth rate is coupled with the electric field intensity, flow field velocity distribution and heat conduction conditions. When the bubble reaches a certain size, the sum of the buoyancy and fluid shear force is compared with the interfacial adhesion force. When the former is greater than or equal to the latter, it is determined that the bubble has detached from the electrode surface. The gas-liquid two-phase distribution is dynamically updated through the volume fraction method to correct the local electrical conductivity, thermal conductivity, density and viscosity, thereby realizing the feedback coupling of bubbles to the electric field, thermal field and flow field.

4. The method for simulating consequences of hydrogen production based on water electrolysis according to claim 1, characterized in that: S3 specifically: The bubble interfacial tension is mapped as a local load acting on the electrode surface and dynamically updated with the bubble radius and contact angle; The thermal expansion effect caused by Joule heat and reaction heat is introduced into the stress calculation, and the thermal stress distribution of the material is corrected in combination with the linear expansion coefficient; use The stress criterion calculates the equivalent stress and identifies the stress concentration area; the crack initiation threshold is determined based on the material yield strength, notch coefficient, and temperature and corrosion environment correction factors, and microcrack formation is determined when the equivalent stress exceeds the threshold and the duration reaches the set conditions.

5. The method for simulating consequences of hydrogen production based on water electrolysis according to claim 4, characterized in that: S3 also includes: crack solidification at the crack initiation location using virtual crack units or phase field methods, re-dividing the local grid to ensure the accuracy of crack tip stress calculation, and calibrating the model based on experimental results.

6. The consequence simulation method based on water electrolysis hydrogen production according to claim 1, characterized in that: S4 is specifically: A cyclic load spectrum is established based on current density, temperature, and electrolyte flow rate and used as stress environment input; the cumulative damage variable is calculated using the Miner linear damage criterion to determine the material fatigue state; Paris' law is introduced in the crack propagation stage to calculate the crack growth rate based on the stress intensity factor amplitude, and the evolution of crack length with the number of cycles is determined in combination with material constants; An acoustic source term is introduced at the crack tip to characterize the instantaneous release of local strain energy, and is input into the acoustic wave equation as a driving term to calculate the propagation process of the stress wave generated by crack extension.

7. The method for simulating consequences of hydrogen production based on water electrolysis according to claim 6, characterized in that: S4 also includes: setting a virtual acoustic sensor at the electrode back plate or shell position to collect the displacement or acceleration signal and extract the peak value, root mean square value, count rate and spectrum centroid acoustic characteristic parameters; The acoustic emission signals predicted by simulation are verified by experimental comparison with the actual observation results to achieve calibration and correction of the model.

8. The consequence simulation method based on water electrolysis hydrogen production according to claim 1, characterized in that: S5 is specifically: based on The electrode kinetic model was established using the equations, and the parameters were calibrated using experimental polarization curves. The dynamic model is modified by combining the crack evolution curve and acoustic emission characteristics, and the effective reaction area factor and impedance parameters are used to reflect the electrode performance attenuation; In the modified model, the hydrogen evolution rate is calculated by combining the electric field, flow field and temperature field distribution, and the Faraday efficiency is further calculated; Construct a multi-dimensional consequence assessment indicator system, including hydrogen production decay rate, life prediction, risk level and triggering conditions; Simulate typical operating conditions including steady-state operation, high-load operation, low-temperature startup, high-temperature continuous operation, and insufficient flow, and output hydrogen production performance curves and risk result matrices. The Latin hypercube sampling method is used to perform uncertainty analysis on key input parameters and output confidence intervals for hydrogen production rate and lifetime predictions.

Images