Pressure balance control method, device and system for compressible multiphase flow

Abstract

This invention provides a pressure balance control method, apparatus, and system for compressible multiphase flow. The method includes calculating the flux and non-conservative terms contributing to adjacent cells through the interface based on the flow field states on both sides of the mesh interface, thereby updating the conserved quantities of each mesh cell. For each mesh cell, based on the characteristics of the state equations of each phase and the updated conserved quantities, an adaptive pressure solution path is selected to obtain an intermediate pressure value. Constraints based on the reference pressures of each phase are introduced onto the intermediate pressure value to output a corrected cell equilibrium pressure. The cell equilibrium pressure is used to update the thermodynamic state variables within the cell to advance the flow field calculation. This invention addresses the technical problem of global computational collapse caused by the lack of physical constraints in the pressure balance solution process under extreme conditions in existing numerical simulation methods for compressible multiphase flow, leading to the generation and diffusion of non-physical processes.

Description

Technical Field

[0001] This invention relates to the field of computational fluid dynamics, and more particularly to a method, apparatus, and system for pressure balance control of compressible multiphase flow. Background Technology

[0002] In the field of numerical simulation of compressible multiphase flows, especially in extreme engineering scenarios involving extremely high pressure-density ratios such as underwater explosions and high-speed impacts, the five-equation model based on the assumption of velocity-pressure balance has become the mainstream theoretical framework. This model significantly reduces computational complexity while maintaining physical plausibility by introducing phase fractional equations to characterize interface evolution. However, applying this model to stable numerical calculations under realistic three-dimensional complex conditions still faces fundamental challenges. The core contradiction lies in the fact that extreme physical conditions easily induce minute errors at the numerical level, and existing mainstream technical solutions typically consist of two parts: first, a three-dimensional flux solver based on coordinate system rotation, which introduces significant computational overhead on unstructured meshes; second, a fixed, non-adaptive pressure solution strategy. This combined approach lacks effective fault tolerance and correction mechanisms when dealing with extreme physical states.

[0003] This systemic vulnerability stems from and is most prominently manifested in the crucial step of pressure equilibrium solution. Existing technological bottlenecks are concentrated in this critical step; both purely analytical and purely iterative methods are purely mathematical solutions, outputting only numerical solutions without establishing a mandatory connection to physical reality. When faced with extreme conditions such as strong shock waves and large deformations, these methods either fail due to limitations imposed by the state equations or diverge due to sensitivity to initial conditions, easily generating non-physical solutions such as negative pressure and negative density. More seriously, once such non-physical solutions occur, they rapidly spread throughout the computational domain through flux calculations, forming a chain reaction from single-point failure to global computational collapse, causing simulation interruptions.

[0004] Therefore, how to design a pressure balance solution mechanism with inherent fault tolerance under extreme conditions that can ensure the production of physically feasible solutions, thereby guaranteeing the global robustness of the entire multiphase flow simulation system, has become a core and urgent technical problem restricting the application of this technology in engineering. Summary of the Invention

[0005] The purpose of this invention is to provide a pressure balance control method, device and system for compressible multiphase flow, so as to solve the technical problem that the lack of physical constraints in the pressure balance solution process of existing compressible multiphase flow numerical simulation methods under extreme conditions leads to the generation and diffusion of non-physical solutions and causes global calculation collapse.

[0006] To achieve the above objectives, the present invention is implemented as follows: In a first aspect, the present invention provides a pressure balance control method for compressible multiphase flow, the method comprising: The flux and non-conservative terms at the interface are calculated based on the flow field state on both sides of the mesh interface to contribute to adjacent cells, so as to update the conserved quantities of each mesh cell. For each grid cell, the pressure solution path is adaptively selected based on the state equation characteristics of each phase within it and the updated conservation quantities to obtain the intermediate pressure value. Constraints based on the reference pressure of each phase are introduced on the intermediate pressure value to output the corrected cell equilibrium pressure. The thermodynamic state variables within the unit are updated using the unit's equilibrium pressure to advance the flow field calculation.

[0007] Furthermore, the compressible multiphase flow includes a phase described by the ideal gas law and a phase described by the rigid gas law. This defines the boundary of extreme and common physical problems, such as underwater explosions (high-pressure air / water) and fuel injection (fuel vapor / liquid), the core of which is the significant difference in compressibility and pressure threshold between the phases. Ideal gas phases (such as air) have no lower pressure limit, while rigid gas phases (such as water) have a definite and very large negative infinity pressure parameter. The adaptive pressure solution path proposed in this invention is highly effective in this scenario: for pure ideal gas regions or simple mixtures, analytical methods can be used efficiently; while for regions involving rigid gases or complex mixing interfaces, a stable iterative method is automatically activated. This intelligent scheduling ensures that the optimal strategy can be used for pressure solution throughout the entire complex flow field containing both gas and liquid phases, thus balancing computational efficiency and universality overall.

[0008] More importantly, this invention provides an indispensable and explicit physical basis for mandatory physical constraints. The negative infinity pressure parameter in the rigid gas equation of state assigns specific values ​​to the reference pressures of each phase, transforming the introduction of constraints based on the reference pressures of each phase from an abstract principle into an executable and physically meaningful operation: ensuring that the pressure of the rigid gas phase (water) is never lower than its cavitation threshold under any circumstances. It is this rigid constraint on the rigid gas phase that directly prevents non-physical negative pressure or cavitation errors in the liquid region caused by strong shock waves or expansion waves in examples such as underwater explosions, thereby cutting off the main cause of computational collapse. In summary, this scheme not only points out the typical application of this invention, but also highlights the synergistic strengthening effect of its two core mechanisms—adaptive solution and physical constraints—in solving the most challenging multiphase flow systems composed of ideal and rigid gases, ultimately ensuring the success and stability of simulations for such extreme conditions.

[0009] Furthermore, the flux calculation based on the HLLC format is performed in a three-dimensional global coordinate system. The intermediate wave velocity is used to synchronously calculate the non-conservative terms.

[0010] Furthermore, the direct and sole pressure input source for updating the thermodynamic state variables within the unit using the unit equilibrium pressure is the corrected unit equilibrium pressure after the constraint step of introducing reference pressures based on each phase.

