Intelligent collaborative control system and method for deep-sea hydrogen production and storage based on TPLBFS

Abstract

The application discloses a kind of deep-sea hydrogen extraction and storage intelligent collaborative control system and method based on TPLBFS, belong to deep-sea energy development and intelligent control technical field.System includes the perception layer, intelligent collaborative computing engine layer and feedback control layer connected in turn, and the perception layer passes through multi-source sensor array real-time acquisition deep-sea extraction and storage working condition data;Computing engine layer uses improved three-phase lattice boltzmann flux solver TPLBFS, through five special optimization for hydrogen extreme property, realize the stable simulation of high-density three-phase flow, mass conservation error is reduced to below 0.01%;Feedback control layer is based on simulation prediction result and generates closed-loop control instruction.The application first utilizes long-chain hydrocarbon spreading effect to improve hydrogen transmission efficiency by 40%, break the existing technical prejudice, realize 2300 times calculation acceleration by GPU parallel acceleration, system response time is less than 0.5 second, can be widely applied to deep-sea natural hydrogen efficient exploitation, submarine depleted oil and gas reservoir hydrogen storage safety evaluation and operation optimization.

Description

Technical Field

[0001] This invention relates to the field of intelligent control technology for deep-sea hydrogen extraction and geological hydrogen storage, and in particular to an intelligent collaborative control system and method for deep-sea hydrogen extraction and storage based on an improved three-phase lattice Boltzmann flux solver. Background Technology

[0002] With the advancement of the global energy transition, natural hydrogen discovered in deep-sea hydrothermal vents has become an important strategic clean energy source, and depleted seabed oil and gas reservoirs have become ideal geological carriers for large-scale hydrogen storage and transportation. However, the migration environment of deep-sea hydrogen is extremely complex, usually involving immiscible three-phase flow processes of hydrogen, seawater, and long-chain saturated hydrocarbons, and the high-pressure environment of the deep sea brings extreme differences in physical properties.

[0003] In a deep-sea environment of 3000 meters (30 MPa, 4℃), the density of hydrogen is only about 80 kg / m³. 3 Hydrogen's density ratio to seawater is as high as 1:13.75, and under normal pressure, the ratio reaches 1:12200; meanwhile, hydrogen's dynamic viscosity is only 8.9 μPa·s, two orders of magnitude lower than that of water. This extreme difference in physical properties presents insurmountable technical bottlenecks for traditional numerical simulation and engineering control. Numerical stability is poor under extreme physical property conditions. When dealing with multiphase flow systems with high density ratios, traditional volume function (VOF) and level-set methods are prone to generating severe numerical oscillations and non-physical spurious velocities at the interface. The magnitude of spurious velocities can reach 10%-50% of the real physical velocities. The condition number of the pressure Poisson equation increases exponentially with the density ratio, causing the iterative solver to fail to converge and the simulation process to collapse in a short time.

[0004] Mass conservation fails in long-term simulations. In the simulation and prediction of hydrogen storage facilities operating for decades, traditional methods often result in a cumulative mass loss of 3%-5%. For hydrogen storage projects with a capacity of millions of tons, this can lead to economic assessment errors of tens of millions of yuan, as well as incorrect estimates of hydrogen escape, causing serious safety hazards. Although existing phase-field methods can capture interface evolution, their fourth-order Cahn-Hilliard equation is prone to negative volume fractions and mass drift during long-term integration, failing to meet the requirements of long-term simulations.

[0005] The lack of dynamic correction for model-physical deviations. Traditional simulation software adopts an open-loop working mode of "offline simulation + manual analysis", which cannot automatically correct the initial geological model based on real-time monitoring data from the seabed. When the actual working conditions deviate from the design conditions, the deviation between the simulation results and the physical entity continues to accumulate, eventually leading to the complete failure of the prediction and making it unable to cope with the dynamic changes in the deep-sea environment.

[0006] The lack of real-time closed-loop control capability is a significant issue. In the 3,000-meter deep-sea environment, manual intervention is extremely costly and slow to respond. When there is a risk of hydrogen fingering or escape, control commands need to be generated and executed within hundreds of milliseconds. However, current technologies require several hours to several days from the discovery of an anomaly to the formulation of an adjustment plan, which is completely inadequate to cope with dynamic risks and cannot achieve complete closed-loop control of "perception-decision-execution".

[0007] The Lattice Boltzmann Flux Solver (LBFS) is an effective scheme for numerical simulation of high-density-ratio multiphase flows. Existing techniques have extended it to isothermal flow simulation of three-component fluids. However, for extreme conditions in deep-sea hydrogen-seawater-long-chain hydrocarbon three-phase systems, existing LBFS models suffer from poor numerical stability, inability to adapt to extreme hydrogen properties, and lack of real-time data assimilation and closed-loop control interfaces, making them unsuitable for direct application in deep-sea hydrogen production and storage engineering. Related technologies can be found in CN114562345A, which discloses a simulation method for deep-sea natural gas hydrate extraction, but does not address three-phase flow coupling. The paper "Research Progress in Numerical Simulation of Deep-Sea Hydrogen Extraction" (Ocean Engineering, 2023) reviews existing numerical methods, pointing out that high-density-ratio three-phase flow simulation remains a gap in the industry. Summary of the Invention

[0008] The purpose of this invention is to overcome the aforementioned deficiencies of the prior art and provide a smart collaborative control system and method for deep-sea hydrogen production and storage based on TPLBFS. TPLBFS stands for Three-phase Lattice Boltzmann Flux Solver. The core inventive concept and innovation of this invention (FluxSolver) lies in the following: Firstly, addressing the extreme operating conditions and existing technological bottlenecks in deep-sea hydrogen extraction and storage, an improved three-phase lattice Boltzmann flux solver (TPLBFS) is proposed and constructed. This solver unifies and couples the Navier-Stokes (NS) equation, the Cahn-Hilliard phase-field equation, and the energy equation. Through five specific optimizations targeting extreme hydrogen properties, stable and high-precision simulation of three-phase flow under extremely high density ratios is achieved. Based on this, a three-phase flow physical model adapted to deep-sea porous / fractured media and a method for quantifying hydrogen transport dynamics are established, filling the industry gap in pore-scale simulation of deep-sea hydrogen transport. Secondly, the mechanism of long-chain hydrocarbon interface spreading lubrication enhancement is revealed and engineered for the first time, overturning traditional understanding and forming a standardized hydrogen transport enhancement regulation strategy. Finally, a five-dimensional closed-loop collaborative control architecture of "sensing-assimilation-prediction-decision-execution" is constructed to achieve real-time intelligent regulation and proactive risk prevention throughout the entire process of deep-sea hydrogen extraction and geological hydrogen storage. To achieve the above objectives, this invention provides a method for intelligent collaborative control of deep-sea hydrogen production and storage based on TPLBFS, comprising the following steps: S1. Establish a physical model for deep-sea hydrogen production and storage, determine the computational domain, initialize the physical parameters and boundary conditions of the grid cells within the computational domain, and complete the initialization of the physical parameters of the grid cell interface through linear interpolation. S2. Calculate the density distribution function at the current time step based on the improved three-phase lattice Boltzmann flux solver TPLBFS. Temperature distribution function Phase distribution function The TPLBFS is coupled with the Navier-Stokes equation, the Cahn-Hilliard phase field equation, and the energy equation, and five specific optimizations are performed for the extreme properties of hydrogen. S3. Calculate the chemical potential function of the three-phase fluid, and obtain the potential function values ​​of the hydrogen phase, seawater phase, and long-chain hydrocarbon phase. m 1 ,m 2 、 m 3. Simultaneously calculate surface tension and volume force; S4. By deploying a multi-source sensor array in the deep-sea operation area, the measured data of the current working conditions are collected in real time, and the measured data is assimilated into the TPLBFS model to complete the dynamic update of the model boundary conditions and physical property parameters. S5. Using the optimized TPLBFS calculation interface, the macroscopic flux and mass correction source terms are used to solve the evolution equations through the fourth-order Runge-Kutta method to obtain the macroscopic physical parameters of the current time step, and to predict the phase distribution, velocity field and pressure field under future working conditions. S6. Calculate the dynamic relative penetration rate and spreading coefficient based on the prediction results, determine whether the current working condition deviates from the optimal operating window, and if it deviates, generate the corresponding control command and send it to the actuator to complete the closed-loop control. S7. Repeat steps S2 to S6 until the mining and storage operation cycle ends. Furthermore, in step S1, the physical model is established based on the actual geological characteristics, well logging data, and seismic exploration data of the deep-sea hydrothermal production area or the depleted oil reservoir hydrogen storage area, restoring the porous media structure, formation distribution, and well location of the target area; the computational domain is divided into three-dimensional structured grids, with the upper boundary being the caprock closed boundary, the left, right, front, and rear boundaries being periodic boundaries, and the lower boundary being the production well / injection well boundary. Navier slip boundary conditions are set for both the well wall and the formation wall; the initialized physical parameters include formation porosity, permeability, initial pressure, initial temperature, volume fraction, density, viscosity, and surface tension of each phase fluid.

