A two-phase flow simulation method based on density-based implicit solution

Abstract

The application relates to a two-phase flow simulation method based on density-based implicit solving. The method comprises the following steps: adopting a double-time stepping method to discretize a physical time into a pseudo-time steady-state problem, and linearizing to construct a Jacobian matrix. In spatial discretization, a high-order format is adopted to calculate a numerical flux and a first-order upwind format is adopted to calculate the Jacobian matrix through a decoupling strategy for the same interface, so that a full-coupling sparse linear system is formed. A parallel algebraic multigrid solver is used for efficient solving and variable updating. After convergence, a chemical relaxation model is introduced for phase change correction to ensure thermodynamic equilibrium. Finally, the loop is advanced to a preset physical time, and the result is output. The method can effectively support high-fidelity numerical simulation of a compressible gas-liquid two-phase flow process in a liquid engine thrust chamber.

Description

Technical Field

[0001] This application relates to the field of engine numerical simulation technology, and in particular to a two-phase flow simulation method based on density-based implicit solution. Background Technology

[0002] The injection, atomization, mixing, and combustion process of propellant in the thrust chamber of an underwater liquid engine is essentially an interaction process between gas-liquid two-phase flow and shock waves under highly compressible conditions. During this process, the high-speed injected propellant forms a large number of micron-sized droplets and bubbles. The strong shock waves generated under the high-pressure environment of the combustion chamber interact with these bubbles / droplets, inducing droplet breakup, interface deformation, and the generation and evolution of vortex structures. This process directly determines the propellant atomization efficiency and the uniformity of fuel-oxidizer mixing, thus affecting the engine's combustion stability, specific impulse performance, and operational reliability.

[0003] High-fidelity numerical simulation of shock wave-bubble / droplet interaction processes is a core technical means for performance optimization, scheme iteration, and fault diagnosis of underwater liquid engines. This simulation places stringent demands on numerical methods: on the one hand, it needs to accurately capture complex flow characteristics such as strong shock wave propagation, high-density phase interface evolution, and baroclinic vortex generation to ensure the accuracy of simulation results; on the other hand, it needs to adapt to the unstructured mesh of the complex configuration of the engine thrust chamber to meet the large-scale, high-resolution simulation requirements of full-size 3D models at the engineering level. Currently, for simulations of compressible two-phase flow problems, density-based solution methods are the mainstream approach. Among them, the explicit time integration scheme is widely used for simulation verification of benchmark problems such as shock wave-bubble interactions due to its simple implementation logic and strong numerical robustness, and it is also the mainstream scheme for two-phase flow engineering simulation in aerospace propulsion systems.

[0004] However, existing explicit time integration schemes suffer from insurmountable technical bottlenecks, failing to meet the core requirements of underwater liquid engine simulation. Strictly limited by CFL stability conditions, large-scale simulation computations are extremely inefficient, hindering rapid iteration of engineering solutions. In benchmark simulations of shock-bubble interactions, to ensure numerical stability, the CFL number in explicit schemes is typically set to only 0.2–0.3, severely restricting the time step. In contrast, engineering simulations of underwater liquid engines, aiming to accurately capture micrometer-level droplet breakup and interface evolution, require high-resolution unstructured meshes with tens or even hundreds of millions of elements. This reduces the time step of explicit schemes to nanoseconds, requiring millions of iterations to simulate a single millisecond-level engine operation. This results in computation cycles lasting weeks or even months, leading to extremely high computational costs that cannot support rapid iteration and multi-condition optimization of advanced engine solutions. Furthermore, in high-density two-phase flow scenarios, an inherent contradiction exists between numerical stability and simulation accuracy. The density ratio of gas-liquid two-phase flows (such as liquid oxygen-fuel gas, kerosene-high-temperature fuel gas) in underwater liquid engines can reach thousands of times, far exceeding the gas-to-gas density ratio in conventional shock-bubble benchmark tests. When shock waves propagate across the phase interface, the strong acoustic impedance contrast induces complex shock wave refraction, reflection, and baroclinic vortex generation. Existing explicit schemes, to avoid numerical divergence, inevitably introduce excessive numerical dissipation, directly leading to blurred phase interface capture, distorted shock wave position, and excessively rapid vortex structure evolution dissipation. This makes it impossible to accurately predict propellant atomization and mixing characteristics, thus affecting the accuracy of engine combustion performance design. Existing improved implicit schemes cannot meet the simulation requirements of underwater liquid engine scenarios. To overcome the limitations of CFL (Continuous Fluid Dynamics) conditions, the industry has conducted research on implicit solution methods for two-phase flows. However, existing implicit solvers generally suffer from core defects: First, under high-order spatial discretization, the analytical derivation of the Jacobian matrix is ​​cumbersome and error-prone. Using automatic differentiation techniques leads to a significant increase in computational and memory overhead for matrix assembly, offsetting the time step advantage of implicit solutions. Second, for large-scale linear systems with unstructured meshes, there is a lack of efficient solution strategies, resulting in slow iterative convergence. Especially in rigid two-phase flow scenarios with strong shock waves and high density ratios, robustness is insufficient, and iterative divergence is prone to occur. Third, most existing solutions only complete two-dimensional benchmark case verification and cannot meet the large-scale parallel simulation requirements of full-size three-dimensional models of underwater liquid engines. Summary of the Invention