[0011] In a second aspect, a pressure balance control device for compressible multiphase flow is provided, comprising: The flux calculation module is used to calculate the flux through the interface and the contribution of non-conservative terms to adjacent cells based on the flow field state on both sides of the mesh interface, so as to update the conserved quantities of each mesh cell. The pressure balance processing module is used to adaptively select a pressure solution path for each grid cell based on the state equation characteristics of each phase within it and the conservation quantities updated by the flux calculation module, to obtain an intermediate pressure value, and to introduce constraints based on the reference pressures of each phase onto the intermediate pressure value, so as to output a corrected cell balance pressure; and, The state update module is used to update the thermodynamic state variables within the unit using the corrected unit balance pressure output by the pressure balance processing module, in order to advance the flow field calculation.

[0012] Thirdly, a robust flow field update system for compressible multiphase flow simulation is provided, comprising sequentially connected components for unidirectional data transmission: The flux update unit is used for flux calculation based on the grid interface and outputs the updated value of the conserved quantity of the grid cell. The pressure processing and verification unit, connected to the flux update unit, is used to receive the updated value of the conserved quantity, perform pressure balance solution and forcibly apply physical constraints, and output the physically verified unit balance pressure. The state advancement unit, connected to the pressure processing and verification unit, is used to receive the physically verified unit balance pressure, update the flow field state, and advance the time step. The flow field data is transmitted from the flux update unit to the state propulsion unit after being processed by the pressure processing and verification unit.

[0013] Fourthly, a computing system for compressible multiphase flow simulation is provided, including a memory and a processor; The memory stores executable instructions; The processor is configured to be controlled by the executable instructions to perform the pressure balance control method for compressible multiphase flow as described in the first aspect.

[0014] Fifthly, the present invention also provides a computer-readable storage medium on which a computer program is stored, wherein the computer program, when executed by a processor, implements the steps of the method described in the first aspect.

[0015] The beneficial effects of this invention are as follows: The pressure balance control method for compressible multiphase flow of this invention fundamentally solves the global computational collapse problem caused by pressure solution failure in existing technologies by introducing a dual mechanism of adaptive solution path selection and mandatory physical constraints based on reference pressure. Specifically, this invention dynamically selects the optimal solution path based on the characteristics of the state equations of each phase, significantly improving the solution success rate for complex materials and extreme conditions while ensuring computational efficiency. Secondly, by applying constraints based on the reference pressure of each phase to the intermediate pressure values ​​obtained from numerical solutions, the calculation results are forcibly corrected to the physically feasible region, thereby completely eliminating the generation and spread of non-physical solutions such as negative pressure at the algorithm level. Through the synergistic effect of these two core measures, the unity of mathematical rationality and physical reality of the pressure balance solution is ensured, enabling the entire multiphase flow simulation system to possess unprecedented numerical robustness and computational stability when facing extreme conditions such as high pressure differentials and strong discontinuities.

[0016] Furthermore, this invention eliminates the local coordinate system rotation and inverse rotation operations required for each mesh surface in traditional algorithms by employing the HLLC format, which is executed directly in the three-dimensional global coordinate system. This significantly reduces the computational overhead and potential interpolation errors in large-scale three-dimensional computations, especially on unstructured meshes. Moreover, by utilizing an integrated process of synchronously calculating non-conservative terms at intermediate wave speeds, it ensures a high degree of coordination between flux, divergence, and non-conservative terms in both data and time, thereby improving the consistency and overall computational accuracy of the discrete format from the source.

[0017] Furthermore, this invention employs an adaptive mechanism for real-time judgment and path selection based on analytical root-finding of the state equations for each phase. This allows for the intelligent use of the optimal solver (fast analytical method or stable iterative method) in different physical regions within the same simulation (such as the pure gas phase region, the pure liquid phase region, and the mixed interface region). This dynamic allocation capability, combined with the Newton iteration method using the pressure value from the previous time step as the initial value, provides high-quality initial values ​​with clear physical meaning for iterative solutions. This significantly reduces the number of iterations and avoids divergence caused by inappropriate initial values, ensuring global stability while achieving efficient utilization of computational resources. Attached Figure Description

[0018] Figure 1 This is a schematic flowchart of a pressure balance control method for compressible multiphase flow according to an embodiment of the present invention. Figure 2 This is a schematic flowchart of a pressure balance control method for compressible multiphase flow according to another embodiment of the present invention; Figure 3This is a schematic flowchart of a pressure balance control method for compressible multiphase flow according to another embodiment of the present invention; Figure 4 This is a schematic flowchart of a pressure balance control method for compressible multiphase flow according to another embodiment of the present invention; Figure 5 This is a schematic flowchart of a pressure balance control method for compressible multiphase flow according to another embodiment of the present invention; Figure 6 This is a schematic flowchart of a pressure balance control method for compressible multiphase flow according to a specific embodiment of the present invention. Figure 7 This is a schematic flowchart of a pressure balance control method for compressible multiphase flow according to another specific embodiment of the present invention; Figure 8 This is a schematic flowchart of a pressure balance control method for compressible multiphase flow according to another specific embodiment of the present invention; Figure 9 This is a schematic diagram of the initial flow field for an underwater explosion simulation. Figure 10 A schematic diagram showing the flow field distribution changes in an underwater explosion simulation example; Figure 11 This is a schematic connection structure diagram of a pressure balance control device and a robust flow field update system for compressible multiphase flow according to an embodiment of the present invention. Figure 12 This is a topology diagram of a computer-readable storage medium disclosed in this invention. Detailed Implementation

[0019] The present invention will now be described in detail with reference to the embodiments shown in the accompanying drawings. However, it should be noted that these embodiments are not intended to limit the present invention. Equivalent changes or substitutions in function, method, or structure made by those skilled in the art based on these embodiments are all within the scope of protection of the present invention.

[0020] The technical solutions provided by the various embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0021] Example 1: like Figure 1 As shown in the figure, this embodiment provides a pressure balance control method for compressible multiphase flow, the method comprising: Step 102. Calculate the contribution of flux and non-conservation terms through the interface to adjacent cells based on the flow field states on both sides of the mesh interface, in order to update the conserved quantities of each mesh cell. The calculation of flux and non-conservation terms through the interface to adjacent cells is performed by calculating flux and simultaneously calculating non-conservation terms using the HLLC format. The flux calculation based on the HLLC format is performed in a three-dimensional global coordinate system.