[0009] The physical parameters of the mesh cell interface are initialized by linear interpolation. Specifically, based on the initialized physical parameters of the center point of the mesh cell, the physical parameter values ​​at the interface of adjacent mesh cells are calculated by bilinear interpolation (two-dimensional model) / trilinear interpolation (three-dimensional model). This ensures that physical quantities such as pressure, velocity, temperature, and phase volume fraction transition smoothly between cells, eliminates abrupt changes in physical quantities at the interface, and improves the numerical stability and calculation accuracy of the simulation.

[0010] Further, in step S2, the density distribution function f i ∗ Temperature distribution function g i ∗ Phase distribution function h i ∗,j The velocity distribution is calculated based on the equilibrium distribution function at the current moment and the distribution function at the previous moment. The lattice model adopts the D3Q19 three-dimensional discrete velocity model, and the lattice sound velocity c is... s =c / 3, where c is the unit velocity of the grid. The specific calculation formula is as follows: The essential difference between the improved three-phase lattice Boltzmann flux solver TPLBFS and the traditional LBFS is as follows: (1) The lattice Boltzmann distribution function is transformed into a macroscopic flux expression of the finite volume method through Chapman-Enskog multi-scale expansion, realizing a direct mapping from the mesoscopic distribution function to the macroscopic flux, avoiding the numerical error introduced by the additional interpolation required by the traditional LBFS when solving the macroscopic flux; (2) For the three-phase fluid system, the Navier-Stokes equation, Cahn-Hilliard phase field equation and energy equation are solved simultaneously, and the mass correction source term is introduced through the Lagrange multiplier method to ensure that the volume fraction of the three phases strictly meets the conservation constraints at each time step; (3) The complete TPLBFS algorithm flow is: distribution function initialization → collision step (BGK approximation) → migration step → flux calculation (reconstructing macroscopic flux from the distribution function) → macroscopic quantity update (density, velocity, temperature, volume fraction of each phase) → mass conservation correction → enter the next time step.

[0011] The TPLBFS method reduces the mass conservation error in long-period simulations from 3%-5% in traditional methods to below 0.01% by introducing a mass correction source term using the Lagrange multiplier λ into the phase-field equations, achieving an improvement of 300-500 times. This correction method determines the value of λ by solving the mass conservation equations at each time step. Specifically, after the collision-migration step at each time step, the deviation of the sum of the current phase volume fractions from 1 is calculated. This deviation is used as a constraint, and the optimal value of λ is obtained by iteratively solving the linear equations, thereby correcting the distribution functions of each phase and ensuring that the sum of the three phase volume fractions is always 1.

[0012] The density equilibrium distribution function at the current moment is: , in, Indicates the current position. Indicates the current moment. Indicates the particle's direction at the current moment. Indicates the weighting factor. Indicates pressure, This represents the mass density of a fluid at a specific location and time. Represents the speed of sound in a lattice. This represents the discrete velocity of the crystal lattice. This indicates the velocity of a fluid at a specific location and time. Temperature equilibrium distribution function at the current moment for: , Where T represents the macroscopic thermodynamic temperature at the current position and time; Phase balance distribution function at the current moment for: , in, Indicates the first j Volume fraction of the phase fluid, Indicates the first j The potential function value of the phase fluid, Indicates the first j Thermal diffusivity of the phase fluid Indicates the relationship with the first j Relaxation time related to the diffusivity of the phase fluid Represents a constant; Density distribution function at the previous moment for: , in, Indicates the time step. Indicates the previous moment; Temperature distribution function at the previous moment for: , Phase distribution function at the previous moment for: , Density distribution function The calculation formula is: , in, Indicates the relationship with the first jRelaxation time related to the kinematic viscosity of the phase fluid; Temperature distribution function The calculation formula is: , in, Indicates the relationship with the first j Relaxation time related to the thermal diffusivity of the phase fluid; Phase distribution function The calculation formula is: , In step S3, the specific calculation process for the potential function value of the three-phase fluid is as follows: The simplified formula for body free energy is: , in, Indicates the first j Volume fraction of the phase fluid; Indicates the first j Surface tension of the phase fluid; , Where W represents the interface thickness; The weighted average surface tension of a three-phase fluid is calculated using the following formula: , Indicates the first j The functional relationship satisfied by the volume fraction of the phase fluid is expressed by the following formula: , Surface tension between two phase interfaces The calculation formula is: , Volume forces throughout the computational domain The calculation formula is: , in, Indicates high density, Indicates low density, It represents the acceleration due to gravity.

[0013] In step S4, the evolution equation of TPLBFS is: , in, Indicates total flux; This represents the flux parallel to the x-axis; This represents the flux parallel to the y-axis; Indicates the volume of an adjustable control element; This represents the enclosed area of ​​the k-th grid cell; Indicates the source term; , Represents a unit vector related to spatial direction; Total flux The calculation formula is: , in, Indicates density, Represents the vector velocity in the x-direction. Represents the vector velocity in the y-direction. This indicates the volume fraction of the first phase in the mixture. Indicates the volume fraction of the second phase in the mixture; Flux parallel to the x-axis The calculation formula is: , in, The discrete velocity of the crystal lattice is represented in x Components along the axial direction; The discrete velocity of the crystal lattice is represented in y Components along the axial direction; Flux parallel to the y-axis The calculation formula is: , The calculation formula is: , in, This represents the component of surface tension along the x-axis. This represents the component of surface tension along the y-axis. This represents the component of the volume force along the x-axis. This represents the component of the volume force along the y-axis. Phase transition source term The calculation formula is: , in, Indicates thermal conductivity. Indicates ambient temperature. Indicates the thermal diffusivity. Indicates phase volume fraction; The phase field equations are modified as follows: , Where M represents thermal diffusivity; The temperature equation is modified as follows: , in, and Both represent constants determined by latent heat.