[0005] Therefore, it is necessary to provide a density-based implicit solution-based two-phase flow simulation method that can effectively support high-fidelity numerical simulation of compressible gas-liquid two-phase flow processes in the thrust chamber of underwater liquid engines, addressing the aforementioned technical problems.

[0006] A two-phase flow simulation method based on density-based implicit solution, the method comprising:

[0007] Step S1: Obtain the unstructured mesh data of the physical domain to be simulated and initialize the conservation variable vector containing temperature, density, velocity, pressure and total energy information in each control volume element; Step S2: Based on the dual-time-step scheme, construct the physical time discrete terms with the unit average values ​​of the conserved variable vector at different physical time levels, reconstruct the nonlinear implicit system after the original physical time discretization into a steady-state problem in the pseudo-time domain, and linearize the residual terms in the steady-state problem to construct a linearized algebraic equation system containing the Jacobian matrix. Step S3: Using a decoupling strategy, numerical flux terms and Jacobian matrix terms are constructed at the spatial discretization level respectively. Specifically, for the same control volume interface, a higher-order spatial discretization scheme is used to calculate the numerical flux terms, and a first-order upwind scheme is used to calculate the Jacobian matrix terms. Step S4: Substitute the numerical flux term into the residual term of the linearized algebraic equation system, and use the Jacobian matrix term as the coefficient matrix of the linearized algebraic equation system to form a fully coupled global sparse linear system. Step S5: Solve the global sparse linear system using the parallel algebraic multigrid method and update the conservation variable vector at the current pseudo-time step; Step S6: Determine whether the pseudo-time iteration has converged according to the preset convergence condition. If it has converged, proceed to step S7; otherwise, return to step S3 to continue the iteration. Step S7: Apply the chemical relaxation model to the converged conserved variable vector to perform phase transition correction, so that the two phases reach thermodynamic equilibrium. Step S8: Determine whether the physical time has advanced to the preset endpoint. If yes, output the simulation result; otherwise, return to step S2 to proceed to the next physical time step.

[0008] The aforementioned density-based implicit solution method for two-phase flow simulation first reconstructs the nonlinear implicit system after physical time discretization using a dual-time-step scheme, and improves the accuracy of Jacobian matrix calculation by combining automatic differentiation techniques. This effectively overcomes the strict limitations of the CFL condition in explicit methods, allowing for larger physical time steps. A decoupling strategy is employed between flux calculation and Jacobian matrix construction, retaining high-order space discretization for accurate flux calculation while using a simplified upwind scheme for Jacobian matrix calculation, reducing the complexity of derivative calculation while maintaining solution accuracy. Furthermore, the GMRES solver, integrating algebraic multigrids, improves the solution efficiency for large-scale sparse linear systems. This application significantly improves computational efficiency while maintaining comparable accuracy to explicit schemes, effectively supporting high-fidelity numerical simulation of compressible gas-liquid two-phase flow processes in the thrust chamber of underwater liquid engines, and meeting the needs of rapid iteration in engineering solutions. Attached Figure Description

[0009] Figure 1This is a flowchart illustrating a two-phase flow simulation method based on density-based implicit solution in one embodiment. Figure 2 This is a schematic diagram illustrating the strategy for decoupling flux assessment and Jacobian matrix construction in one embodiment. Figure 3 This is a schematic diagram of an scalable implicit solver framework in one embodiment; Figure 4 This is a schematic diagram of the initial configuration for simulating the shock wave-bubble interaction in another embodiment; Figure 5 This is a schematic diagram of shock wave-bubble interaction in a helium bubble / air case in one embodiment; Figure 5 (a) is an experimental schlieren image under helium bubble / air conditions. Figure 5 (b) A plot of numerical results obtained at selected times using the explicit Euler scheme in the helium bubble / air case. Figure 5 (c) A plot of numerical results obtained at selected times using the single-step implicit BDF2 scheme in the helium bubble / air case; Figure 6 This is a shock wave-bubble interaction diagram for a krypton bubble / air case in one embodiment; Figure 6 (a) is an experimental schlieren image under krypton bubble / air conditions. Figure 6 (b) A plot of numerical results obtained at selected times using the explicit Euler scheme in the krypton bubble / air case. Figure 6 (c) A plot of numerical results obtained at selected times using the single-step implicit BDF2 scheme in the krypton bubble / air case; Figure 7 In one embodiment = 274 Figure showing the verification results of the two-time-step method for the helium bubble / air condition at point s; Figure 8 This is an initial configuration diagram of the three-dimensional shock-bubble interaction in one embodiment; Figure 9 This is a three-dimensional shock wave-bubble interaction diagram in one embodiment; Figure 9 (a) is a three-dimensional shock-bubble interaction diagram in the yz-plane view. Figure 9 (b) is a three-dimensional shock wave-bubble interaction diagram in the xy-plane view; Figure 10 This is a representative snapshot of the three-dimensional shock-bubble interaction at different times in one embodiment. Detailed Implementation

[0010] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0011] In one embodiment, such as Figure 1 As shown, a two-phase flow simulation method based on density-based implicit solution is provided, including the following steps: Step S1: Obtain the unstructured mesh data of the physical domain to be simulated and initialize the conservation variable vector containing temperature, density, velocity, pressure and total energy information in each control volume element.

[0012] Neglecting viscosity and thermal conductivity effects, the governing equations are expressed in conservation form as follows:

[0013] In the formula , For temperature and density, the conserved variable vector is defined as follows:

[0014] In the formula, For the velocity field, For mixed pressure, For unit tensors, This represents the total energy.