[0022] In this step, all grid interfaces of the computational domain are traversed within each time step. The contributions of mass, momentum, energy flux, and inherent non-conservation terms to adjacent cells through the interface are calculated based on the reconstructed flow field state of the cells on both sides of the interface. This updates the stored conserved quantities (including phase mass, mixed momentum, and total energy) within all grid cells. Thus, on the one hand, the discrete framework of the finite volume method strictly guarantees the conservation of mass, momentum, and energy between cells, laying a mathematical foundation for physically correct flow field evolution. On the other hand, by synchronizing or integrating the calculation of non-conservation terms with flux calculation, two tasks can be completed simultaneously in a single interface traversal and data access. This eliminates the numerical asynchrony errors introduced by intermediate storage or step disconnection in traditional separate algorithms, significantly improving the locality of data processing and overall computational efficiency. This provides a more collaborative and consistent data source of conserved quantities for subsequent high-precision stress solutions.

[0023] like Figure 2 As shown, flux is calculated based on the HLLC scheme, and non-conservative terms are calculated simultaneously, including: Step 201. For each mesh interface, obtain the flow field state of its left and right elements and the interface normal vector. That is, before performing HLLC calculations for each mesh interface, first retrieve all necessary flow field state variables (such as density, velocity, pressure, phase fraction, etc.) of the left and right elements of that interface from storage, and simultaneously obtain the geometric properties of the interface (unit normal vector). This provides complete and accurate input data for the subsequent Riemann solver. Obtaining the normal vector is crucial in 3D calculations, as it defines the orientation of the interface and is the geometric basis for calculating normal flux. Therefore, efficiently and accurately completing data acquisition is a prerequisite for ensuring the correctness and robustness of the entire flux calculation, especially in complex unstructured meshes, where the implementation of this step directly relates to the algorithm's adaptability to complex geometries.

[0024] Step 202. Calculate the flux components and intermediate wave velocity of the interface based on the HLLC format, and calculate the contribution to the element velocity divergence based on the intermediate wave velocity.

[0025] In the three-dimensional global coordinate system, the obtained left and right states and normal vectors are directly substituted into the calculation formula of the HLLC (Harten-Lax-van Leer-Contact) approximate Riemann solver to obtain the normal flux components on the interface and a key intermediate variable—the intermediate wave velocity. S M Furthermore, this implementation utilizes the intermediate wave velocity generated during the HLLC solution process. S M Simultaneously calculate the velocity divergence of the interface with respect to its adjacent elements. u The HLLC format contributes significantly to flux calculation. Because it can accurately capture complex wave systems such as shock waves and contact discontinuities, it ensures high precision in flux calculation. Furthermore, calculations are performed directly in a three-dimensional global coordinate system, completely eliminating the time-consuming local coordinate system rotation and inverse rotation operations at each interface in traditional methods, greatly improving the efficiency of large-scale three-dimensional computation. Finally, this embodiment uses intermediate wave velocities not only to determine the flux state but also to instantly convert them into estimates of velocity divergence, providing timely and matching data for the next step of synchronously calculating non-conservative terms, serving as a key technological link in realizing the integrated process.

[0026] Step 203. Accumulate the flux components into the residuals of the corresponding cells, and simultaneously calculate the contribution of the non-conservative terms based on the contribution of the velocity divergence.

[0027] The calculated flux components, according to their directions (defined by the normal vector), are accumulated with opposite signs into the residual vectors (or right-hand side terms) of the corresponding units on the left and right sides of the interface, thus updating the time derivative of the conserved quantities. Simultaneously, based on the same intermediate wave velocity... S M The calculated velocity divergence contribution is used to simultaneously calculate the contribution of non-conservation terms caused by the interface (in the five-equation model, this mainly refers to the influence of the product of pressure work and velocity divergence on the phase fraction equation), and these contributions are accumulated in the phase fraction residual of the corresponding unit. With this setup, this embodiment of the invention achieves a true integrated calculation loop for flux, divergence, and non-conservation terms. Work that would traditionally require multiple independent loops or subsequent processing steps is compressed into a single, efficient interface traversal. This not only significantly reduces memory accesses and improves computational speed but also ensures that flux, divergence, and non-conservation terms are strictly synchronized in time, space, and logic. This fundamentally avoids numerical errors caused by data asynchrony, enhances the overall consistency and stability of the format, and prepares a highly consistent data environment for subsequent stress solutions facing extreme physical conditions.

[0028] Step 104. For each grid cell, based on the state equation characteristics of each phase within it and the updated conservation quantities, adaptively select the pressure solution path to obtain the intermediate pressure value, and introduce constraints based on the reference pressure of each phase on the intermediate pressure value to output the corrected cell equilibrium pressure.

[0029] This embodiment, for each grid cell that has completed flux updates, intelligently selects the most efficient or stable mathematical path to solve for a temporary intermediate pressure value based on the mathematical characteristics of the equations of state (EOS) of each phase fluid within it (especially its ability to analytically solve the internal energy equation of the mixture), and the updated conserved quantities (total energy, density of each phase, etc.). Subsequently, a mandatory safety rule based on physical laws is applied to this numerical solution, namely, verification and correction are performed using the inherent reference pressure of each phase material as a threshold, and finally, the corrected cell equilibrium pressure is output. Thus, through the dual guarantee of intelligent selection and physical enforcement, the dilemma of incompatibility between efficiency and stability faced by traditional pressure solutions under extreme conditions, as well as the fundamental risk of numerical solutions deviating from physical reality, are fundamentally solved. This ensures that no matter how complex the flow or how poor the initial guess, the pressure value ultimately passed to the updated flow field is always physically feasible, acting like a robust firewall to prevent the chain reaction of global computational collapse caused by the failure of a single-point pressure solution (such as the generation of negative pressure).