[0014] In step S5, the fourth-order Runge-Kutta method is used to calculate the equation. , Obtain total flux ; pressure The calculation formula is: fluid velocity The calculation formula is: Through the first j Volume fraction of phase fluid The calculation formula yields the volume fraction of the first phase. Second phase volume fraction : temperature The calculation formula is:

[0015] The first phase volume fraction obtained by calculation Second phase volume fraction The volume fraction of the third phase is determined using the following formula. : , density The calculation formula is: , Viscosity The calculation formula is: , This completes one iteration process.

[0016] Furthermore, in step S6, the k output in step S5 is used... rH2 Spreading coefficient S, future operating condition prediction results, and preset optimal operating window threshold k rH2 =0.35-0.45, the downhole actuator parameter range, and the hydrogen fingering leading edge safety distance threshold are used as inputs. If there is a deviation, a corresponding control command is generated and sent to the actuator to complete closed-loop control. Based on the judgment result, a quantitative control command is generated according to the preset trigger logic. The specific trigger logic and control action are fixed as follows: Category 1: Injection speed optimization, when krH2 When <0.25, an injection rate increase command is generated, increasing the injection rate from 10⁻⁴ m / s to 10⁻³ m / s; The second type: lubrication effect induced, when the spreading coefficient S>0 and the long-chain hydrocarbon saturation <5%, generates a long-chain hydrocarbon injection command, adding C 16 -C 20 The saturation of long-chain hydrocarbons is maintained at 8-12%; The third type is pressure limiting protection. When the depth of hydrogen fingering at the leading edge exceeds 80% of the safe distance, a pressure reduction command is generated to reduce the injection pressure by 10%-15%, or a seawater reinjection start command is generated. Based on the pre-simulation of the control command execution, verify that the operating conditions after control can return to the optimal operating window without any new safety risks. If the conditions are not met, the control command is re-optimized. The finally confirmed control command is issued to the downhole production pump set, injection device, and reinjection device for execution, and the command issuance and execution feedback data are recorded simultaneously. This step outputs an operating condition deviation diagnosis report, a quantitative closed-loop control command file, command execution feedback records, and a pre-evaluation result of the control effect.

[0017] This invention also provides a deep-sea hydrogen production and storage intelligent collaborative control system based on TPLBFS, used to execute the above-mentioned control method. The system includes a sensing layer, an intelligent collaborative computing engine layer, and a feedback control layer connected in sequence. The sensing layer includes a multi-source sensor array deployed in the deep-sea operation area for real-time acquisition of measured operating data during the production and storage process, and transmission to the intelligent collaborative computing engine layer. The intelligent collaborative computing engine layer has a built-in improved three-phase lattice Boltzmann flux solver TPLBFS, used to receive measured data to complete dynamic model updates, perform flow field numerical simulation and future operating condition prediction, and output performance indicators to the feedback control layer. The feedback control layer has a built-in model prediction control module, used to generate control commands based on predicted performance indicators, and send them to the downhole production pump group, injection device, and reinjection device for execution, completing the closed-loop intelligent collaborative control of the production and storage process.

[0018] Compared with the prior art, the present invention has the following advantages: High-precision and stable simulations under extreme conditions were achieved. By introducing the Lagrange multiplier mass correction source term, the mass conservation error of long-period simulations was reduced from 3-5% in traditional methods to below 0.01%, an improvement of 300-500 times. Through five specific optimizations targeting the extreme properties of hydrogen, the prediction error of TPLBFS in hydrogen systems was reduced from >50% in general methods to below 5%, and the spurious interface velocity was suppressed from 35% to below 3%, solving the problem of numerical instability of traditional methods in systems with high density and low viscosity ratios.

[0019] A complete real-time closed-loop control architecture is constructed, 2300 times of computational acceleration is achieved through GPU parallel acceleration, and the total system response time is less than 0.5 seconds, which realizes the leap from static offline simulation to real-time digital twin. The system can dynamically correct the model and generate regulatory instructions in a predictive manner based on real-time monitoring data. Compared with traditional manual decision-making, the reaction speed is increased by more than 10000 times, which can effectively avoid safety risks such as hydrogen escape and fingering propagation.

[0020] The long-chain hydrocarbon spreading effect is discovered for the first time and applied in engineering, breaking the technical prejudice that "impurity phases necessarily hinder flow". Under the conditions of 5-15% long-chain hydrocarbon saturation and spreading coefficient S>0, a lubricating film is formed by the spontaneous spreading of long-chain hydrocarbons at the hydrogen-seawater interface, which increases the hydrogen transmission efficiency by 40%, and significantly improves the hydrogen production recovery factor and the storage capacity utilization rate of the hydrogen storage reservoir.

[0021] It has significant engineering and economic value. Through accurate prediction and dynamic optimization, the cost of deep-sea hydrogen production can be reduced by 15-25%, and the utilization rate of hydrogen storage capacity can be increased by more than 20%; in the on-site implementation in the Kunlun hydrothermal area of the South China Sea, the daily hydrogen production of a single well is increased from 120kg to 165kg, with an increase of 38%, which has extremely high engineering promotion value. Brief Description of the Drawings

[0022] Figure 1 is an architecture logic diagram of the deep-sea hydrogen production and storage intelligent cooperative control system according to the present invention; Figure 2 is a schematic structural view of the matching stepped hydrogen capture cover of the present invention; Figure 3 is a schematic diagram of the component control logic of the matching stepped hydrogen capture cover of the present invention; Figure 4 is an enlarged schematic diagram showing the details of the three-phase interface in the matching stepped hydrogen capture cover of the present invention; Reference signs in the drawings: 1 - Oil injection pump, 2 - Water injection pump, 3 - Control valve, 4 - Stepped hydrogen capture cover, 5 - Heating and heat insulating layer, 6 - Capture cavity, 7 - Thermoelectric compensation module. Detailed Description of the Embodiments

[0023] The technical solutions of the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0024] The core of the TPLBFS-based deep-sea hydrogen production and storage intelligent cooperative control system and method disclosed by the present invention is to realize high-precision and stable simulation of the deep-sea hydrogen-seawater-long-chain hydrocarbon three-phase flow through the TPLBFS solver optimized for the extreme physical properties of hydrogen, and at the same time construct a five-dimensional closed-loop control architecture to realize real-time intelligent cooperative control of the production and storage process.

[0025] Embodiment 1 This embodiment uses a deep-sea natural hydrogen extraction project in the Kunlun hydrothermal vent area (118°45'E, 15°23'N) in the South China Sea as a case study, and fully implements the control method and system of this invention. The water depth in this area is 3050 meters, the hydrothermal vent temperature is 20°C, the hydrogen concentration is approximately 6.3 mmol / kg, and associated with the growth of chain hydrocarbons (C...). 16 -C 20 The saturation is approximately 10%.