[0015] Step S2: Based on the dual-time-step scheme, construct the physical time discrete terms using the unit average values ​​of the conserved variable vector at different physical time levels. Reconstruct the nonlinear implicit system after the original physical time discretization into a steady-state problem in the pseudo-time domain. Linearize the residual terms in the steady-state problem to construct a linearized algebraic equation system containing the Jacobian matrix.

[0016] The dual-time-step method introduces a pseudo-time derivative term, constructing a physical-time discrete term using the unit average value of the conserved variable vector at different physical time levels. This transforms the nonlinear implicit problem after physical-time discretization into a steady-state problem in the pseudo-time domain. This method effectively handles the nonlinearity and stiffness of the two-phase flow control equations, enhances numerical stability, and allows for larger physical time steps.

[0017] Specifically, considering the homogeneous part of the governing equation in step S1, the governing equation can be written as follows:

[0018] Within the framework of the Finite Volume Method (FVM), the semi-discrete finite volume formula can be written as follows:

[0019] in, For the common surface between units, From The unit normal vector, It represents a set of indices for adjacent units.

[0020] Will - Applying the backward difference formula (BDF) time integral to the semi-discrete finite volume formula, the fully discrete finite volume scheme is expressed as follows:

[0021] In the above formula, by selecting parameter pairs ( , ) =(1,0 ) Recover the first-order backward difference formula (BDF1), and through ( , )=(1,0.5 ) The second-order backward difference formula (BDF2) is obtained. Parameters Time processing of control flux terms: = 1 corresponds to a fully explicit format. = 0 produces a fully implicit format, with an intermediate value of 0. < <1 leads to a mixed explicit-implicit (IMEX) formula.

[0022] The fully discrete finite volume scheme can be solved efficiently using the pseudo-time step (or two-time step) method. The governing equations can be rewritten as:

[0023] in, This represents the BDF time format parameter. This represents the BDF time format parameter. Represents the pseudo-time iterative solution vector. Indicates pseudo-time. Indicates the first The unit average value of the physical time step. Indicates the first The unit average value of the physical time step. Indicates the physical time step. This indicates the degree of implicitness of flux. Represents the residual.

[0024] As pseudo-time →∞ Solution vector Converging to steady-state solution Advancing using implicit integration schemes Pseudo-temporal evolution:

[0025] Here, superscript r is The pseudo-time iteration index is defined as follows: .

[0026] Moreover, for term R( r+1 ) conducted a discussion on r First-order Taylor expansion linearization:

[0027] Rearranging the above equation, the resulting linear system can be written as:

[0028] in, This represents the BDF time format parameter. This represents the BDF time format parameter. Represents the pseudo-time iterative solution vector. Indicates pseudo-time. Indicates the first The unit average value of the physical time step. Indicates the first The unit average value of the physical time step. Indicates the physical time step. This indicates the degree of implicitness of flux. Indicates pseudo time step, Represents the identity matrix. This represents the solution vector at the current pseudo-time step. Indicates the contribution of the flux term. Represents the contribution of the time discrete term. Represents the Jacobian matrix. This represents the solution vector for the next pseudo-time step. This represents the residual at the current pseudo-time step. Represents the residual of the reference state.

[0029] Therefore, the two-time-step method can be interpreted as a relaxation-based iterative process for solving fully implicit systems. The introduction of the pseudo-time derivative transforms the original nonlinear implicit problem into a series of linearized systems that reach a steady state in the pseudo-time domain. In the limiting case, when λ = 0, the pseudo-time step tends to infinity, and the two-time-step formula degenerates into a fully implicit discretized direct Newton iteration. Furthermore, if only one pseudo-time iteration is performed in each physical time step, the method simplifies to a linearized implicit scheme, where the solution is updated based on a single linearization. Thus, the two-time-step method provides a unified framework connecting relaxation-based iterative solvers and fully implicit Newton-type solvers, offering enhanced robustness and flexibility for rigid flow problems.

[0030] The effectiveness and robustness of the aforementioned implicit solution strategy largely depend on the accurate and stable computation of the numerical flux and its corresponding Jacobian matrix. In particular, designing well-coordinated flux evaluation and Jacobian matrix construction strategies is crucial for ensuring numerical stability, convergence efficiency, and accuracy under conditions of strong nonlinearity and discontinuity.

[0031] Step S3: Using a decoupling strategy, numerical flux terms and Jacobian matrix terms are constructed at the spatial discretization level respectively. Specifically, for the same control volume interface, a higher-order spatial discretization scheme is used to calculate the numerical flux terms, and a first-order upwind scheme is used to calculate the Jacobian matrix terms.

[0032] The decoupling strategy separates flux calculation from Jacobian matrix construction: flux evaluation can employ high-order discretization schemes, such as TVD and MUSCL, to ensure solution accuracy; while Jacobian derivative calculation uses a simplified formula based on the upwind direction, significantly reducing computational complexity without sacrificing accuracy. Figure 2 As shown, the numerical flux on the cell surface is calculated using a high-order scheme, while the Jacobian matrix depends only on the conserved variables of the adjacent control volumes, forming a compact template and a sparse matrix structure.