[0030] like Figure 3 As shown, based on the characteristics of the state equations of each phase within it, and the updated conserved quantities, the pressure solution path is adaptively selected, including: Step 301. Determine if the equations of state for each phase within the cell support analytical root-solving of the mixture's internal energy equation. Before solving for pressure in each cell, first query the equations of state objects used by all phases within it, and call a unified interface (such as supportsAnalyticalMixture()) to determine if these equations of state are capable of analytically root-solving the total internal energy equation of the current mixture, thus achieving dynamic optimal allocation of computational resources. The judgment logic in Step 301 does not use a single method, but rather dynamically allocates a suitable solver to each cell based on the mathematical characteristics of the material combinations within the cell. For example, for pure ideal gases or specific combinations, a fast analytical path can be used; for mixtures containing complex equations of state (such as Stiffened Gas), a stable iterative path is automatically adopted. In this way, computational efficiency and algorithm universality are improved at the global level.

[0031] Step 302. If yes, call the analytical solver submodule to directly calculate the intermediate pressure value using an explicit formula.

[0032] Step 303. If not, call the iterative solution submodule to perform numerical iterations starting from the preset initial pressure value until convergence to obtain the intermediate pressure value. The iterative solution submodule uses Newton's iteration method, and the preset initial pressure value is the pressure value of the previous time step for that unit.

[0033] It should be understood that when the judgment result in step 301 is yes, an analytical solution submodule is invoked to directly substitute the updated conserved quantities of the unit and the parameters of the state equations of each phase into a pre-derived explicit mathematical formula. The intermediate pressure value is obtained directly through finite algebraic operations, without any iteration, thus achieving extremely high computational efficiency while ensuring stability. For phase combinations that support analytical solutions (such as certain gas-gas mixtures), the overhead and convergence problems of iterative loops can be avoided. The calculation speed is extremely fast and highly deterministic, compressing calculations that might otherwise require dozens of iterations into a few formula evaluations. This is particularly suitable for regions with a large proportion of the flow field and favorable material properties, making a key contribution to reducing overall computational costs. If the judgment result in step 301 is no, an iterative solution submodule is invoked. This module uses a preset initial pressure value (usually set to the successful pressure value of the previous time step of the unit). p old Starting with the pressure at the previous time step, the algorithm iteratively solves the nonlinear equations, which are the coupling of the internal energy equation of the mixture and the state equations of each phase, using numerical methods such as Newton's iteration method, until the solution meets the convergence tolerance, thus obtaining the intermediate pressure value. This provides a universal and stable solution capability. Because Newton's iteration method has a second-order convergence speed and high efficiency, and because using the pressure at the previous time step as the initial value is very close to the true solution, it can significantly reduce the number of iterations and greatly avoid iteration divergence caused by poor initial values. Therefore, it can robustly handle any complex material combination and extreme state, which is an important pillar of the algorithm's robustness.

[0034] The intermediate pressure values ​​obtained from the solution are constrained based on the reference pressures of each phase, which is achieved through comparison and weighted averaging. For example... Figure 4 As shown, it specifically includes: Step 401. Compare the intermediate pressure value with the reference pressure of each phase, and determine the local pressure for synthesis of each phase based on the comparison results. The reference pressure is usually the negative infinity pressure in the equation of state parameter of each phase, which represents the lowest physically permissible pressure threshold of that phase material (such as the cavitation pressure of water).

[0035] For the intermediate pressure value obtained in steps 302 or 303, each phase within the unit is iterated and compared with a reference pressure defined in the phase's equation of state, which has a clear physical meaning. Based on the comparison result, a local pressure is determined for each phase: if the intermediate pressure value is higher than the phase's reference pressure, the local pressure is set to the intermediate value; if it is lower than or equal to the phase's reference pressure, the local pressure is forcibly set to the reference pressure, thus achieving a mandatory mapping from the mathematical solution to the physical solution. Step 401 is no longer a purely mathematical calculation but introduces an inviolable physical law: the pressure of any phase cannot be lower than its physical lower limit. Therefore, step 401 actively detects and corrects non-physical negative pressures that may result from numerical errors.

[0036] Step 402. Calculate the corrected unit equilibrium pressure by weighted averaging of the local pressures based on the volume fraction of each phase. After determining the physically consistent local pressure for each phase, a weighted average is calculated according to the volume fraction of each phase in the current unit. p final =Σ( α k p local,k This process yields the modified unit equilibrium pressure used to drive the flow field update, thereby reconstructing the overall pressure that conforms to the multiphase mixture equilibrium assumption after applying physical constraints. This ensures that the final output pressure strictly satisfies the physical constraints of each phase (thus achieving global stability) while maintaining consistency with the mixture model.

[0037] Step 106. Update the thermodynamic state variables within the unit using the unit equilibrium pressure to advance the flow field calculation. That is, using the unit equilibrium pressure, update the thermodynamic state variables of all phases within the unit in reverse order through independent equations of state for each phase, including but not limited to the density, temperature, internal energy, and sound velocity of each phase. Subsequently, based on these updated and consistent thermodynamic states, recalculate and finally confirm the updated values ​​of conserved variables (such as total energy) for this time step, thereby completing the solution for all flow field variables in the current time step and reliably advancing the entire flow field state to the next time step. This ensures that the pressure data flowing into the time-propagation engine has undergone both adaptive solving and physical constraints, making it physically absolutely safe. Therefore, even under extreme physical conditions, the entire system can perform iterative cycles based on a reliable thermodynamic foundation, thus achieving a transition from single-step stability to full-process stability.

[0038] like Figure 5 As shown, before performing calculations based on the flow field states on both sides of the mesh interface, the following steps are also included: Step 101. Initialize and clear the residual vectors of all mesh cells to zero.

[0039] It should be understood that at the beginning of each time step or at the beginning of a new round of nonlinear iteration, clearing all values ​​of the vector or array storing all grid cell residuals (i.e., the time derivative terms composed of flux and non-conservative contributions) to provide a clean and consistent starting point for the calculation of flux and contribution in step 102 is a key prerequisite for achieving a highly robust and deterministic framework. This ensures that the residual accumulation calculation at each time step is not contaminated by residual values ​​from the previous time step or iteration, preventing the accidental accumulation and propagation of errors, and laying a reliable data foundation for all subsequent accurate calculations based on residual contributions (such as flux accumulation and divergence calculation).

[0040] Therefore, the pressure balance control method for compressible multiphase flow in this embodiment fundamentally solves the global computational collapse problem caused by pressure solution failure in existing technologies by introducing a dual mechanism of adaptive solution path selection and mandatory physical constraints based on reference pressure. It dynamically selects the optimal solution path based on the characteristics of the state equations of each phase, significantly improving the success rate of solutions for complex materials and extreme conditions while ensuring computational efficiency. By applying constraints based on the reference pressure of each phase to the intermediate pressure values ​​obtained from the numerical solution, the calculation results are forcibly corrected to the physically feasible region, thereby completely eliminating the generation and spread of non-physical solutions such as negative pressure at the algorithm level. Thus, the synergistic effect of these two core measures ensures the unity of mathematical rationality and physical reality in the pressure balance solution, enabling the entire multiphase flow simulation system to possess unprecedented numerical robustness and computational stability when facing extreme conditions such as high pressure differentials and strong discontinuities. This solves the technical problem of existing numerical simulation methods for compressible multiphase flow causing global computational collapse due to the generation and spread of non-physical solutions caused by the lack of physical constraints in the pressure balance solution process under extreme conditions.

[0041] The method in this embodiment is based on the five-equation model framework of compressible multiphase flow. This model, by assuming instantaneous pressure and velocity equilibrium between each phase, solves the following set of governing equations at the mixture level to describe the flow involving shock wave and interface interactions: 1. Mass conservation equations for each phase

[0042] The two equations describe the independent conservation of mass for each phase. Among them, Representing the k The volume fraction of the phase, for a two-phase flow, satisfies ; The density of each phase (the first phase) k Phase density); This represents the velocity field (common velocity vector) of the mixture.