[0026] Step 1: Physical model establishment and initialization; Based on the actual geological exploration data of the target hydrothermal zone, a three-dimensional porous media physical model was established, with the computational domain size determined to be 500m × 300m × 200m. A structured mesh was used for mesh generation, with a minimum mesh size Δx. min =0.5m, total grid number 1.2×10 7 The computational domain boundary conditions were set as follows: the upper boundary was a caprock closed boundary, the left, right, front, and rear boundaries were periodic boundaries, the lower boundary was the production well boundary, and Navier slip boundary conditions were set for both the well wall and the formation wall. The hydrogen phase slip length was set to 500 nm. The physical parameters of the grid cells were initialized based on well logging data: formation porosity 0.22, permeability 120 mD, initial pressure 30.5 MPa, temperature 4 °C, hydrogen phase volume fraction 0.06, seawater phase volume fraction 0.84, and long-chain hydrocarbon phase volume fraction 0.10. The density, viscosity, and surface tension of each phase were also initialized simultaneously. Linear interpolation was used to initialize the physical parameters of all grid cell interfaces to ensure a smooth transition of physical quantities between cells.

[0027] Step 2: Initialize and calculate the distribution function; Based on the improved TPLBFS solver, the density distribution function of the initial time step is calculated. Temperature distribution function Phase distribution function The density distribution function of the current node Phase distribution function and temperature distribution function It is based on the equilibrium distribution function at the current moment. , , and the distribution function of the previous time step , , We can approximate the result.

[0028] and The main descriptions are the density distribution and temperature distribution of the fluid, which are derived based on the fluid velocity u and corresponding physical quantities such as pressure and temperature. It describes the distribution of phase field, chemical potential, etc., and is used for capturing three-phase interfaces and obtaining three-phase distribution.

[0029] Density equilibrium distribution function at the current moment for: , The above function represents the distribution function of particles along direction i at position x and time t when local thermodynamic equilibrium is reached. Indicates the current position. Indicates the current moment. Indicates the weighting factor. Indicates pressure, This represents the mass density of a fluid at a specific location and time. Represents the speed of sound in a lattice. This represents the discrete velocity of the crystal lattice. It indicates the velocity of a fluid at a specific location and time.

[0030] Temperature equilibrium distribution function at the current moment for: , Where T represents the macroscopic thermodynamic temperature at the current position and time.

[0031] Phase balance distribution function at the current moment for: , in, Indicates the first j Volume fraction of the phase fluid, Indicates the first j The potential function value of the phase fluid, Indicates the first j Thermal diffusivity of the phase fluid Indicates the relationship with the first j Relaxation time related to the diffusivity of the phase fluid It is a constant, and is related to j A phase fluid-related factor used to correct or adjust the physical effects of the distribution function.

[0032] Density distribution function at the previous moment for: , in, Indicates the time step. It indicates the previous moment.

[0033] Temperature distribution function at the previous moment for: , Phase distribution function at the previous moment for: , Therefore, the density distribution function The calculation formula is: , in, Indicates the relationship with the first j Relaxation time related to the kinematic viscosity of the phase fluid.

[0034] Temperature distribution function The calculation formula is: , in, Indicates the relationship with the first j Relaxation time related to the thermal diffusivity of the phase fluid.

[0035] Phase distribution function The calculation formula is: , The lattice model adopts the D3Q19 three-dimensional discrete velocity model, with weighting factors... Discrete velocity According to the D3Q19 model standard, the lattice sound velocity c s =c / √3, where c is the unit velocity per grid cell. The specific calculation formula for the distribution function is shown in claim 3, where the flow field relaxation time... Temperature field relaxation time Phase relaxation time The value is determined based on the kinematic viscosity, thermal diffusivity, and diffusion coefficient of each phase. The calculation formula is as follows: The essential difference between the improved three-phase lattice Boltzmann flux solver TPLBFS and the traditional LBFS is as follows: (1) The lattice Boltzmann distribution function is transformed into a macroscopic flux expression of the finite volume method through Chapman-Enskog multi-scale expansion, realizing a direct mapping from the mesoscopic distribution function to the macroscopic flux, avoiding the numerical error introduced by the additional interpolation required by the traditional LBFS when solving the macroscopic flux; (2) For the three-phase fluid system, the Navier-Stokes equation, Cahn-Hilliard phase field equation and energy equation are solved simultaneously, and the mass correction source term is introduced through the Lagrange multiplier method to ensure that the volume fraction of the three phases strictly meets the conservation constraints at each time step; (3) The complete TPLBFS algorithm flow is: distribution function initialization → collision step (BGK approximation) → migration step → flux calculation (reconstructing macroscopic flux from the distribution function) → macroscopic quantity update (density, velocity, temperature, volume fraction of each phase) → mass conservation correction → enter the next time step.

[0036] The TPLBFS method reduces the mass conservation error in long-period simulations from 3%-5% in traditional methods to below 0.01% by introducing a mass correction source term using the Lagrange multiplier λ into the phase-field equations, achieving an improvement of 300-500 times. This correction method determines the value of λ by solving the mass conservation equations at each time step. Specifically, after the collision-migration step at each time step, the deviation of the sum of the current phase volume fractions from 1 is calculated. This deviation is used as a constraint, and the optimal value of λ is obtained by iteratively solving the linear equations, thereby correcting the distribution functions of each phase and ensuring that the sum of the three phase volume fractions is always 1.

[0037] , in, Let be the kinematic viscosity of the j-th phase fluid. Here is the thermal diffusivity, and the initial time step is... d t Take 0.001s.

[0038] Step 3: Calculation of chemical potential and forces; In multiphase flow simulations, in order to accurately describe the interface information of two-phase flows, in addition to the flow equations, interface trapping equations are also needed. In diffusion interface methods, the Cahn-Hilliard (CH) equations are a widely used and important tool for calculating the trapping phase interface.

[0039] The formula for volume free energy can be simplified to the following form: , in, Indicates the first j Volume fraction of the phase fluid; Indicates the first j Surface tension of the phase fluid.

[0040] , Where W represents the interface thickness. The weighted average surface tension of a three-phase fluid is calculated using the following formula: , in, Indicates the first in the mixture j The functional relationship satisfied by the volume fraction of the phase fluid is expressed by the following formula: , Surface tension between two phase interfaces The calculation formula is: , Volume forces throughout the computational domain The calculation formula is: , in, Indicates high density, Indicates low density, It represents the acceleration due to gravity.

[0041] When one of the phases disappears and its order parameter When it is 0, such as =0, the phase field equations of the three phases become the classic two-phase CH equations, satisfying .

[0042] Calculate the chemical potentials of the hydrogen phase, seawater phase, and long-chain hydrocarbon phase. m 1. m 2. m 3. The interface thickness W is taken as 5 times the grid size, and the weighted average surface tension δ is calculated based on the measured values ​​of the three-phase interfacial tension. s HW =0.072 N / m, s HO =0.025N / m, s OW =0.036 N / m. The interfacial surface tension F is calculated based on chemical potential. s Calculate the gravitational body force F based on the density difference between the two phases b And substitute it into the source term vector Q.