[0033] Specifically, to balance accuracy and efficiency, this application employs a decoupled flux evaluation and Jacobian matrix construction strategy, such as... Figure 2 As shown. Decoupling occurs at both the spatial discretization and algebraic solver levels. The vector R( on the right-hand side of the linear system formula...) Various spatial discretization schemes (including high-resolution formulas) can be used for evaluation. Depending on the local flow characteristics and accuracy requirements, first-order upwind, QUICK, higher-order TVD schemes, or higher-order reconstruction methods, such as the MUSCL-Hancock scheme, can be employed. This flexibility can accurately address discontinuous and smooth flow characteristics even under high CFL number conditions.

[0034] Conversely, the Jacobian matrix in the formula for linear systems It uses a first-order upwind discrete construction. For example... Figure 2 As shown in the left figure, the numerical flux on the cell surface depends only on the conserved variables of the two adjacent control volumes. Therefore, the resulting Jacobian matrix contains non-zero entries that are only related to neighboring cells, resulting in a compact template and a sparse matrix structure.

[0035] This decoupling formula makes the numerical accuracy of flux calculation independent of the complexity of the Jacobian matrix. Although higher-order schemes introduce extended templates and strongly nonlinear terms, these effects are deliberately excluded from the construction of the Jacobian matrix. Therefore, the linear system maintains good conditions and computational efficiency, improving the overall spatial accuracy of the solution.

[0036] Step S4: Substitute the numerical flux term into the residual term of the linearized algebraic equation system, and use the Jacobian matrix term as the coefficient matrix of the linearized algebraic equation system to form a fully coupled global sparse linear system.

[0037] Step S5: Solve the global sparse linear system using the parallel algebraic multigrid method and update the conservation variable vector at the current pseudo-time step.

[0038] like Figure 3 As shown, this method integrates the YHAMG algebraic multigrid library and uses the AMG-preprocessed GMRES solver to solve large-scale sparse linear systems. YHAMG supports MPI / OpenMP hybrid parallelism and is suitable for large-scale parallel architectures.

[0039] Step S6: Determine whether the pseudo-time iteration has converged according to the preset convergence condition. If it has converged, proceed to step S7; otherwise, return to step S3 to continue the iteration.

[0040] The convergence condition for pseudo-time iteration is usually set to the residual decreasing to a preset threshold, or the relative change of the solution vector being less than a set tolerance.

[0041] Step S7: Apply a chemical relaxation model to the converged conserved variable vector to perform phase transition correction, so that the two phases reach thermodynamic equilibrium.

[0042] The chemical relaxation model achieves thermodynamic equilibrium between the two phases through the Gibbs free energy equilibrium condition. The density and total internal energy of the mixture remain unchanged before and after relaxation, and the pressure and temperature at equilibrium are obtained by solving algebraic equations.

[0043] The calculation process of the chemical relaxation model specifically includes: Solve After the homogeneous part, the phase transition problem is solved by achieving thermodynamic equilibrium between the two phases under conditions of Gibbs free energy non-equilibrium; this process is called relaxation. Essentially, it involves solving a partial differential equation: (1) However, if we consider a relaxation coefficient ≥ ∞, meaning that each relaxation process reaches equilibrium... 1= 2. Then, the problem can be simplified to solving an algebraic equation using conservation relations. In this case, the density of the mixture remains unchanged before and after relaxation. (2) The total internal energy of the mixture remains constant. (3) Therefore, we can obtain information about ( , The quadratic equation of )) (4) The coefficients can be calculated as follows: (5) (6) (7) At the same time, pressure and temperature satisfy the saturation relationship: (8) Among them, coefficient , , , for: (9) By combining the quadratic equation with the saturation curve, the Newton-Raphson iterative method can be used to solve the algebraic equations of the above formulas (1)-(9).

[0044] Step S8: Determine whether the physical time has advanced to the preset endpoint. If yes, output the simulation result; otherwise, return to step S2 to proceed to the next physical time step.

[0045] If the current physical time is less than the preset total simulation time, update the physical time step index and return to step S2 to continue the calculation of the next physical time step; otherwise, output the simulation results.

[0046] The aforementioned density-based implicit solution method for two-phase flow simulation first reconstructs the nonlinear implicit system after physical time discretization using a dual-time-step scheme, and improves the accuracy of Jacobian matrix calculation by combining automatic differentiation techniques. This effectively overcomes the strict limitations of the CFL condition in explicit methods, allowing for larger physical time steps. A decoupling strategy between flux calculation and Jacobian matrix construction is employed, retaining high-order space discretization for accurate flux calculation while using a simplified upwind scheme for Jacobian matrix calculation, reducing the complexity of derivative calculation while maintaining solution accuracy. Furthermore, the GMRES solver, integrating algebraic multigrid preprocessing, improves the solution efficiency for large-scale sparse linear systems. This application significantly improves computational efficiency while maintaining comparable accuracy to explicit schemes, effectively supporting high-fidelity numerical simulation of compressible gas-liquid two-phase flow processes in the thrust chamber of underwater liquid engines, and meeting the needs of rapid iteration in engineering solutions.

[0047] In one embodiment, based on a two-time-step scheme, the nonlinear implicit system after discretization of the original physical time is reconstructed into a steady-state problem in the pseudo-time domain, including: Based on the two-time-step scheme, the nonlinear implicit system after discretization of the original physical time is reconstructed into the governing equations of a steady-state problem in the pseudo-time domain.

[0048] in, This represents the BDF time format parameter. This represents the BDF time format parameter. Represents the pseudo-time iterative solution vector. Indicates pseudo-time. Indicates the first The unit average value of the physical time step. Indicates the first The unit average value of the physical time step. Indicates the physical time step. This indicates the degree of implicitness of flux. Represents the residual.