[0043] 2. Momentum Conservation Equation for Mixtures

[0044] This equation describes the overall momentum conservation of the mixture, where, For the density of the mixture, For pressure, It is a unit tensor.

[0045] 3. Energy Conservation Equation for Mixtures

[0046] This equation ensures the conservation of the total energy of the mixture, where, Total energy per unit mass, including internal energy per unit mass And kinetic energy.

[0047] 4. Phase fractional transport equations

[0048] By introducing an interface compression item Accurately describe the evolution dynamics of phase interfaces in compressible flows. K This is the interface compression factor.

[0049] In the actual numerical solution of the above model, pressure p It is not a direct variable, but a crucial parameter that must be inversely derived from thermodynamic relations through the equation of state. Specifically, after updating the conserved quantities at each step, the pressure must be determined by solving the following internal energy relationship of the mixture derived from the energy equation:

[0050] , It is a function defined by the independent state equations of each phase (such as the ideal gas equation or the rigid gas equation). Traditional single-solution strategies are prone to computational failure or non-physical solutions (such as negative pressure) at this stage, leading to global collapse. To solve the problem of global computational collapse caused by pressure solution failure under extreme conditions, in a specific embodiment of the present invention, the pressure balance control method for compressible multiphase flow solves the above-mentioned technical problem by combining adaptive pressure solution with mandatory physical constraints. Figure 6 The implementation process is explained as follows: I. Initialization and Residual Clearing At the beginning of each physical time step Δ t Before the calculation, the vectors of the residuals (i.e. time derivative terms) of the conservation equations stored in all grid cells are cleared to provide a clean starting point for the subsequent accumulation of interface contributions.

[0051] II. Simultaneous Calculation of Flux and Non-Conservation Terms in the Five-Equation Multiphase Flow System under Global Coordinates The convection term is calculated using the HLLC approximate Riemann solver, and the interface flux is determined based on the wave structure model. This step updates the conserved quantities for each grid cell. U The flux is calculated directly in the three-dimensional global coordinate system, and the contribution of non-conservative terms is completed simultaneously, avoiding the overhead of traditional coordinate rotation. The conserved state vector U is shown in Equation 1: (1) Combination Figure 7 To illustrate, for each internal mesh interface in the computational domain, the following steps are performed: 1. Obtain interface information: Read the left-side unit of the interface. L and the right unit R Reconstructed flow field state U L , U R and the unit normal vector of the interface. n= ( n x ,n y ,n z ) and interface area A.

[0052] 2. Calculate the intermediate wave velocity and pressure: (1) Calculate the interface normal velocity (inverter velocity):

[0053] (2) Estimate the freezing sound velocity on the left and right sides. .

[0054] (3) Calculate the left and right wave velocity estimates:

[0055] (4) Calculate the intermediate wave velocity s M and pressure in the middle region p :

[0056] 3. Determine flux and accumulate contribution: Select the flux formula based on the wave velocity relationship, and accumulate the flux and divergence contribution in one cycle.

[0057] (1) Calculate the contribution to the element velocity divergence: The contribution of this interface to the velocity divergence of the left and right elements is - s M A and s M A. This contribution will be temporarily stored for the calculation of non-conservative terms.