[0043] Step 4: Real-time data acquisition and model assimilation; Field-measured data were collected using a multi-source sensor array in the sensing layer: three in-situ Raman spectroscopy probes were deployed at the hydrothermal vent outlet to measure hydrogen concentration and long-chain hydrocarbon mass fraction with a measurement accuracy of 0.1 mmol / kg; eight MEMS pressure and temperature sensors were deployed every 50 meters along the extraction pipeline, with a pressure measurement accuracy of 0.1%FS; Coriolis flow meters were installed on the main pipeline with a flow measurement accuracy of 0.2%. All sensors collected data at a frequency of 1 Hz and transmitted it in real time to the intelligent collaborative computing engine layer of the surface control platform via submarine fiber optic cable. The data assimilation module used an ensemble Kalman filter (EnKF) to assimilate the measured pressure, temperature, component concentration, and flow data into the TPLBFS model, dynamically updating the model's boundary conditions, porosity, permeability, and other physical properties to eliminate deviations between the initial geological model and the actual strata, ensuring consistency between the simulation results and the physical entity.

[0044] The EnKF-TPLBFS coupling method is as follows: (1) Mapping of measured data to grid cells: Kriging interpolation is used to map the sparsely distributed sensor measured data (pressure, temperature, component concentration) to all grid cells in the three-dimensional computational domain through spatial statistical interpolation, generating an assimilation dataset that matches the TPLBFS grid; (2) Assimilation time step matching: EnKF data assimilation is performed once every 10 TPLBFS computation time steps (corresponding to 0.01 seconds of actual time) to ensure that the assimilation frequency matches the time scale of the flow field evolution characteristics; (3) Model parameter update method: After assimilation, based on the deviation between the measured data and the simulated predicted values, the porosity φ, permeability K and other physical property parameters of the formation are dynamically adjusted through the EnKF update equation. Specifically, the relaxation time τ of each phase fluid is adjusted to achieve continuous update of physical property parameters. The relaxation time adjustment formula is: Where α is the relaxation factor, with a value of 0.3-0.7; (4) The engineering significance of the data assimilation error threshold being less than 5% is that the threshold corresponds to the acceptable range of engineering errors in deep-sea mining and storage operations. When the prediction error exceeds 5%, the confidence level of the hydrogen leakage risk assessment will be less than 95%, and assimilation needs to be re-executed to ensure the reliability of safety decisions.

[0045] Step 5: TPLBFS solution and load condition prediction; To address the extreme physical properties of deep-sea hydrogen, five specific optimizations were performed on TPLBFS, as described in claim 6, wherein: Dynamic state equation coupling: At each time step, the pressure P and temperature T measured by the sensor are read, and the compressibility factor Z is solved iteratively through the modified Peng-Robinson state equation to dynamically update the hydrogen phase density. This addresses the issue of density prediction errors exceeding 200% under high pressure in the deep sea, as determined by the ideal gas law. Adaptive time step control: Each time step is based on the maximum velocity of the flow field. Dynamically update time step To ensure the accuracy of interface capture when hydrogen flows at high speed and to avoid numerical explosion; Fine control of interface thickness: The diffusion interface thickness ε is set to 10 times the grid size to suppress the numerical oscillation of the interface under a high density ratio of 1:13.75; Dynamic evaluation of spreading coefficient: The spreading coefficient is calculated every 100 time steps. In the initial operating conditions of this embodiment =+0.011N / m>0, which satisfies the lubrication effect activation condition; Mass conservation correction: Navier slip boundary conditions, with the wall boundary using the slip formula: , As the source term for mass correction, the corrected phase-field equation is: , in The mass conservation equations are determined by solving them at each time step. The mass conservation error is controlled to be below 0.01%.

[0046] By using the optimized TPLBFS computational interface to calculate macroscopic flux and source terms, the evolution equations are solved using the fourth-order Runge-Kutta method to obtain macroscopic physical parameters such as pressure field, velocity field, and volume fraction distribution of each phase at the current time step. The calculation then evolves forward 10 seconds to predict the hydrogen phase distribution, fingering expansion trend, and dynamic relative permeability changes under future operating conditions. The computational engine utilizes an NVIDIA A10080GB GPU for parallel acceleration, combined with an intrinsic orthogonal decomposition (POD) order reduction model, achieving a computation speed of 0.38 seconds per time step, a speedup exceeding 2300 times, meeting real-time requirements.

[0047] The hardware architecture of the intelligent collaborative computing engine layer is as follows: it adopts an NVIDIA A10080GB GPU cluster, and achieves large-scale parallel computing through the Domain Decomposition Method, with each GPU responsible for approximately 1×10 6 The computational tasks for each grid cell are performed using high-speed interconnection between GPUs via NVLink, achieving a speedup of 2300 times and a total system response time of less than 0.5 seconds. The software modules include: a data preprocessing module (responsible for cleaning, filtering, and standardizing sensor data), a TPLBFS core solver module (responsible for flow field numerical simulation), an EnKF data assimilation module (responsible for dynamic model updates), a model predictive control module (responsible for condition prediction and control decisions), and a results output and visualization module. The sensing layer transmits sensor data via submarine fiber optic cable using the Modbus protocol. The computation engine layer and the feedback control layer communicate bidirectionally via Ethernet TCP / IP protocol. Control commands are issued from the feedback control layer to the downhole actuators via submarine fiber optic cable.

[0048] Step 6: Closed-loop control decision-making and execution; Based on the prediction results, the dynamic relative permeability k of the hydrogen phase was calculated. rH2 In the initial operating condition of this embodiment, k rH2 =0.25, lower than the optimal operating window of 0.35-0.45, indicating that hydrogen exists in the form of isolated bubbles, with a phase trapping fraction of 35%, resulting in a low recovery rate. The model predictive control module of the feedback control layer generates a dual regulation strategy based on the prediction results: The optimal operating window threshold k rH2The determination of 0.35-0.45 is based on: 100 sets of different operating conditions (injection pressure 20MPa-30MPa, injection speed 0.5×10⁻⁶). -4 m / s - 2.0 × 10 -3 TPLBFS simulation experiments (m / s, long-chain hydrocarbon saturation 0-15%) revealed that the hydrogen transport efficiency was highest and the fingering risk was below 1% when the dynamic relative permeability of the hydrogen phase, krH2, was in the range of 0.35-0.45. Below 0.25, the phase trapping fraction exceeded 30%, and above 0.45, the fingering risk increased significantly. The logic for determining the spreading coefficient calculation cycle of "every 100 time steps" is as follows: 100 TPLBFS time steps correspond to 0.1 seconds of actual time (time step δt=0.001s). This cycle can capture the dynamic changes in the spreading state of long-chain hydrocarbons in a timely manner (the characteristic time of interfacial tension change is about 0.05-0.2 seconds) while ensuring computational efficiency and not significantly increasing additional computational overhead. The quantitative relationship between the optimal saturation of long-chain hydrocarbons (8-12%) and transport efficiency is as follows: When the saturation of long-chain hydrocarbons is 8-12%, the interfacial lubricating film thickness is 8-12 μm, and the hydrogen transport efficiency is 40% higher than that without long-chain hydrocarbons. When the saturation is below 5%, the lubricating film is incomplete, and the improvement is less than 15%. Above 12%, the lubricating film thickness tends to saturate, and the transport efficiency no longer increases significantly. Moreover, excessively high long-chain hydrocarbon saturation will occupy effective pore space and reduce hydrogen storage capacity. The engineering significance of a data assimilation error threshold of less than 5% is that this threshold corresponds to the acceptable engineering error range in deep-sea mining and storage operations. When the prediction error exceeds 5%, the confidence level of the hydrogen leakage risk assessment is less than 95%, and assimilation needs to be re-performed to ensure the reliability of safety decisions.