[0049] Specifically, the governing equations reconstruct the nonlinear implicit system discretized from the original physical time into a steady-state problem in the pseudo-time domain. When the pseudo-time iteration converges, and The equation degenerates into the original physical-time discrete scheme. This dual-time-step method allows for time steps much larger than those limited by the explicit CFL method, effectively improving computational efficiency, while ensuring the solution accuracy of nonlinear systems through pseudo-time iteration.

[0050] In one embodiment, a system of linearized algebraic equations containing the Jacobian matrix is ​​constructed as follows:

[0051] in, This represents the BDF time format parameter. This represents the BDF time format parameter. Represents the pseudo-time iterative solution vector. Indicates pseudo-time. Indicates the first The unit average value of the physical time step. Indicates the first The unit average value of the physical time step. Indicates the physical time step. This indicates the degree of implicitness of flux. Indicates pseudo time step, Represents the identity matrix. This represents the solution vector at the current pseudo-time step. Indicates the contribution of the flux term. Represents the contribution of the time discrete term. Represents the Jacobian matrix. This represents the solution vector for the next pseudo-time step. This represents the residual at the current pseudo-time step. Represents the residual of the reference state.

[0052] Specifically, the linearized algebraic equations are obtained by performing a first-order Taylor expansion on the residual terms. The coefficient matrix on the left side contains the pseudo-time derivative contribution and the Jacobian matrix contribution, while the vector on the right side contains the flux term and the time-discrete term. The solution to this system of equations is obtained through pseudo-time iteration, gradually approximating the solution of the physical time-discrete system. Automatic differentiation techniques are used to calculate the Jacobian matrix, ensuring that the derivative accuracy reaches machine precision and avoiding the tediousness and errors of manual derivation.

[0053] In one embodiment, the numerical flux term is computed using a high-order spatial discretization scheme, including: The numerical flux term is calculated using an HLLC-type approximate Riemann solver; the HLLC-type approximate Riemann solver is...

[0054] in, , For left and right wave velocities, For contact discontinuity wave velocity, , For inviscid flux in the left and right states, , These are the conserved variables for the intermediate state.

[0055] In one embodiment, the numerical flux term is calculated using an HLLC-type approximate Riemann solver as follows:

[0056] in, The phase density represents the intermediate state. Representing different stages. Indicates left or right. Density of a mixture in its intermediate state Indicates the contact discontinuity wave velocity. The pressure representing the intermediate state. This represents the total energy density of the intermediate state.

[0057] In one embodiment, the contact discontinuity wave velocity is

[0058] in, Indicates pressure on the right side. Indicates pressure on the left side. This indicates the density of the mixture on the left. Indicates the velocity on the left. For left and right wave velocities, For right wave velocity, This indicates the density of the mixture on the right. Indicates the velocity on the right; The expressions for the phase density, total energy density, and pressure of the intermediate state are as follows:

[0059] in, Represents phase density, Indicates left and right wave velocities. Represents the total energy in the initial state. Indicates the density of the mixture on the left and right. The pressure is on the left and right sides.

[0060] Specifically, this HLLC-type approximate Riemannian solver can accurately capture complex wave systems such as shock waves, contact discontinuities, and rarefied waves. In two-phase flow simulations, due to the large density ratio and strong acoustic impedance contrast on both sides of the phase interface, the HLLC scheme can effectively suppress numerical dissipation and maintain the clarity of the phase interface and the accuracy of the shock wave location. In the above formula, wave velocity... , and It is set to a constant value during the construction of the Jacobian matrix to simplify the calculation of the derivative.

[0061] In one embodiment, when constructing the Jacobian matrix terms in the first-order upwind scheme, the Jacobian matrix of each control volume element contains only the off-diagonal elements corresponding to the control volume elements directly adjacent to the control volume element.

[0062] Specifically, when constructing the Jacobian matrix using a first-order upwind scheme, the numerical flux on the cell surface depends only on the conserved variables of two adjacent control volumes. Therefore, each row of the Jacobian matrix contains only off-diagonal elements related to the cells directly adjacent to the target cell. This compact template structure ensures the diagonal dominance of the linear system, improving the convergence speed and numerical stability of the iterative solution. Compared to the extended template introduced by higher-order discretization schemes, this strategy significantly reduces matrix storage overhead and solution complexity.

[0063] In one embodiment, the global sparse linear system is stored and assembled using a block sparse row matrix format.

[0064] Specifically, the block sparse row matrix scheme combines variables with the same sparsity pattern into blocks, optimizing memory access patterns using the underlying coupling structure. For block-structured linear systems generated by the four-equation two-phase flow model, the BSR scheme can maintain good memory locality and parallel scalability, reduce indirect addressing overhead in matrix-vector multiplication operations, and improve solution efficiency.

[0065] In one embodiment, the parallel algebraic multigrid method includes using an algebraic multigrid as a preprocessor and combining it with the generalized minimum residual method to solve the global sparse linear system.

[0066] Specifically, such as Figure 3 As shown, this method integrates the YHAMG algebraic multigrid library, which is designed for large-scale parallel architectures and supports MPI / OpenMP hybrid parallelism. This is existing technology and will not be elaborated upon further in this application. The GMRES solver, preprocessed with AMG, effectively handles rigid linear systems in two-phase flow simulations caused by large CFL numbers, strong nonlinearity, and severe density contrast, ensuring robustness and efficiency in the convergence of linear systems.