[0058] (2) Calculate and sum the HLLC flux Based on wave speed, select one of the following four cases ( F L , F R , , One of the fluxes F ( U The flux is calculated using the formula. The branching is used as shown in Equation 2: (2) in, ,in, and These represent the inverter speeds of the left and right unit interfaces, respectively. ,in, For unit normal vector, , and These represent the velocities in the three directions of the three-dimensional flow field. For intermediate wave speeds, and These are the left and right wave velocities, respectively. ,in , It is the freezing speed of sound of the left and right units. If Then the flux is completely determined by the left-hand state, through the formula F L Calculate flux; if Then the flux is completely determined by the right-hand state, through the formula F R Calculate flux; if The flux is then calculated using the left satellite region value, through... Flux calculation, the left satellite region state calculation is shown in Equation 3; if The flux is then calculated using the right satellite region value, through... Flux is calculated, and the right satellite region state is calculated as shown in Equation 4.

[0059] The state calculation for the left star region is shown in Equation 3: (3) The right satellite region state is calculated as shown in Equation 4: (4) Flux accumulation: Calculate the interface flux according to the above formula. F Then, multiply it by the interface area A, and add it to the residual of adjacent elements according to the finite volume method rules: RL -= ·A, R R += • A (Mass, Momentum, and Energy Equations).

[0060] 4. Simultaneously calculate and accumulate non-conservative terms. (1) Calculate the velocity divergence contribution: using the calculated intermediate wave velocity The interface contributes to the velocity divergence of the left and right elements respectively. and ( V (This is for the unit volume). This contribution is directly used in the calculation of non-conservative terms.

[0061] (2) Calculate and sum the non-conservation terms: The non-conservation terms in the phase fractional equation are ,in H ( U As shown in Equation 5: (5) Using the velocity divergence contribution generated by this interface, its contribution to the phase fraction residual is calculated in real time and accumulated to the corresponding cell.

[0062] 5. After traversing all interfaces, each cell accumulates the flux contribution and non-conservation term contribution from all adjacent interfaces, thus obtaining the complete residual. R Used to update conserved quantities: U n+1 = U n +(Δ t ) R (in R (This already includes the effect of volume averaging).

[0063] III. Adaptive Hybrid Pressure Solution and Physical Constraint Correction Combination Figure 8 To explain, for each unit, based on the updated conserved quantities... U n+1 Solve for the physically feasible equilibrium pressure.

[0064] 1. Adaptive path selection: Query all related paths within the unit. k The state equations are used to call a unified interface (such as supportsAnalyticalMixture()) to determine whether all analytical root finding of the internal energy equations of mixtures (including synchronous calculation of non-conservation terms) is fully supported.

[0065] 2. Analytical Method Branch: If analytical methods are supported, the analytical solution submodule is invoked. For example, for a mixture of an ideal gas and a rigid gas, the solution can be obtained by solving for... p The intermediate pressure value is obtained directly from the quadratic equation. p calc :

[0066] In the formula, α k Volume fraction of each phase, ρ k For density, γ k For specific heat ratio, p ∞ , k For pressure constant, b k For co-volume, For reference internal energy and internal energy of mixtures e mix .

[0067] 3. Iterative method branch: If analytical methods are not supported, the Newton iteration submodule is called.

[0068] (1) Set the initial value for iteration: p 0 =p n (Previous time step pressure).

[0069] (2) Perform Newton iteration: For the th m The next iteration, based on current pressure valuation. p m The internal energy of each phase is calculated from the equations of state for each phase. e k ( p m ,ρ k and speed of sound c k .

[0070] (3) Calculate the residual function of internal energy of the mixture .

[0071] (4) Calculate the residual derivative .

[0072] (5) Update pressure: .

[0073] (6) Repeat until | F ( p m )|< ε Set the final pressure to bep calc .

[0074] 4. Physical constraint correction: Read the minimum pressure threshold (reference pressure) defined by the equation of state for each phase, such as the minimum pressure for a rigid gas phase. p ∞,k .

[0075] (1) Initialization: phase index Ultimate pressure p final =0.

[0076] (2) Traverse all phases k : Compare p calc and- p ∞,k (Reference pressure): If p calc >- p ∞,k The local pressure of that phase p local,k = p calc ;like p calc ≤ - p ∞,k Then force correction p local,k =- p ∞,k .

[0077] (3) Weighted composite final pressure: p final =Σ( α k p local,k ), p final This refers to the unit equilibrium pressure, which has been corrected for physical constraints to ensure that the pressure of each phase is not lower than its physical lower limit.

[0078] IV. State Variable Updates and Time Step Progression 1. Update the thermodynamic state: using a determined... p final Solve the equations of state for each phase to update the internal energy of each phase. e k ,temperature T k and the actual speed of sound under current pressure. c k All thermodynamic state variables.

[0079] 2. Synchronous Update of Conserved Quantities: To ensure thermodynamic consistency, based on the updated internal energy, the total energy in the conserved quantities is recalculated and synchronously corrected. ρE .

[0080] 3. Advance the time step: Set the current time... t Increase Δ t The updated flow field state U n+1 Use the initial value for the next time step, return to the first step, and repeat until the simulation ends.

[0081] In summary, the robust compressible multiphase flow pressure balance control method of this embodiment dynamically matches different materials and physical states by adaptively selecting the pressure solution path, balancing computational efficiency and universality. It also forcibly introduces physical constraints based on phase state reference pressure to mathematically ensure the physical feasibility of the pressure solution, fundamentally preventing the generation and spread of non-physical solutions such as negative pressure. Simultaneously, by utilizing HLLC flux calculation directly executed in a three-dimensional global coordinate system, the overhead of traditional coordinate rotation is eliminated, improving computational efficiency. Finally, through a compact process of integrated calculation of flux, divergence, and non-conservative terms, a high degree of data synergy is achieved across all stages, thereby constructing a complete numerical simulation scheme with a robust chain from flux updates and pressure solutions to state propagation, capable of stably coping with extreme high pressure differential conditions such as underwater explosions.

[0082] Example: Underwater Explosion Simulation To verify the robustness and effectiveness of the aforementioned method in extreme high pressure differentials and strongly discontinuous multiphase flow problems, a typical underwater explosion simulation example is provided below. This example simulates the explosive expansion of high-pressure air bubbles in still water, involving huge initial pressure ratios (approximately 1e9 Pa / 1e5 Pa = 10000:1) and density ratios, which poses a severe test to the stability of numerical methods and their ability to capture strong shock waves.

[0083] 1. Case setup and initial conditions The computational domain is a rectangular region 4.0 meters long and 3.0 meters high, discretized using a 400×300 unstructured quadrilateral mesh. A high-pressure air bubble with an initial radius of 0.12 meters, centered at (2.0 meters, 1.2 meters), is surrounded by still water. The calculation employs explicit time-margining with a fixed time step of 8.0e-7 seconds. All boundaries are set as reflection-free supersonic exit boundary conditions to simulate an open water environment. The initial flow field configuration is as follows. Figure 9 As shown (the high-pressure air bubble-phase 1 in the center of the computational domain is surrounded by water-phase 2, and different initial phase fractions are set in the upper and lower halves of the domain with y=1.5 meters as the boundary to define the bubble position), the specific parameters are detailed in Tables 1 and 2.

[0084] Table 1 Initial flow field parameters for underwater explosion example

[0085] Table 2 Fluid state parameters for underwater explosion case study

[0086] 2. Calculation Results and Analysis Figure 10 The flow field evolution results at different times (t=0.00032s, 0.00096s, 0.00256s) calculated using the method of this invention are shown, where the left column (i.e. Figure 10 (a) Figure 10 (c) Figure 10 In the diagram (e), the phase fraction (Alpha1) distribution is represented, and the right column (i.e.) Figure 10 (b) Figure 10 (d) in Figure 10 In the diagram, (f) represents the pressure (P) distribution. (From...) Figure 10 The flow field simulation results shown at different times demonstrate that the simulation using the method of this embodiment can stably calculate up to 3.0 milliseconds of physical time without crashing. In contrast, if the physical constraint correction function is disabled, the program will diverge and crash at approximately 0.2 milliseconds due to non-physical solutions such as negative pressure. This directly proves the crucial role of the physical constraint module in maintaining the robustness of global computation.

[0087] The method in this embodiment successfully simulated the core physical process of an underwater explosion: high-pressure bubble expansion → shock wave propagation → bubble contraction. The phase fraction throughout the entire field remained strictly within the range of [0,1], and both density and pressure were positive, conforming to physical laws. The bubble dynamics and pressure wave propagation characteristics met theoretical expectations, verifying the method's excellent ability to capture strong shock waves, large-deformation gas-liquid interfaces, and maintain long-term computational stability. Therefore, the pressure balance control method in this embodiment, through its core mechanism of adaptive pressure solution and mandatory physical constraints, can effectively cope with the numerical challenges brought about by extreme conditions such as underwater explosions, fundamentally preventing computational collapse and achieving highly robust and high-fidelity numerical simulation of compressible multiphase flow.

[0088] Example 2: by Figure 11This embodiment provides a pressure balance control device 112 for compressible multiphase flow, comprising: a flux calculation module 1120, used to calculate the flux through the interface and the contribution of non-conservative terms to adjacent cells based on the flow field state on both sides of the grid interface, so as to update the conserved quantities of each grid cell; a pressure balance processing module 1121, used for each grid cell, based on the state equation characteristics of each phase and the conserved quantities updated by the flux calculation module, to adaptively select the pressure solution path to obtain an intermediate pressure value and introduce constraints based on the reference pressure of each phase on the intermediate pressure value, so as to output the corrected cell equilibrium pressure; and a state update module 1122, used to update the thermodynamic state variables within the cell using the corrected cell equilibrium pressure output by the pressure balance processing module, so as to advance the flow field calculation.

[0089] The pressure balance control device 112 provided in this embodiment uses modular design to hardware-encode and solidify the complex numerical process, ensuring that the calculation method is executed rigorously and reliably, avoiding logical deviations that may occur in software implementation, thereby fundamentally improving the overall stability and determinism of the simulation system. Furthermore, the adaptive selection function of the pressure balance processing module 1121 can intelligently call the optimal solution path (fast analytical method or stable iterative method) based on the real-time state equation characteristics of each phase material within the unit, achieving dynamic optimal allocation of computational resources for complex flow regions, greatly improving computational efficiency while ensuring versatility. The pressure balance processing module 1121 forcibly introduces a function based on the reference pressure constraints of each phase, forming a built-in physical safety valve that can actively correct the intermediate pressure values ​​obtained from mathematical solutions to a physically feasible range, preventing the generation and propagation of non-physical solutions such as negative pressure. This is crucial for solving the global computational collapse caused by the failure of a single-point pressure solution. In addition, the state update module 1122 must rely on the physically verified balance pressure output by the preceding module to advance the flow field, forming a strict data-driven-constraint verification closed loop. The coordinated operation of the modules in the pressure balance control device 112 enables the device to continuously output physically reasonable and stable results in engineering simulations involving extreme pressure and density ratios, such as underwater explosions and high-speed impacts, achieving high robustness and high reliability at the level of compressible multiphase flow numerical simulation devices.

[0090] It should be noted that the scheme or principle involved in the pressure balance control device 112 for compressible multiphase flow in this embodiment is the same as that described in Embodiment 1, and the same or similar contents will not be described in detail.

[0091] Example 3: Continue with Figure 11This embodiment provides a robust flow field update system 113 for compressible multiphase flow simulation, comprising: a flux update unit 1130, used to calculate and output the conservation update values ​​of grid cells based on the flux of the grid interface; a pressure processing and verification unit 1131, connected to the flux update unit, used to receive the conservation update values, perform pressure balance solutions and forcibly apply physical constraints, and output the physically verified cell equilibrium pressure; and a state advancement unit 1132, connected to the pressure processing and verification unit, used to receive the physically verified cell equilibrium pressure, update the flow field state and advance the time step; wherein, the flow field data is transferred from the flux update unit 1130 to the state advancement unit 1132 after being processed by the pressure processing and verification unit 1131.

[0092] The robust flow field update system 113 provided in this embodiment achieves fundamental system-level hardening by constructing a mandatory serial data processing architecture of flux update → pressure processing and verification → state advancement. This configuration elevates the key mechanisms ensuring computational physics reliability (pressure balance solution and physical constraints) from an optional algorithmic step to an indispensable and mandatory verification step in the entire data flow. This architecture mandates that all flow field data, after completing the flux update, must and can only flow through the pressure processing and verification unit 1131, which has physical constraint functions, before entering the final state advancement unit 1132. This completely blocks, at the system level, the path of potential non-physical solutions (such as negative pressure) that have not undergone physical verification from directly contaminating and damaging the subsequent flow field. This eliminates the possibility of global computational collapse caused by local pressure solution failure at its source, giving the entire simulation system unprecedented inherent stability and resilience against extreme conditions.

[0093] It should be noted that the robust flow field update system 113 for compressible multiphase flow simulation in this embodiment involves the same scheme or principle as that described in Embodiment 1, and the same or similar contents will not be described in detail.

[0094] Example 4: This invention also provides a computing system 111, which may include a processor, a memory, and a computer program stored in the memory and executable on the processor, the computer program implementing the above when executed by the processor. Figure 1-8 The various processes shown in the embodiment of the pressure balance control method for compressible multiphase flow achieve the same technical effect, and will not be described again here to avoid repetition.

[0095] Combination Figure 12As shown, this embodiment also discloses a specific implementation of a computer-readable storage medium 1200. This computer-readable storage medium 1200 can be configured wholly or partially in a physical computer, server, cluster server, or data center.

[0096] In this embodiment, the computer-readable storage medium 1200 stores computer program instructions 1201, which are read and executed by a processor 1202 to perform the steps in the pressure balance control method for compressible multiphase flow as disclosed in Embodiment 1.

[0097] Optionally, the computer-readable storage medium 1200 can be configured as a server, and the server runs on a physical device used to build a private cloud, hybrid cloud, or public cloud. The computer-readable storage medium 1200 can also be configured as random access memory (RAM), read-only memory (ROM), programmable read-only memory (PROM), erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), etc.

[0098] The computer-readable storage medium 1200 is used to store a program, and the processor 1202, upon receiving an execution instruction, executes the pressure balance control method for compressible multiphase flow disclosed in Embodiment 1.

[0099] Meanwhile, the processor 1202 disclosed in this embodiment may be an integrated circuit chip with signal processing capabilities. The processor 1202 can be a general-purpose processor, including a central processing unit (CPU), a network processor (NP), etc.; it can also be a digital signal processor (DSP), an application-specific integrated circuit (ASIC), a field-programmable gate array (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components. It can implement or execute the various methods, steps, and logic block diagrams disclosed in the embodiments of this invention. The general-purpose processor can be a microprocessor or any conventional processor.

[0100] The technical solution of the computer-readable storage medium 1200 disclosed in this embodiment that is the same as that in Embodiment 1 and / or Embodiment 2 is described in Embodiment 1 and / or Embodiment 2, and will not be repeated here.

[0101] The detailed descriptions listed above are merely specific descriptions of feasible embodiments of the present invention, and are not intended to limit the scope of protection of the present invention. All equivalent embodiments or modifications made without departing from the spirit of the present invention should be included within the scope of protection of the present invention.

[0102] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the invention can be implemented in other specific forms without departing from its spirit or essential characteristics. Therefore, the embodiments should be considered in all respects as exemplary and non-limiting, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of equivalents of the claims are intended to be included within the present invention. No reference numerals in the claims should be construed as limiting the scope of the claims.

[0103] Furthermore, it should be understood that although this specification describes embodiments, not every embodiment contains only one independent technical solution. This narrative style is merely for clarity. Those skilled in the art should consider the specification as a whole, and the technical solutions in each embodiment can also be appropriately combined to form other embodiments that can be understood by those skilled in the art.

Claims

1. A pressure balance control method for compressible multiphase flow, characterized by, The method includes: The flux and non-conservative terms at the interface are calculated based on the flow field state on both sides of the mesh interface to contribute to adjacent cells, so as to update the conserved quantities of each mesh cell. For each grid cell, the pressure solution path is adaptively selected based on the state equation characteristics of each phase within it and the updated conservation quantities to obtain the intermediate pressure value. Constraints based on the reference pressure of each phase are introduced on the intermediate pressure value to output the corrected cell equilibrium pressure. The thermodynamic state variables within the unit are updated using the unit's equilibrium pressure to advance the flow field calculation.

2. The method according to claim 1, characterized in that, The calculation of the flux and non-conservation terms through the interface contributes to adjacent cells, and the flux is calculated based on the HLLC format while the non-conservation terms are calculated simultaneously.

3. The method according to claim 2, characterized in that, The calculation of flux based on the HLLC format and the simultaneous calculation of non-conservative terms include: For each mesh interface, obtain the flow field state of its left and right cells and the interface normal vector; The flux components and intermediate wave velocity of the interface are calculated based on the HLLC format, and the contribution to the unit velocity divergence is calculated based on the intermediate wave velocity. The flux components are accumulated into the residuals of the corresponding units, and the contribution of the non-conservative terms is calculated synchronously based on the contribution of the velocity divergence.

4. The method according to claim 3, characterized in that, The flux calculation based on the HLLC format is performed in a three-dimensional global coordinate system.

5. The method according to claim 1, characterized in that, The adaptive selection of the pressure solution path based on the state equation characteristics of each phase and the updated conserved quantities includes: Determine whether the state equations of each phase within the unit support analytical root-finding of the internal energy equation of the mixture; If so, the analytical solver submodule is invoked to directly calculate the intermediate pressure value using an explicit formula; If not, the iterative solution submodule is invoked to perform numerical iterations starting from a preset initial pressure value until convergence to obtain an intermediate pressure value.

6. The method according to claim 1 or 5, characterized in that, The constraint on the intermediate pressure value obtained by the solution, based on the reference pressure of each phase, is achieved by comparison and weighted averaging.

7. The method according to claim 6, characterized in that, The method of comparison and weighted averaging includes: The intermediate pressure value is compared with the reference pressure of each phase, and the local pressure for synthesis of each phase is determined based on the comparison results. The corrected unit equilibrium pressure is obtained by weighted averaging of the local pressure based on the volume fraction of each phase.

8. The method according to claim 5, characterized in that, The iterative solution submodule uses the Newton-Raphson method, and the preset initial pressure value is the pressure value of the unit at the previous time step.

9. A pressure balance control device for compressible multiphase flow, characterized in that, include: The flux calculation module is used to calculate the flux through the interface and the contribution of non-conservative terms to adjacent cells based on the flow field state on both sides of the mesh interface, so as to update the conserved quantities of each mesh cell. The pressure balance processing module is used to adaptively select the pressure solution path for each grid cell based on the state equation characteristics of each phase within it and the conservation quantity updated by the flux calculation module, to obtain an intermediate pressure value, and to introduce constraints based on the reference pressure of each phase on the intermediate pressure value, so as to output the corrected cell balance pressure. as well as, The state update module is used to update the thermodynamic state variables within the unit using the corrected unit balance pressure output by the pressure balance processing module, in order to advance the flow field calculation.

10. A robust flow field update system for compressible multiphase flow simulation, characterized in that, include: The flux update unit is used to calculate the conservation update value of the output grid cell based on the flux of the grid interface; The pressure processing and verification unit, connected to the flux update unit, is used to receive the updated value of the conserved quantity, perform pressure balance solution and forcibly apply physical constraints, and output the physically verified unit balance pressure. The state advancement unit, connected to the pressure processing and verification unit, is used to receive the physically verified unit balance pressure, update the flow field state, and advance the time step. The flow field data is transmitted from the flux update unit to the state propulsion unit after being processed by the pressure processing and verification unit.

Images