[0049] Injection speed optimization: The seawater injection speed was increased from 10... -4 m / s increased to 8.5×10 -4 m / s, by breaking up isolated hydrogen bubbles through viscous forces, promoting bubble coalescence to form a continuous gas phase channel; Lubrication effect induced: C is supplemented by injection through downhole injection device. 16 -C 20 For long-chain hydrocarbons, maintaining their saturation within the optimal window of 8-12% maximizes the lubricating effect of the spreading effect and reduces interfacial flow resistance.

[0050] Control commands are transmitted in real time to the downhole mining pump unit and injection device via submarine optical cable. After execution, new measured data are collected through sensor array, and the next closed loop of "sensing-assimilation-prediction-decision-execution" is initiated.

[0051] Implementation effect verification Forty-eight hours after the implementation of the above control strategy, on-site measurement data showed that the dynamic relative permeability of hydrogen, k... rH2The hydrogen recovery rate stabilized in the range of 0.33-0.37, with an error of less than 5% compared to the TPLBFS prediction. The daily hydrogen production per well increased from 120 kg to 165 kg, an increase of 38%. The phase capture fraction decreased from 35% to 14%, and the hydrogen recovery rate was significantly improved, verifying the technical effectiveness of the present invention.

[0052] Example 2 This embodiment uses a hydrogen storage project in a depleted subsea oil reservoir in the South China Sea as a case study to implement the system and method of this invention and verify its application effect in hydrogen storage safety assessment and risk control. The depleted oil reservoir has a reservoir size of 2km × 1.5km × 200m, an average porosity of 0.18, an average permeability of 50mD, a designed hydrogen injection capacity of 1 million standard cubic meters, an injection pressure of 26MPa, a caprock fracture pressure of 30MPa, and a residual long-chain hydrocarbon saturation of approximately 5%.

[0053] Following the steps in Example 1, the physical model was established, initialized, the TPLBFS solver was deployed, and the sensor array was installed. The system first ran TPLBFS to simulate the hydrogen injection process on a one-year timescale based on the initial seismic interpretation data. The baseline prediction results showed that under an injection pressure of 26 MPa, the hydrogen gas was distributed in a stable layered manner with no fingering risk and a hydrogen leakage probability of <1%.

[0054] After hydrogen injection began, the system collected wellhead pressure and flow rate data in real time at a frequency of 1Hz. In the third month of hydrogen injection, measured data showed that the bottom hole pressure rose to 27.8 MPa, exceeding the design value of 26 MPa, and the injection rate reached 5000 m³ / h. 3 / day. After the EnKF data assimilation module assimilates the measured data into the TPLBFS model, it identifies an unexplored high-permeability zone in the northwest corner of the reservoir, with a permeability as high as 200mD, which is 4 times the average permeability. Hydrogen exhibits severe viscous fingering at this location.

[0055] TPLBFS forward evolution simulation early warning: If injection continues at the current pressure, the fingering front will reach the capping layer in 45 days, increasing the probability of hydrogen leakage to 13%, posing a serious safety risk. The model predictive control module generates a control scheme after multi-objective optimization (maximizing hydrogen storage capacity vs. minimizing leakage risk). Reduce the injection pressure to 25MPa and reduce the finger extension driving force; A new monitoring well was added near the well location in the high-permeability zone to monitor changes in the pressure field in real time. The periodic "inject-stop-test" operation mode is activated, with 2 days of injection and 1 day of stoppage. During the injection and stoppage period, the pressure field diffusion is monitored to stabilize the gas-liquid interface.

[0056] After 30 days of operation using the new scheme, TPLBFS simulation results showed that the hydrogen fingering front stabilized at 50 meters from the capping layer, the leakage risk was reduced to <2%, and an 82% hydrogen recovery rate was ultimately achieved. The hydrogen storage capacity utilization rate increased by 22%, and there were zero leakage incidents throughout the entire lifecycle. Compared to the traditional quarterly offline simulation and manual decision-making model, the reaction speed of this invention is increased by more than 10,000 times, significantly improving the safety of subsea hydrogen storage.

[0057] Example 3 This embodiment is a supporting solution for Embodiment 1. The system is equipped with a stepped hydrogen capture hood 4, the structure of which is as follows: Figure 2 As shown, the main body of the trapping hood is a three-stage stepped expansion trapping chamber 6. The inner wall of the chamber is provided with a heating insulation layer 5, and a thermoelectric compensation module 7 is installed at the bottom. The thermoelectric compensation module is communicatively connected to the feedback control layer 8.

[0058] When the intelligent collaborative computing engine layer predicts that the seawater ambient temperature is below 10℃, which leads to increased viscosity of long-chain hydrocarbons and suppressed spreading effect, the feedback control layer automatically adjusts the output power of the thermoelectric compensation module to raise the local ambient temperature inside the capture hood to the optimal lubrication temperature range of 15-20℃, ensuring that long-chain hydrocarbons are always in a low viscosity state, maintaining a stable interface spreading effect, and ensuring hydrogen transmission efficiency.