[0067] In one embodiment, the chemical relaxation model achieves phase transition correction by solving algebraic equations under Gibbs free energy equilibrium conditions, where the density of the mixture and the total internal energy of the mixture remain unchanged before and after relaxation, and solves a set of quadratic equations about pressure and temperature through Newton-Raphson iteration.

[0068] In a specific embodiment, the numerical solver was validated based on experimental data from shock-bubble interactions conducted by Layers and Tayer, which provides a comprehensive description of the experimental setup. Specifically, helium and krypton were chosen as the bubble-filling gases, representing the light and heavy gas configurations, respectively. The ambient medium was air at atmospheric pressure and temperature. For each gas configuration, a plane shock wave with a Mach number of 1.5 was used to trigger the shock-bubble interaction. Time reference. 0 = 0 This corresponds to the moment when the shock wave interacts with the upstream interface of the bubble.

[0069] The thermodynamic closure used in the numerical framework corresponds to the ideal gas law. For all cases considered in this application, the initial conditions for the numerical simulation are as follows: Figure 4 As shown. To systematically evaluate the accuracy, robustness, and computational efficiency of the proposed numerical framework under different flow conditions, three different numerical simulations were performed: (i) helium bubble / air case: evaluating the accuracy and efficiency of the implicit solver under the light gas configuration; (ii) krypton bubble / air case: evaluating the accuracy and efficiency of the implicit solver under the heavy gas configuration; (iii) validating the dual-time-stepping method with relaxed CFL constraints, highlighting its ability to significantly improve computational efficiency while maintaining accuracy. The computational domain used is approximately 2.4 × 10⁻⁶. 6 Discretization is performed using a hexahedral mesh composed of individual elements, ensuring sufficient spatial resolution to accurately capture the shock dynamics and subsequent bubble deformation in all test cases. In this configuration, the bubble is filled with helium, which is significantly lighter than the surrounding air. Due to the strong contrast in density and acoustic impedance across the interface, baroclinic vorticity is generated during shock wave propagation and refraction. Consequently, the shock wave propagating into the helium bubble propagates faster than the incident shock wave in the air, leading to complex interface deformation and vortex dynamics. The current numerical simulations aim to reproduce these experimentally observed features and evaluate the accuracy and efficiency of the proposed implicit solver. As shown in Figure 5, Figure 5 (a) is an experimental schlieren image under helium bubble / air conditions. Figure 5 (b) A plot of numerical results obtained at selected times using the explicit Euler scheme in the helium bubble / air case. Figure 5 (c) A plot of numerical results obtained at selected times using the single-step implicit BDF2 scheme in the helium bubble / air case. Figure 5Experimental schlieren images at several representative moments were compared with numerical results obtained using the explicit Euler scheme and the single-step implicit BDF2 scheme. First-order spatial discretization was employed in both numerical methods to ensure consistency of the comparison. Experiments showed that following the compression phase was the development of vorticity along the bubble interface, which evolved over a longer timescale. The numerical results successfully captured the formation process of the penetrating jet along the axis of symmetry of the flow direction. This high-speed jet penetrated the helium bubble, producing two consecutive flow structures. The vorticity concentrated in the downstream vortex ring. The second structure was a ring with lower vorticity. Subsequently, the two rings separated, indicating that they had different velocities. Simultaneously, the vortex ring increased in size, while the upstream vortex ring decreased in size. These key flow characteristics are qualitatively consistent with experimental observations. Quantitatively, both the explicit and implicit solvers demonstrated considerable accuracy in predicting global shock propagation, interface deformation, and vortex evolution, indicating that the implicit formula does not compromise the fidelity of the solution even at large time steps. The explicit Euler solver was run at a CFL number of 0.298 to ensure numerical stability in this simulation, while the single-step implicit BDF2 solver remained stable and accurate at a larger CFL number of 5.953.

[0070] From a computational efficiency perspective, the implicit solver, as described in this application, offers significant advantages. For this helium bubble / air case, the total computation time required by the explicit solver is approximately 6.22 times that of the implicit solver. This substantial reduction in computational cost, achieved without sacrificing accuracy, clearly demonstrates the effectiveness of the implicit time integration strategy for impact-bubble interaction problems involving strong material property contrasts.

[0071] In the krypton bubble / air case, the bubble is filled with krypton, which is significantly heavier than the surrounding air. Similar to the case of helium, the strong contrast in density and acoustic impedance across the material interface leads to baroclinism during the shock-bubble interaction. However, due to the higher density and acoustic impedance of krypton, the overall flow evolution exhibits distinct characteristics compared to the lighter gas configuration. In the early stages of the interaction, the flow dynamics are dominated by compression. The shock wave propagating into the krypton bubble travels slower than the incident shock wave in the air, consistent with experimental schlieren observations. Eddy generation is present in this initial stage, but remains secondary. Over time, the eddy effect becomes increasingly significant and controls the subsequent deformation of the bubble interface.