[0059] The above description is merely a preferred embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for intelligent collaborative control of deep-sea hydrogen extraction and storage based on TPLBFS, characterized in that, Includes the following steps: S1. Physical Model Establishment and Initialization: Using seismic exploration data, well logging data, wellbore parameters, formation static physical properties, hydrogen-seawater-long-chain hydrocarbon three-phase fluid reference physical properties, and numerical simulation grid configuration parameters as the sole inputs, perform the following operations: 1.1 Based on the actual geological characteristics of the target work area, a three-dimensional porous media physical model was established, the computational domain size was determined, a three-dimensional structured mesh was used to divide the computational domain, and the minimum mesh size Δx was defined. min ; 1.2 Fixed computational domain boundary conditions: The upper boundary is set as a caprock closed boundary, the left, right, front, and rear boundaries are set as periodic boundaries, the lower boundary is set as the production well / injection well boundary, and Navier slip boundary conditions are set for both the well wall and the formation wall. 1.3 Initialize the physical parameters of the grid cells in the entire computational domain, including formation porosity, permeability, initial pressure, initial temperature, volume fraction, density, viscosity and surface tension of hydrogen phase / seawater phase / long-chain hydrocarbon phase, and force the constraint that the sum of the volume fractions of the three phases is 1. 1.4 Using the initial physical parameters of the center points of the grid cells as a reference, bilinear interpolation is used for the two-dimensional model and trilinear interpolation is used for the three-dimensional model to complete the initialization of the physical parameters of the interface between adjacent grid cells. The interpolation method is based on the physical parameter values ​​of the center points of four (two-dimensional) / eight (three-dimensional) adjacent grid cells, and the parameter values ​​on the interface are calculated by linear weighting to ensure a smooth transition of physical quantities between cells. This step outputs the three-dimensional physical model file, the structured mesh file, the boundary condition configuration table, and the initial physical parameter matrix of the grid cells and interfaces in the entire computational domain, which serve as the input reference for all subsequent steps. S2. Based on the improved three-phase lattice Boltzmann flux solver TPLBFS, using the initial physical parameter matrix of the mesh elements and interfaces output in step S1, the standard parameters of the D3Q19 three-dimensional discrete velocity lattice model, and the initial time step... δ t The density distribution function at the current time step is calculated using the kinematic viscosity, thermal diffusivity, and diffusion coefficient of each phase fluid as input. f i ∗ Temperature distribution function g i ∗ Phase distribution function h i ∗j ; ， The TPLBFS is coupled with the Navier-Stokes equation, the Cahn-Hilliard phase field equation, and the energy equation, and five specific optimizations are performed for the extreme properties of hydrogen. S3, the phase distribution function matrix and three-phase volume fraction output in step S2. ϕ j Matrix, measured value of three-phase interfacial tension σ HW / σ HO / σ OW Using interface thickness parameter W, gravitational acceleration g, and density parameters of each phase fluid as inputs, the chemical potential function of the three-phase fluid is calculated, and the potential function values ​​of the hydrogen phase, seawater phase, and long-chain hydrocarbon phase are obtained. μ 1. μ 2. μ 3. Simultaneously calculate surface tension and volume force. This step outputs the three-phase chemical potential matrix, surface tension field matrix, volume force field matrix, and evolution equation source term vector Q matrix of the entire computational domain. S4. Measured data acquisition and dynamic model assimilation: Taking the configuration parameters of the multi-source sensor array in the deep-sea operation area, the initial physical model output in step S1, the current physical field matrix output in step S3, and the configuration parameters of the ensemble Kalman filter (EnKF) algorithm as input, perform the following operations: 4.1 Start the multi-source sensor array and collect on-site measured data at a fixed frequency of 1Hz; the sensor array includes an in-situ Raman spectroscopy sensor, a MEMS pressure sensor, a temperature sensor, and a Coriolis flow meter, and the collected data items are fixed as hydrogen concentration, long-chain hydrocarbon mass fraction, pressure, temperature, and flow rate; 4.2 Perform outlier removal, filtering and smoothing, and timestamp alignment preprocessing on the collected raw data to generate a standardized measured dataset; 4.3 The EnKF algorithm of ensemble Kalman filtering is adopted to assimilate the standardized measured dataset into the TPLBFS solver and dynamically update the boundary conditions, formation porosity, and permeability parameters of the model. 4.4 Verify that the relative error between the model predictions and the field measurements after assimilation is less than 5%. If the error is not met, the assimilation process should be repeated. This step outputs a standardized measured dataset with timestamp alignment, the updated model boundary condition configuration file, the formation physical parameter matrix, and a data assimilation error assessment report. S5, TPLBFS Solution and Operating Condition Prediction: Taking the model parameters and boundary conditions assimilated and updated in step S4, the three-phase distribution function matrix output in step S2, the source term vector Q matrix output in step S3, the preset five specific optimization parameters, the fourth-order Runge-Kutta solver configuration, and the prediction duration parameters as input, the evolution equation is solved using the fourth-order Runge-Kutta method to obtain the macroscopic physical parameters of the current time step. The system then evolves forward to predict the phase distribution, velocity field, and pressure field under future operating conditions, outputting the macroscopic physical field matrix of the full computational domain at the current time step, the future operating condition prediction dataset, and the dynamic relative permeability k of the hydrogen phase. rH2 Compare the calculation results of the spreading coefficient S with the numerical simulation accuracy verification report; S6. Closed-loop control decision-making and execution, using the output k from step S5. rH2 Spreading coefficient S, future operating condition prediction results, and preset optimal operating window threshold k rH2 =0.35-0.45, the downhole actuator parameter range, and the hydrogen fingering leading edge safety distance threshold are the inputs. If there is a deviation, a corresponding control command is generated and sent to the actuator to complete the closed-loop control. 6.1 Based on the judgment result, a quantitative control command is generated according to the preset trigger logic. The specific trigger logic and control action are fixed as follows: Category 1: Injection speed optimization, when k rH2 When <0.25, an injection rate increase command is generated, increasing the injection rate from 10⁻⁴ m / s to 10⁻³ m / s; The second type: lubrication effect induced, when the spreading coefficient S>0 and the long-chain hydrocarbon saturation <5%, generates a long-chain hydrocarbon injection command, adding C 16 -C 20 The saturation of long-chain hydrocarbons is maintained at 8-12%; The third type is pressure limiting protection. When the depth of hydrogen fingering at the leading edge exceeds 80% of the safe distance, a pressure reduction command is generated to reduce the injection pressure by 10%-15%, or a seawater reinjection start command is generated. 6.2 Based on the control commands, perform pre-simulation to verify that the operating conditions after control can return to the optimal operating window without any new safety risks. If the conditions are not met, the control commands should be re-optimized. 6.3 Issue the final confirmed control commands to the downhole production pump unit, injection unit, and reinjection unit for execution, and record the command issuance and execution feedback data simultaneously; this step outputs a working condition deviation diagnosis report, a quantitative closed-loop control command document, command execution feedback records, and a pre-evaluation result of the control effect; S7. Loop Iteration and Operation Termination: Using the measured data after execution output in step S6, the physical field parameters of the current time step, and the preset end-of-cycle judgment conditions of the mining and storage operation as inputs, perform the following operations: 7.1 Use the updated physical field parameters and model parameters of the current time step as the initial inputs for the next time step. 7.2 Determine whether the end point of the production and storage operation cycle has been reached. If not, repeat steps S2 to S6 to enter the next closed-loop cycle. If the design recovery rate, total hydrogen injection, or end point conditions such as the end of the operation cycle are reached, terminate the cycle. This step outputs the initial input parameter matrix for the next cycle, the full production and storage operation data set, and the operation termination summary report.

2. The intelligent collaborative control method for deep-sea hydrogen extraction and storage based on TPLBFS according to claim 1, characterized in that, In step S1, the physical model is established based on the actual geological characteristics of the deep-sea hydrothermal extraction area or the depleted oil reservoir hydrogen storage area. The upper boundary of the computational domain is the caprock closed boundary, the left and right boundaries are periodic boundaries, the lower boundary is the production well / injection well boundary, and the wall is set with Navier slip boundary conditions. The initialized physical parameters include pressure, temperature, volume fraction of each phase fluid, density and viscosity.

3. The intelligent collaborative control method for deep-sea hydrogen extraction and storage based on TPLBFS according to claim 1, characterized in that, In step S2, the density distribution function f i ∗ Temperature distribution function g i ∗ Phase distribution function h i ∗j The density equilibrium distribution function is calculated based on the equilibrium distribution function at the current moment and the distribution function at the previous moment. The specific calculation formula is: density equilibrium distribution function; where the constant Γ is taken as 1.0, which is obtained by fitting three-phase interfacial tension experimental data: ， in, Indicates the current position. Indicates the current moment. This indicates the lattice orientation of the particle at the current moment. Indicates the weighting factor. Indicates pressure, This represents the mass density of a fluid at a specific location and time. Represents the speed of sound in a lattice. This represents the discrete velocity of the crystal lattice. This indicates the velocity of a fluid at a specific location and time. Temperature equilibrium distribution function at the current moment for: ， Where T represents the macroscopic thermodynamic temperature at the current position and time; Phase balance distribution function at the current moment for: ， in, Indicates the first j Volume fraction of the phase fluid, Indicates the first j The potential function value of the phase fluid, Indicates the first j Thermal diffusivity of the phase fluid Indicates the relationship with the first j Relaxation time related to the diffusivity of the phase fluid Represents a constant; Iterative formula for distribution function: 。 4. The intelligent collaborative control method for deep-sea hydrogen extraction and storage based on TPLBFS according to claim 1, characterized in that, In step S3, the chemical potential function is calculated based on a simplified formula for volume free energy, specifically as follows: The intelligent collaborative control method for deep-sea hydrogen extraction and storage based on TPLBFS according to claim 1 is characterized in that, in step S3, the specific calculation process of the potential function value of the three-phase fluid is as follows: The simplified formula for body free energy is: ， in, Indicates the first j Volume fraction of the phase fluid; Indicates the first j Surface tension of the phase fluid; ， Where W represents the interface thickness; The weighted average surface tension of a three-phase fluid is calculated using the following formula: ， Indicates the first j The functional relationship satisfied by the volume fraction of the phase fluid is expressed by the following formula:

5. Surface tension between two phase interfaces The calculation formula is: ， Volume forces throughout the computational domain The calculation formula is: ， in, Indicates high density, Indicates low density, It represents the acceleration due to gravity.

6. The intelligent collaborative control method for deep-sea hydrogen extraction and storage based on TPLBFS according to claim 1, characterized in that, In step S4, the multi-source sensor array includes an in-situ Raman spectroscopy sensor, a MEMS pressure sensor, a temperature sensor, and a Coriolis flow meter. It collects measured data of hydrogen concentration, long-chain hydrocarbon mass fraction, pressure, temperature, and flow rate at a frequency of 1 Hz. The measured data are assimilated into the TPLBFS model through an ensemble Kalman filter (EnKF) to dynamically update the boundary conditions and physical property parameters of the model.

7. The intelligent collaborative control method for deep-sea hydrogen extraction and storage based on TPLBFS according to claim 1, characterized in that, In step S5, the evolution equation of the TPLBFS is: ， Where W is the total flux vector, W=( ) T , For the volume of the grid control unit, Let k be the area of ​​the k-th control surface. , , It is a unit vector in spatial direction. , , These represent the macroscopic fluxes along the x, y, and z axes, respectively, and Q is the source term vector. The TPLBFS performs five specific optimizations for extreme deep-sea hydrogen conditions, including: Optimization 1: Dynamic state equation coupling. At each time step, the pressure P and temperature T are measured by the sensor, and the Peng-Robinson state equation with temperature correction term is called to dynamically update the hydrogen phase density ρ. H2 The corrected equation of state is: , Where R is the gas constant, b is the gas volume correction parameter, a(T) is the gas property correlation parameter, ξ(T) is the temperature correction factor, with a value of 0.01-0.1; fc(ρ,T) is the density-temperature correction function, fitted with a cubic polynomial, fc(ρ,T)=aρ 3 +bρ 2 +cρ+d, the coefficient ad was obtained by fitting experimental data of hydrogen density in the range of 30MPa and 0-20℃; Optimization 2: Adaptive time step control, dynamically updating the calculation step size at each time step. The calculation formula is: The CFL number is 0.

1. Minimum grid size, This represents the maximum velocity of the flow field. Optimization 3: Fine control of interface thickness, with the diffusion interface thickness ε set to 8-12 times the mesh size to suppress interface numerical oscillations under high density ratios; Optimization 4: Dynamic evaluation of the spreading coefficient, calculating the spreading coefficient every 100 time steps. ,in The interfacial tension between hydrogen and seawater. For the interfacial tension between hydrogen and long-chain hydrocarbons, For long-chain hydrocarbon-seawater interfacial tension, when S>0, the lubrication effect-induced mode is activated; Optimization 5: Navier slip boundary conditions, with the wall boundary using the slip formula: , As the source term for mass correction, the corrected phase-field equation is: , in The mass conservation equations are determined by solving them at each time step. The mass conservation error is controlled to be below 0.01%.

8. The intelligent collaborative control method for deep-sea hydrogen extraction and storage based on TPLBFS according to claim 1, characterized in that, In step S6, the optimal operating window is the dynamic relative permeability k of the hydrogen phase. rH2 =0.35-0.45, based on the dynamic relative permeability and spreading coefficient, three types of control commands are generated, and the specific triggering logic and control actions are as follows: Category 1: Injection speed optimization, when k rH2 When <0.25, the injection rate is increased from 10⁻⁴ m / s to 10⁻³ m / s to break up isolated hydrogen bubbles through viscous force and improve gas connectivity; Category 2: Lubrication effect induced, when the spreading coefficient S>0 and the long-chain hydrocarbon saturation <5%, C is injected. 16 -C 20 Long-chain hydrocarbons, with their saturation maintained at 8-12%, spread at the interface to form a lubricating film, reducing interfacial flow resistance; The third type is pressure limiting protection. When the depth of the hydrogen fingering edge exceeds 80% of the safe distance, the injection pressure will be reduced by 10%-15%, or seawater reinjection will be initiated to suppress fingering expansion.

9. A deep-sea hydrogen extraction and storage intelligent collaborative control system based on TPLBFS, characterized in that, For executing the control method according to any one of claims 1-7, the system comprises a sensing layer, an intelligent collaborative computing engine layer, and a feedback control layer that are sequentially connected via submarine optical cables and Ethernet; the sensing layer includes a multi-source sensor array deployed in the deep-sea operating area for real-time acquisition of measured operating data during the production and storage process, and transmission to the intelligent collaborative computing engine layer; the intelligent collaborative computing engine layer has a built-in improved three-phase lattice Boltzmann flux solver TPLBFS for receiving measured data to complete dynamic model updates, performing flow field numerical simulation and future operating condition prediction, and outputting performance indicators to the feedback control layer; the feedback control layer has a built-in model prediction control module for generating control commands based on predicted performance indicators, issuing them to the downhole production pump group, injection device, and reinjection device for execution, and completing closed-loop intelligent collaborative control of the production and storage process; the intelligent collaborative computing engine layer adopts a GPU parallel computing architecture.

10. The intelligent collaborative control system for deep-sea hydrogen extraction and storage based on TPLBFS according to claim 8, characterized in that, The system is also equipped with a stepped hydrogen capture hood, which has a built-in heating and insulation layer and a thermoelectric compensation module. The thermoelectric compensation module is communicatively connected to the feedback control layer. When the intelligent collaborative computing engine layer predicts that the ambient temperature is below 10°C, the feedback control layer automatically adjusts the thermoelectric compensation power to maintain the local ambient temperature inside the capture hood in the optimal lubrication temperature range of 15-20°C. The power of the thermoelectric compensation module increases by 2% for every 1°C decrease, ensuring that long-chain hydrocarbons are always in a low viscosity state.

Images