[0072] like Figure 6 As shown, Figure 6 (a) is an experimental schlieren image under krypton bubble / air conditions. Figure 6 (b) A plot of numerical results obtained at selected times using the explicit Euler scheme in the krypton bubble / air case. Figure 6(c) A plot of numerical results obtained at selected times using the single-step implicit BDF2 scheme in the krypton bubble / air case. Figure 6 Experimental schlieren images at several representative moments were compared with numerical results obtained using explicit Euler schemes and single-step implicit BDF2 schemes. The numerical simulations successfully reproduced the main flow characteristics observed experimentally. In particular, along the flow direction, filaments first appear near the top and bottom regions of the bubble. These filaments then coil and merge to form vortex rings, the size and intensity of which increase as the interactions develop. The overall bubble deformation, vortex ring formation, and shock wave structure predicted by the numerical simulations qualitatively agree very well with the experimental observations.

[0073] Both numerical methods demonstrate comparable accuracy in capturing shock propagation, interface evolution, and vortex dynamics. In the current simulation, both solvers employ first-order spatial discretization to ensure a fair comparison. The explicit Euler solver maintains numerical stability at a CFL number of 0.2021, while the single-step implicit BDF2 solver remains stable and accurate at a CFL number of 4.1454. From a computational efficiency perspective, the implicit solver again shows a significant advantage. For the krypton / air case, the explicit solver requires approximately 6.23 times the total computation time of the implicit solution. These results further confirm that the proposed implicit time integration strategy can achieve considerable efficiency improvements while maintaining solution accuracy, even for impact-bubble interactions involving heavy gas configurations and strong material property contrasts.

[0074] In a specific embodiment, to further relax the time step restrictions imposed by the stability constraints, this embodiment evaluates the dual time step formula. Validation is performed using a helium bubble / air configuration. = 274 The flow field at physical time s is compared. This instant corresponds to a stage where penetrating flow and vortex structures are well-developed, making it suitable for evaluating accuracy and robustness.

[0075] Figure 7 Numerical results are presented as the number of flow fluxes (CFLs) increases. The reference solution, calculated using a single-step implicit BDF2 scheme, yields a CFL number of 5.953. Using a two-time-step formula, the CFL number can be further increased to 14.882 while maintaining the main flow characteristics. Interface deformation, jet structure, and vortex morphology remain essentially unchanged, indicating that the solution accuracy is well maintained even at this much larger time step.

[0076] In the dual-time-step simulation, the relaxation parameter is set to λ = 0.5, and when the residual decreases to 10... 7The pseudo-time iteration terminates at this point. Although the introduction of pseudo-time iteration increases the number of linear system solutions per physical time step, resulting in a slight decrease in efficiency compared to the single-step BDF2 scheme, the overall computational performance remains good.

[0077] For the same physics simulation time, the two-timestep method still achieves a speedup of approximately 2.37 compared to the explicit Euler solver. These results demonstrate that the two-timestep method is an effective approach that can further expand the allowed number of CFLs while maintaining acceptable accuracy and providing considerable efficiency gains on explicit time integration.

[0078] In this embodiment, the proposed implicit solver is evaluated through a large-scale 3D simulation of spherical impact-bubble interaction. This case represents the first extension of the implicit formula to a full 3D configuration, aiming to evaluate its robustness and computational performance on large-scale meshes. The initial problem settings are as follows: Figure 8 As shown, a spherical helium bubble is embedded in still air and impacted by a plane shock wave with a Mach number of 1.2 propagating from right to left. Due to the significant density and acoustic impedance contrast between helium and air, this interaction generates complex three-dimensional flow characteristics during the subsequent deformation and collapse of the bubble.

[0079] like Figure 9 As shown, Figure 9 (a) is a three-dimensional shock-bubble interaction diagram in the yz-plane view. Figure 9 (b) shows the 3D shock-bubble interaction diagram in the xy-plane view. The computational domain is discretized using an unstructured tetrahedral mesh with a control volume of approximately 10,718,760, which is necessary to fully address the 3D shock structure, interface deformation, and vortex motion generated during the interaction. The simulation was performed using the implicit time integration framework developed in this work, demonstrating its applicability to large-scale 3D problems. A single-step implicit BDF2 scheme was used, with first-order spatial discretization and a CFL number of 4. Representative snapshots of the shock-bubble interaction at selected time instants are shown below. Figure 10 As shown.

[0080] In the early stages, when the incident shock wave impacts the interface, the helium bubble is subjected to intense compression. Due to its high density and large acoustic impedance contrast, the transmitted shock wave propagates rapidly within the bubble, causing significant interface deformation. As the interaction continues, the bubble collapses and gradually twists into a characteristic crescent-shaped structure, accompanied by the formation of a shear layer and the development of vortex motion near the interface.

[0081] In the later stages, the flow becomes increasingly three-dimensional, with complex structures related to bubble collapse, jetting phenomena, and vortex dynamics. The smooth evolution of shock waves and phase interfaces, coupled with the absence of spurious oscillations or numerical instabilities, clearly demonstrates the robustness of the proposed implicit solver in large-scale three-dimensional simulations.

[0082] It should be understood that, although Figure 1 The steps in the flowchart are shown sequentially as indicated by the arrows, but these steps are not necessarily executed in the order indicated by the arrows. Unless otherwise specified herein, there is no strict order in which these steps are executed, and they can be performed in other orders. Figure 1 At least some of the steps in the process may include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these sub-steps or stages is not necessarily sequential, but can be executed in turn or alternately with other steps or at least some of the sub-steps or stages of other steps.

[0083] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0084] The embodiments described above are merely illustrative of several implementation methods of this application, and while the descriptions are specific and detailed, they should not be construed as limiting the scope of this application. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of this application, and these modifications and improvements all fall within the protection scope of this application. Therefore, the protection scope of this application should be determined by the appended claims.

Claims

1. A two-phase flow simulation method based on density-based implicit solution, characterized in that, Numerical simulation of compressible gas-liquid two-phase flow processes in the thrust chamber of an underwater liquid engine, the method comprising: Step S1: Obtain the unstructured mesh data of the physical domain to be simulated and initialize the conservation variable vector containing temperature, density, velocity, pressure and total energy information in each control volume element; Step S2: Based on the dual-time-step format, construct physical time discrete terms using the unit average values ​​of the conserved variable vector at different physical time levels, reconstruct the nonlinear implicit system after the original physical time discretization into a steady-state problem in the pseudo-time domain, and linearize the residual terms in the steady-state problem to construct a linearized algebraic equation system containing the Jacobian matrix. Step S3: Using a decoupling strategy, the numerical flux term and the Jacobian matrix term are constructed at the spatial discretization level respectively. Specifically, for the same control volume interface, the numerical flux term is calculated using a high-order spatial discretization scheme, and the Jacobian matrix term is calculated using a first-order upwind scheme. Step S4: Substitute the numerical flux term into the residual term of the linearized algebraic equation system, and use the Jacobian matrix term as the coefficient matrix of the linearized algebraic equation system to form a fully coupled global sparse linear system. Step S5: Solve the global sparse linear system using the parallel algebraic multigrid method and update the conservation variable vector at the current pseudo-time step; Step S6: Determine whether the pseudo-time iteration has converged according to the preset convergence condition. If it has converged, proceed to step S7; otherwise, return to step S3 to continue the iteration. Step S7: Apply the chemical relaxation model to the converged conserved variable vector to perform phase transition correction, so that the two phases reach thermodynamic equilibrium. Step S8: Determine whether the physical time has advanced to the preset endpoint. If yes, output the simulation result; otherwise, return to step S2 to proceed to the next physical time step.

2. The method according to claim 1, characterized in that, Based on the two-time-step scheme, the nonlinear implicit system after discretization of the original physical time is reconstructed into a steady-state problem in the pseudo-time domain, including: Based on the two-time-step scheme, the nonlinear implicit system after discretization of the original physical time is reconstructed into the governing equations of a steady-state problem in the pseudo-time domain. in, This represents the BDF time format parameter. This represents the BDF time format parameter. Represents the pseudo-time iterative solution vector. Indicates pseudo-time. Indicates the first The unit average value of the physical time step. Indicates the first The unit average value of the physical time step. Indicates the physical time step. This indicates the degree of implicitness of flux. Represents the residual.

3. The method according to claim 1, characterized in that, Construct a system of linearized algebraic equations containing the Jacobian matrix as follows: in, This represents the BDF time format parameter. This represents the BDF time format parameter. Represents the pseudo-time iterative solution vector. Indicates pseudo-time. Indicates the first The unit average value of the physical time step. Indicates the first The unit average value of the physical time step. Indicates the physical time step. This indicates the degree of implicitness of flux. Indicates the pseudo time step. Represents the identity matrix. This represents the solution vector at the current pseudo-time step. Indicates the contribution of the flux term. Represents the contribution of the time discrete term. Represents the Jacobian matrix. This represents the solution vector for the next pseudo-time step. This represents the residual at the current pseudo-time step. Represents the residual of the reference state.

4. The method according to claim 1, characterized in that, The numerical flux term is calculated using a high-order spatial discretization scheme, including: The numerical flux term is calculated using an HLLC-type approximate Riemann solver; the HLLC-type approximate Riemann solver is... in, , For left and right wave velocities, For contact discontinuity wave velocity, , For inviscid flux in the left and right states, , These are the conserved variables for the intermediate state.

5. The method according to claim 4, characterized in that, The numerical flux term is calculated using an HLLC-type approximate Riemann solver as follows: in, The phase density represents the intermediate state. Representing different stages. Indicates left or right. Density of a mixture in its intermediate state Indicates the contact discontinuity wave velocity. The pressure representing the intermediate state. This represents the total energy in the intermediate state.

6. The method according to claim 5, characterized in that, The method further includes: The contact discontinuity wave velocity is in, Indicates pressure on the right side. Indicates pressure on the left side. This indicates the density of the mixture on the left. Indicates the velocity on the left. For left and right wave velocities, For right wave velocity, This indicates the density of the mixture on the right. Indicates the velocity on the right; The expressions for the phase density, total energy density, and pressure of the intermediate state are as follows: in, Represents phase density, Indicates left and right wave velocities. Represents the total energy in the initial state. Indicates the density of the mixture on the left and right. The pressure is on the left and right sides.

7. The method according to claim 1, characterized in that, When constructing the Jacobian matrix terms using the first-order upwind scheme, the Jacobian matrix of each control volume unit contains only the off-diagonal elements corresponding to the control volume units directly adjacent to the control volume unit.

8. The method according to claim 1, characterized in that, The global sparse linear system is stored and assembled using a block sparse row matrix format.

9. The method according to claim 1, characterized in that, The parallel algebraic multigrid method includes using an algebraic multigrid as a preprocessor and combining it with the generalized minimum residual method to solve the global sparse linear system.

10. The method according to claim 1, characterized in that, The chemical relaxation model achieves phase transition correction by solving algebraic equations under Gibbs free energy equilibrium conditions, where the density and total internal energy of the mixture remain unchanged before and after relaxation, and the system of quadratic equations about pressure and temperature is solved by Newton-Raphson iteration.

Images