Universal multi-moment finite volume method for unstructured grid compressible multiphase flow

Abstract

The invention relates to a universal multi-moment finite volume method for compressible multiphase flow of an unstructured grid. According to the method, two types of independent variables, namely a volume integral average value and a point value, are respectively defined on a grid unit and a vertex of the grid unit; the core innovation of the method lies in an updating strategy of a point value: firstly, constructing a vertex reconstruction polynomial based on a reconstruction template; then, on a flux calculation template formed by geometric midpoints of the template, directly adopting a Roe, an HLL or an AUSM and other general Riemann solver to calculate flux; and finally, flux divergence is obtained through least square fitting so as to complete updating. According to the design, the dependence of a traditional method on a specific Riemann solver based on a derivative is broken through, and the compatibility with various mainstream solvers is realized, so that a high-order multi-moment framework is successfully expanded to compressible multiphase flow simulation. According to the method, flux calculation of the volume integral average value and the point value is unified, the algorithm is simple, excellent geometric adaptability is achieved, and high-order precision and high robustness are shown in single-phase and two-phase flow simulation.

Description

Technical Field

[0001] This invention relates to the field of computational fluid dynamics, and in particular to a general multi-moment finite volume method for compressible multiphase flow in unstructured meshes. Background Technology

[0002] Compressible multiphase flow refers to the complex physical phenomenon of mixing and flowing substances of different phases (such as gas-liquid and liquid-liquid) under compressible conditions. It is widely found in cutting-edge engineering fields such as fuel injection systems, high-speed underwater vehicles, and inertial confinement fusion. High-precision numerical simulation is a key means to study and predict such flows. Its core problem and technical difficulty lies in how to accurately capture the dynamically evolving phase interfaces accompanied by strong discontinuities.

[0003] To address the interface capture problem, two main categories of numerical methods have been developed: the Sharp Interface Method (SIM) and the Diffused Interface Method (DIM). The Sharp Interface Method treats the phase interface as a geometric discontinuity, employing moving meshes or interface tracing techniques for explicit characterization. While these methods offer advantages in interface clarity, they often face challenges such as difficult mesh reconstruction, high algorithmic complexity, and poor computational stability when dealing with complex topological changes like large deformations, fragmentation, and fusion of interfaces. The Diffused Interface Method, by introducing a physically or numerically finite thickness transition, implicitly captures the interface, avoiding complex interface topology management. This makes the algorithm simpler to implement and more suitable for complex flow simulations. Within the Diffused Interface Method, the five-equation model proposed by Allaire et al. is widely used due to its ability to effectively describe interphase mechanical equilibrium. However, traditional Diffused Interface Methods suffer from inherent numerical dissipation, causing the interface to gradually blur over time during simulation, limiting its long-term computational accuracy.

[0004] To improve interface resolution, researchers have developed high-order, high-precision numerical schemes. Among them, the Weighted Essentially Oscillated Non-Oscillation (WENO) scheme and the Discontinuous Galerkin (DG) method are the most representative. The WENO scheme maintains stability near discontinuities and high-order accuracy in smooth regions through a nonlinear weighted smooth template, but its large template size leads to low parallel computation efficiency, high memory overhead, and difficulty in implementation on complex unstructured meshes. The DG method achieves high-order accuracy by increasing the polynomial degree within the cells and has the advantage of a compact template, but its computational cost is high, time step is strictly limited, and robust limiter design remains a challenge for multiphase flow problems with strong discontinuities.

[0005] In recent years, the Multi-Moment Finite Volume Method (VPM) has attracted attention as a novel high-order method. This method achieves high-order reconstruction on a compact template by simultaneously defining two types of variables within the grid cells: the average value of the volume integral and the point values ​​at the cell vertices, effectively balancing computational accuracy and efficiency. However, existing VPM methods have a fundamental limitation: the updating of their point values ​​depends on solving a Riemann problem based on spatial derivatives, which necessitates the use of specific types of Riemann solvers (such as the linearized Roe scheme). This mechanism makes them incompatible with more robust and efficient Riemann solvers such as HLLC and AUSM series in multiphase flow simulations. This "compatibility" deficiency makes it difficult for the traditional VPM framework to effectively handle complex interphase interactions and high density ratios in compressible multiphase flows, greatly limiting its application in complex flow simulations with significant practical applications.

[0006] Therefore, developing a new general-purpose multi-moment finite volume method that can break through the limitations of traditional point value update mechanisms and is compatible with various Riemann solvers has become an urgent technical need to extend its application to high-precision simulation of compressible multiphase flows and promote the engineering application of this advanced numerical method. Summary of the Invention

[0007] The purpose of this invention is to provide a general multi-moment finite volume method for compressible multiphase flow in unstructured meshes, aiming to solve the key problem of incompatibility between the traditional multi-moment finite volume method and the general Riemann solver, and to achieve a smooth extension from single-phase flow to compressible multiphase flow simulation.

[0008] To achieve the above objectives, this invention provides a general multi-moment finite volume method for compressible multiphase flow in unstructured meshes. This method defines the average volume integral as a first-type variable on the grid cells of a computational mesh, defines point values ​​at the vertices of the grid cells as a second-type variable, and iteratively executes the following steps:

[0009] S1. Calculate the point value at the center of each grid cell based on the average value of the volume integral and the point value;

[0010] S2. Based on the nodes and adjacent mesh cells on the boundary surface, calculate the values ​​and gradients on each boundary surface using the least squares method;

[0011] S3. Calculate the first and second spatial derivatives within each grid cell and derive the coefficients of each term in the cell reconstruction polynomial;

[0012] S4. By forcing the average volume integral of the cell reconstruction polynomial on its corresponding grid cell to be equal to the current average volume integral, calculate the constant term coefficients of the cell reconstruction polynomial.

[0013] S5. Apply a MOGN limiter to the unit reconstruction polynomial;

[0014] S6. Apply physical correction conditions to the constrained unit reconstruction polynomial;

[0015] S7. At each boundary surface of each grid cell, the flux is calculated using a Riemann solver;

[0016] S8. Based on the calculated flux, update the volume integral average value using the integral form of the control equation;

[0017] S9. For each target vertex, construct a vertex reconstruction polynomial based on its reconstruction template using the least squares method.

[0018] S10. Apply a MOGN limiter to the vertex reconstruction polynomial;

[0019] S11. Apply physical correction conditions to the constrained vertex reconstruction polynomial;

[0020] S12. The flux is calculated at each point in the flux calculation template using a Riemann solver; wherein the flux calculation template is composed of the midpoints of the line segments connecting the target vertex and each point in the reconstruction template;

[0021] S13. Based on the flux value obtained on the flux calculation template, calculate the flux derivative at the target vertex using the least squares method;

[0022] S14. Update the point value of the target vertex using the differential control equation based on the flux derivative.

[0023] S15. Determine if the predetermined calculation time has been reached; if not, return to S1 to update the next time step; if it has been reached, end the simulation.

[0024] Optionally, the computational grid is any unstructured grid, including at least one grid cell type selected from triangles, quadrilaterals, tetrahedrons, hexahedrons, triangular prisms, or pyramids.

[0025] Optionally, in S3, the first and second derivative values ​​of the computational variables within each grid cell are obtained according to the Gaussian divergence theorem.

[0026] Optionally, in S3, the cell reconstruction polynomial of each grid cell is constructed according to the Taylor series, and the coefficients of the cell reconstruction polynomial are calculated.

[0027] Optionally, the reconstruction template of the target vertex includes the target vertex, the vertices directly adjacent to the target vertex, and the body center of the neighboring mesh cells.

[0028] Optionally, the Riemann solver is at least one of the Roe, HLL, and AUSM series.

[0029] Optionally, the physical protection conditions include: ensuring that the density and pressure are positive in single-phase flow simulation, and ensuring that the volume fraction is in the range of [0, 1] and the component density and internal energy are positive in multiphase flow simulation.

[0030] Optionally, the flux derivative calculation in S13 specifically includes: solving the flux polynomial coefficients by least squares approximation, and then obtaining the flux divergence at the target vertex.

[0031] Optionally, the general multi-moment finite volume method is time-stepped using a third-order Runge-Kutta method.

[0032] Based on the same inventive concept, the present invention also provides a readable storage medium having a computer program stored thereon, which, when executed, enables the implementation of the general multi-moment finite volume method for unstructured mesh compressible multiphase flow as described above.

[0033] The universal multi-moment finite volume method for compressible multiphase flow in unstructured meshes provided by this invention has at least one of the following advantages compared with existing multi-moment finite volume methods:

[0034] (1) It has broad compatibility with various Riemann solvers.

[0035] Existing multi-moment finite volume methods, in compressible applications, rely on derivative-based Riemannian solvers (such as the Roe scheme) for point-value flux calculations, which significantly limits the range of solvers that can be selected. This invention, however, redesigns the point-value flux calculation strategy, making the method compatible with various mainstream Riemannian solvers such as Roe, HLL series, and AUSM series. Users can flexibly select the most suitable Riemannian solver based on the characteristics of the specific simulation problem (such as single-phase / multi-phase, shock wave intensity, etc.).

[0036] (2) Successfully extended to compressible multiphase flow calculation

[0037] Because traditional derivative-based Riemann solvers are difficult to extend to multiphase flow models involving complex interactions, existing multi-moment finite volume methods are mainly limited to single-phase flow simulations in compressible applications. The point value update strategy proposed in this invention enables the direct use of Riemann solvers (such as HLLC) already maturely applied in the multiphase flow field, successfully extending the application scope of high-order multi-moment finite volume methods to the compressible multiphase flow domain, greatly expanding the applicability and practicality of the multi-moment finite volume method.

[0038] (3) The algorithm is unified, the code reuse rate is high, and the computational efficiency is improved.

[0039] In existing multi-moment finite volume methods for compressible applications, the flux calculation of the volume integral average and point values ​​requires two different Riemann solver codes, increasing the algorithm's complexity and implementation difficulty. In this invention, the flux calculation of the volume integral average and point values ​​uses the exact same Riemann solver, enabling the encapsulation of the Riemann solver code, improving code utilization and computational efficiency, and reducing algorithm complexity.

[0040] (4) It is compatible with unstructured meshes of arbitrary shapes and has stronger geometric adaptability.

[0041] Existing multi-moment finite volume methods are typically limited to handling regular-shaped mesh elements such as triangles, quadrilaterals, tetrahedrons, and hexahedrons in compressible applications. This invention, by employing a reconstruction strategy based on a global coordinate system and least squares, can be stably and accurately applied to arbitrary polygonal (two-dimensional) and polyhedral (three-dimensional) meshes. This significantly enhances its adaptability to complex computational geometry. Attached Figure Description

[0042] Those skilled in the art will understand that the accompanying drawings are provided to better understand the invention and do not constitute any limitation on the scope of the invention. Wherein:

[0043] Figure 1 A flowchart of a general multi-moment finite volume method for compressible multiphase flow in unstructured meshes provided in an embodiment of the present invention;

[0044] Figure 2 This is a schematic diagram illustrating various mesh types provided in one embodiment of the present invention;

[0045] Figure 3 A schematic diagram of a template for linear reconstruction of triangular and tetrahedral mesh cell surfaces provided in an embodiment of the present invention;

[0046] Figure 4 This is a schematic diagram of VIA evaluation points in a triangular mesh provided in an embodiment of the present invention;

[0047] Figure 5 A schematic diagram of the target vertex reconstruction template used for the hexagonal mesh and hexahedral mesh provided in an embodiment of the present invention;

[0048] Figure 6 Density schlieren images of a triple point problem at different times provided in an embodiment of the present invention;

[0049] Figure 7 This invention provides a density schlieren image at t=3s for processing the triple point problem using the method proposed in this invention, as provided in an embodiment of the invention.

[0050] Figure 8To process the density schlieren map of the triple point problem at t=3s using the existing third-order WENO-THINC method;

[0051] Figure 9 Density schlieren image at t=378μs for processing R22 bubble problem using the method proposed in this invention and the existing fifth-order WENO method, provided as an embodiment of the present invention;

[0052] Figure 10 Density schlieren image at t=378μs, provided as an embodiment of the present invention, for processing the R22 bubble problem using the method proposed in the present invention and the existing fifth-order MUSCL-THINC-BVD method;

[0053] Figure 11 This is a comparison diagram of numerical simulation results and experimental schlieren plots at different times for processing a three-dimensional helium bubble problem using the method proposed in this invention, provided as an embodiment of the present invention. Detailed Implementation

[0054] To make the objectives, technical solutions, and advantages of the present invention clearer, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0055] Therefore, the following detailed description of the embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0056] In the description of this invention, it should be understood that the terms "center," "upper," "lower," "left," "right," "vertical," "horizontal," "inner," and "outer," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings, or the orientation or positional relationship commonly used when the product is in use, or the orientation or positional relationship commonly understood by those skilled in the art. They are only used to facilitate the description of this invention and to simplify the description, and are not intended to indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on this invention.

[0057] Furthermore, relational terms such as "first" and "second" are used merely to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Moreover, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such an article or apparatus. Without further limitations, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the article or apparatus that includes said element. Those skilled in the art will understand the specific meaning of the above terms in this invention based on the specific circumstances.

[0058] Please refer to Figure 1 This invention provides a general multi-moment finite volume method for compressible multiphase flow in unstructured meshes. It can achieve third-order accuracy based on a compact template on any unstructured mesh, is compatible with various Riemann solvers, and can solve compressible flow problems involving shock waves, material interfaces, etc.

[0059] The generalized multi-moment finite volume method defines the average volume integral as a first-type variable on the grid cells of a computational mesh, and the point values ​​at the vertices of the grid cells as second-type variables, and iteratively executes the following steps:

[0060] S1. Calculate the point value at the volume center of each grid cell based on the average value of the volume integral and the point value;

[0061] S2. Based on the nodes and adjacent mesh cells on the boundary surface, calculate the values ​​and gradients on each boundary surface using the least squares method;

[0062] S3. Calculate the first and second spatial derivatives within each grid cell and derive the coefficients of each term in the cell reconstruction polynomial;

[0063] S4. By forcing the average volume integral of the element reconstruction polynomial on its corresponding mesh element to be equal to the current average volume integral, calculate the constant term coefficients of the element reconstruction polynomial.

[0064] S5. Apply a MOGN limiter to the unit reconstruction polynomial;

[0065] S6. Apply physical correction conditions to the constrained unit reconstruction polynomial;

[0066] S7. At each boundary surface of each grid cell, the flux is calculated using a Riemann solver;

[0067] S8. Based on the calculated flux, update the volume integral average value using the integral form of the control equation;

[0068] S9. For each target vertex, construct a vertex reconstruction polynomial based on its reconstruction template using the least squares method.

[0069] S10. Apply a MOGN limiter to the vertex reconstruction polynomial;

[0070] S11. Apply physical correction conditions to the constrained vertex reconstruction polynomial;

[0071] S12. The flux is calculated at each point in the flux calculation template using a Riemann solver; wherein, the flux calculation template is composed of the midpoints of the line segments connecting the target vertex and each point in the reconstruction template;

[0072] S13. Based on the flux value obtained from the flux calculation template, calculate the flux derivative at the target vertex using the least squares method.

[0073] S14. Update the point value of the target vertex using the differential control equation based on the flux derivative.

[0074] S15. Determine if the predetermined calculation time has been reached; if not, return to S1 to update the next time step; if it has been reached, end the simulation.

[0075] This invention is applicable to various mesh types, including two-dimensional triangles, quadrilaterals, pentagons, and arbitrary polygons, and three-dimensional tetrahedrals, hexahedrals, triangular prisms, pyramids, and arbitrary polyhedral meshes, such as... Figure 2 As shown in the figure, This represents the control volume, where I represents the total number of grid cells. Control Volume boundary surface It indicates that the control body Points used Indicates. Control body body and mind It means, face The filling of the dough express.

[0076] At the same time, an indexing system for nodes is defined. Represents a node, where This represents the total number of nodes. Directly connected nodes express, The surrounding grid is used The multi-moment finite volume method simultaneously calculates and updates the volume integral average (VIAs) and point values ​​(point values ​​s) on each mesh element, using... If the variables to be calculated are defined, then the average value of the volume integral and the point value are defined as follows:

[0077]

[0078] The essence of this method is to realize the flow field evolution by iteratively executing the following two core processes: the volume integral average value update process (S1-S8 above) and the point value update process (S9-S14 above).

[0079] The update process for the average volume integral is as follows:

[0080] First, execute S1 to calculate the point value at the center of each mesh element based on the average value of the volume integral and the point value. For two-dimensional triangular and quadrilateral meshes, and three-dimensional tetrahedral, hexahedral, triangular prism, and pyramid meshes, the calculation is based on the conversion relationship between the center value and the average value of the volume integral and the point value. This embodiment provides the conversion relationship for different mesh element types as shown in Table 1 below.

[0081] Table 1

[0082]

[0083] in, Indicates the variable in the grid cell The value at the body center, Represents grid cells The average value of the volume integral, Represents grid cells Values ​​at each point.

[0084] For any other two-dimensional polygon and any other three-dimensional polyhedron, the volume center value of the target mesh is calculated using the least squares method based on the average volume integral and the point values. Specifically, a quadratic polynomial is assumed to satisfy the following condition:

[0085]

[0086] The volume integral is obtained through the third-order symmetric integration point. This yields an overdetermined system of equations, which can be written in matrix form as follows:

[0087]

[0088] Where M is a matrix that depends only on geometric position, A is the unknown matrix to be solved, and B is the variable matrix. A is solved using the following formula:

[0089]

[0090] The matrix inversion process is calculated using the SVD method.

[0091] Next, execute S2, based on the boundary surface. The vertices (nodes) and adjacent mesh cells on the boundary are used to calculate the values ​​and gradients on each boundary surface using the least squares method. Please refer to [reference needed]. Figure 3 Using two-dimensional triangular meshes and three-dimensional tetrahedral meshes as examples, a template for surface linear reconstruction is shown.

[0092] The linear interpolation polynomial is:

[0093]

[0094] in, Indicates the variable to be calculated. Indicates the variable in the cell The The value at each face. The coordinates are the coordinates of the face center.

[0095] The system of equations formed by substituting the template can be written in matrix form:

[0096]

[0097] It is worth mentioning that this step can be performed before the time loop, meaning that the linear least squares matrix for each face only needs to be calculated once and saved for later use.

[0098] Solve the linear least-squares overdetermined equations for each surface to obtain the average value and first derivative value of the calculated variables for each boundary surface. Substitute the values ​​of the constituent template into the resulting overdetermined equation system and solve it using the following formula:

[0099]

[0100] Then execute S3 to calculate the first and second spatial derivatives within each grid cell and derive the coefficients of each term in the cell reconstruction polynomial;

[0101] In this embodiment, the first and second derivatives of the computational variables within each grid cell are obtained according to the Gaussian divergence theorem. Specifically, the values ​​are as follows:

[0102]

[0103]

[0104] in, Representation unit volume, Representation unit No. The area of ​​each face. This represents the out-of-plane normal vector.

[0105] Based on the Taylor series, construct the element reconstruction polynomial for each mesh element and calculate its coefficients. The element reconstruction polynomial is in the form of:

[0106]

[0107] Among them, except The calculation methods for the external coefficients are as follows:

[0108]

[0109] Next, S4 is executed, which forces the average volume integral of the element reconstruction polynomial on its corresponding mesh element to be equal to the current average volume integral, and then calculates the coefficients of the constant term of the element reconstruction polynomial. That is, the body's core value. In other words:

[0110]

[0111] Once this step is completed, the complete element reconstruction polynomial for the average volume integral will be obtained.

[0112] Then, S5 is executed, applying MOGN restrictors to the element reconstruction polynomials. The MOGN higher-order restrictor combines the MOGV restrictor with the restriction functions proposed by Nishikawa. A restriction function is defined for each element. The calculation method is as follows:

[0113]

[0114] in, Indicates the first The slope limiter is measured at [number] assessment points. The total number of these assessment points is [number]. This includes all Gaussian points on all boundary surfaces of the mesh, as well as the midpoint of the line connecting the mesh volume center and a vertex. Taking a triangular mesh as an example, all evaluation points are as follows: Figure 4 As shown.

[0115] The limiter function calculation method for each evaluation point is as follows:

[0116]

[0117] in, and This represents the maximum and minimum values ​​of the computed variables in the reconstruction template of the current mesh cell. This represents the value of the computed variable at the evaluation point of the polynomial computation reconstructed using unrestricted units in S4.

[0118] To maintain high-order convergence, the limiter needs to be deactivated in smooth regions, and the limiter needs further modification:

[0119]

[0120] in, ,and Indicates the feature grid size.

[0121] Applying a limiter to the element reconstruction polynomial yields the following restricted element reconstruction polynomial:

[0122]

[0123] Next, S6 is executed to apply physical protection conditions to the constrained reconstructed polynomials. These physical protection conditions include protection conditions for the Euler equations of single-phase flow and guarantee boundary conditions for the five-equation model of two-phase flow, applied to the evaluation points mentioned in S5. For the Euler equations of single-phase flow, it is assumed that the reconstruction variables at each evaluation point after S5 constraints are... To ensure that the pressure and density are positive, a small quantity is introduced. Then it is necessary to ensure , .

[0124] First, adjust the density:

[0125]

[0126] in,

[0127]

[0128] Next, adjust the pressure:

[0129]

[0130] in,

[0131]

[0132] The corrected reconstructed variables are This ensures that the pressure and density are positive numbers, which will be used as inputs to the Riemann solver in subsequent flux calculation steps.

[0133] For the five-equation model of two-phase flow, assume that the reconstructed variables at each evaluation point after S5 constraint are: To ensure the volume fraction Within the range [0, 1], while ensuring component density and internal energy It is a positive number.

[0134] First, adjust the volume fraction:

[0135]

[0136] in,

[0137]

[0138] Next, adjust the component density:

[0139]

[0140] in,

[0141]

[0142] The variables obtained after these two steps of correction are:

[0143] ,

[0144] Next, we adjust the internal energy. Let... ,make sure Correct and refactor variables

[0145]

[0146] in,

[0147]

[0148] Then, S7 is executed, and the flux is calculated at each boundary surface of each grid cell using a Riemann solver. The governing equations are:

[0149]

[0150] Where U is the vector of computational variables, Let S be the flux vector and S be the source term vector.

[0151] In the Euler equations for single-phase flow,

[0152]

[0153] In the five-equation model of two-phase flow,

[0154]

[0155] When calculating the flux at each surface, various Riemann solvers such as Roe, HLL, and AUSM can be used, namely:

[0156]

[0157] The left side of the equation represents the dot product of flux and surface flux at the boundary surface, while the right side represents the Riemann solver used. The left and right values ​​of the input are obtained by the element reconstruction polynomial calculation of the upstream and downstream meshes.

[0158] Then, S8 is executed. Based on the calculated flux, the volume integral average value is updated using the integral form of the governing equation, as shown in the following formula:

[0159]

[0160] The above steps describe the operations performed to update the average value of the volume integral. Next, we will process the point value updates, as follows:

[0161] First, execute S9. For each target vertex, construct a vertex reconstruction polynomial using the least squares method based on its reconstruction template. Assume the vertex reconstruction polynomial at each point is:

[0162]

[0163] The least squares method based on compact templates is used to solve for the coefficients of the element reconstruction polynomial. Figure 5 The reconstruction templates for the target vertices are illustrated using two-dimensional hexagonal meshes and three-dimensional hexahedral meshes as examples. These templates include the target vertex, vertices directly adjacent to the target vertex, and the volume centers of neighboring mesh cells. Substituting the coordinates and values ​​of the reconstruction template points into the vertex reconstruction polynomial yields an overdetermined system of equations. Solving this system provides the coefficients of the vertex reconstruction polynomial.

[0164] Next, S10 is executed, applying a MOGN limiter to the vertex reconstruction polynomial. A limit function is defined for each target vertex. The calculation method is as follows:

[0165]

[0166] in, Indicates the first The slope limiter is measured at [number] assessment points. The total number of these assessment points is [number]. This includes the midpoint of the line connecting the target vertex to its adjacent vertex, and the midpoint of the line connecting the target vertex to the center of the adjacent mesh.

[0167] The limiter function calculation method for each evaluation point is as follows:

[0168]

[0169] in, and This represents the maximum and minimum computed variable values ​​in the reconstruction template of the current target vertex. This represents the value of the computed variable at the evaluation point of the polynomial computation reconstructed using unrestricted units in S9.

[0170] To maintain high-order convergence, the limiter needs to be deactivated in smooth regions. Similar to the treatment of the volume integral averaging, the limiter is further modified to obtain the limiter function. .

[0171] Applying a limiter to the unit reconstructed polynomial yields the restricted unit reconstructed polynomial:

[0172]

[0173] Compared to the constraint of the volume integral average, point values ​​require an additional step to restrict the point values ​​to between the maximum and minimum values ​​of the surrounding VIAs to improve computational stability:

[0174]

[0175] Next, S11 is executed to apply physical protection conditions to the constrained vertex reconstruction polynomial. The physical protection conditions include protection conditions for the Euler equations for single-phase flow and guarantee boundary conditions for the five-equation model for two-phase flow, which are applied to the evaluation points mentioned in S10.

[0176] For the Euler equations for single-phase flow, assume that the reconstructed variables at each evaluation point after the S10 constraint are... To ensure that the pressure and density are positive, a small quantity is introduced. Then it is necessary to ensure , .

[0177] First, adjust the density:

[0178]

[0179] in,

[0180]

[0181] Next, adjust the pressure:

[0182]

[0183] in,

[0184]

[0185] The corrected reconstructed variables are This ensures that the pressure and density are positive numbers, which will be used as inputs to the Riemann solver in subsequent flux calculation steps.

[0186] For the five-equation model of two-phase flow, assuming the reconstructed variables at each evaluation point after S10 constraint are: To ensure the volume fraction Within [0,1], while ensuring component density and internal energy It is a positive number.

[0187] First, adjust the volume fraction:

[0188]

[0189] in,

[0190]

[0191] Next, adjust the component density:

[0192]

[0193] in,

[0194]

[0195] The variables obtained after these two steps of correction are:

[0196] ,

[0197] Next, we adjust the internal energy. Let... ,make sure Correct and refactor variables:

[0198]

[0199] in,

[0200]

[0201] Then, step S12 is executed, calculating the flux at each point in the flux calculation template using a Riemann solver. The flux calculation template is composed of the midpoints of the line segments connecting the target vertex to each point in the reconstruction template; more precisely, the midpoints of the lines connecting the target vertex to adjacent vertices, and the midpoints of the lines connecting the target vertex to the center of the adjacent mesh. Various Riemann solvers can be used to calculate the flux at the points in the flux calculation template, including Roe, HLL, and AUSM-type Riemann solvers.

[0202]

[0203] Here, the subscript mh represents the h-th flux calculation template point at the m-th point. The left side of the equation is the dot product of the flux and the direction vector at the flux calculation template point, and the right side is the Riemann solver used. The left-hand side of the input is calculated by the element reconstruction polynomial of the target vertex, and the right-hand side is calculated by the element reconstruction polynomial of adjacent vertices or adjacent meshes. The direction vector here will be defined in the derivation in S13 below.

[0204] Next, step S13 is executed. Based on the flux value obtained from the flux calculation template, the flux derivative at the target vertex is calculated using the least squares method. Specifically, after calculating the flux from the flux calculation template, the flux at the target point is calculated using the least squares method. Assume a quadratic polynomial distribution for the flux:

[0205]

[0206] Substituting the coordinates and flux values ​​of each point in the flux calculation template into the reconstructed polynomial of the above unit yields an overdetermined system of equations. Solving this system yields the coefficients. The coefficients of the three linear terms can be written as:

[0207]

[0208] In the formula, , and It only depends on the location of the flux calculation template point and can be calculated before the time loop.

[0209] Comparing the Taylor expansion, we can see that the coefficients of these three linear terms correspond to the first derivatives, therefore:

[0210]

[0211] Then the divergence of flux can be written as:

[0212]

[0213] in,

[0214]

[0215] This is the definition of the direction vector in S12, which can be calculated and stored before the time loop, and does not need to be calculated repeatedly in each update.

[0216] Then, S14 is executed, and the point value of the target vertex is updated according to the flux derivative using the differential governing equation, which is as follows:

[0217]

[0218] Finally, step S15 is executed to determine if the predetermined computation time has been reached. If not, the process returns to S1 to update the next time step; if it has been reached, the simulation ends. In this embodiment, time stepping is performed using the third-order Runge-Kutta method. The third-order Runge-Kutta formula is as follows:

[0219]

[0220] Specifically, regarding the update of the average volume integral, Let S8 be the derivative of the average volume integral with respect to time; for updating point values, The derivative of the point value obtained for S14 with respect to time.

[0221] The universal multi-moment finite volume method (VPM-GF) for compressible multiphase flow in unstructured meshes proposed in this invention differs from existing multi-moment finite volume methods in that:

[0222] (1) The point value reconstruction methods are different. In most existing multi-moment finite volume methods, the reconstruction of point values ​​requires projection onto a local coordinate system and the least squares method based on the volume average value, node values, and volume center derivative to achieve quadratic polynomial reconstruction. However, in the method proposed in this invention, the point value reconstruction is achieved by the least squares method based on the vertex values ​​directly connected to the target vertex and the volume center values ​​of adjacent meshes in the global coordinate system.

[0223] (2) The flux calculation methods for point values ​​differ. In the previous multi-moment finite volume method, the flux of point values ​​was calculated using a derivative-based Roe Riemann solver, and the flux divergence was transformed into the dot product of the Jacobian matrix and the derivative of the variable using the chain rule. However, in the method proposed in this invention, the flux at the flux calculation template point needs to be calculated first. The flux calculation template includes the midpoint of the line connecting the target vertex and its directly adjacent vertex, as well as the midpoint of the line connecting the target vertex and the center point of the adjacent grid. Then, the flux of the target point is calculated using the least squares method.

[0224] (3) Since the derivative-based Roe Riemann solver is difficult to extend to multiphase flow, the existing multi-moment finite volume method is limited to single-phase flow in compressible applications. However, the Riemann solver proposed in this invention is more simply and directly extended to multiphase flow, so the proposed method can be applied to compressible multiphase flow.

[0225] (4) In the previous point value reconstruction process, no gradient limiter and correction condition were applied, while the method proposed in this invention also adds MOGN limiter and physical correction condition to the point value reconstruction.

[0226] Therefore, the universal multi-moment finite volume method for compressible multiphase flow in unstructured meshes proposed in this invention has at least the following advantages compared with existing multi-moment finite volume methods:

[0227] (1) It has broad compatibility with various Riemann solvers.

[0228] Existing multi-moment finite volume methods, in compressible applications, rely on derivative-based Riemannian solvers (such as the Roe scheme) for point-value flux calculations, which significantly limits the range of solvers that can be selected. This invention, however, redesigns the point-value flux calculation strategy, making the method compatible with various mainstream Riemannian solvers such as Roe, HLL series, and AUSM series. Users can flexibly select the most suitable Riemannian solver based on the characteristics of the specific simulation problem (such as single-phase / multi-phase, shock wave intensity, etc.).

[0229] (2) Successfully extended to compressible multiphase flow calculation

[0230] Because traditional derivative-based Riemann solvers are difficult to extend to multiphase flow models involving complex interphase interactions, existing multi-moment finite volume methods are mainly limited to single-phase flow simulations in compressible applications. The point value update strategy proposed in this invention enables the direct use of Riemann solvers (such as HLLC) already maturely applied in the multiphase flow field, successfully extending the application scope of high-order multi-moment finite volume methods to the compressible multiphase flow domain, greatly expanding the applicability and practicality of the multi-moment finite volume method.

[0231] (3) The algorithm is unified, the code reuse rate is high, and the computational efficiency is improved.

[0232] In existing multi-moment finite volume methods for compressible applications, the flux calculation of the volume integral average and point values ​​requires two different Riemann solver codes, increasing the algorithm's complexity and implementation difficulty. In this invention, the flux calculation of the volume integral average and point values ​​uses the exact same Riemann solver, enabling the encapsulation of the Riemann solver code, improving code utilization and computational efficiency, and reducing algorithm complexity.

[0233] (4) It is compatible with unstructured meshes of arbitrary shapes and has stronger geometric adaptability.

[0234] Existing multi-moment finite volume methods are typically limited to handling regular-shaped mesh elements such as triangles, quadrilaterals, tetrahedrons, and hexahedrons in compressible applications. This invention, by employing a reconstruction strategy based on a global coordinate system and least squares, can be stably and accurately applied to arbitrary polygonal (two-dimensional) and polyhedral (three-dimensional) meshes. This significantly enhances its adaptability to complex computational geometry.

[0235] The following examples demonstrate the superior performance of the proposed Universal Multi-Move Finite Volume Method (VPM-GF).

[0236] (1) Euler equation solution test: used to demonstrate the numerical error and convergence order of the method in solving the Euler equation, and to prove the compatibility of the method with various Riemann solvers.

[0237] Table 2

[0238]

[0239] Table 2 shows the solutions used with three widely used Riemannian solvers: Roe, AUSM, etc. + The calculation results for -up and HLLC demonstrate the compatibility of the proposed method with various Riemann solvers. Tests were conducted on quadrilateral, triangular, and polygonal meshes, proving the compatibility of the proposed method with arbitrary polygons. Furthermore, the results show that the proposed method can achieve third-order convergence accuracy.

[0240] (2) Compressible three-phase point problem: used to demonstrate the good performance of the method proposed in this invention in compressible two-phase flow problem.

[0241] Please refer to Figure 6 , Figure 6 For density schlieren images of the triple point problem at different times, the method proposed in this invention successfully captured the complete temporal evolution of all key flow features: the initial contact discontinuity, the leftward propagating rarefaction wave and the rightward moving shock wave, and then reproduced the formation of vortices near the triple point, the shear layer developing along the contact discontinuity, and the complex shock wave reflection and interaction patterns. Figure 6 The sharp interface and the clear small-scale structure resulting from the Kelvin–Helmholtz instability demonstrate the excellent simulation performance of the proposed method.

[0242] Please refer to Figures 7-8 , Figures 7-8 To address the density schlieren and density contour plots of the triple point problem on a refined mesh using the method proposed in this invention, and to compare it with the existing third-order WENO-THINC method... [1] The calculation results were compared. With the increase in mesh resolution, more and more small-scale structures were reproduced, and the interface capture was sufficiently clear. In contrast, the method proposed in this invention significantly achieved higher resolution and preserved more refined vortex structures along contact discontinuities. This indicates that the method proposed in this invention can reproduce complex flow field structures with significantly improved resolution and successfully capture minute vortex features.

[0243] (3) The problem of the interaction between shock waves in the air and R22 bubbles: This is used to demonstrate the good performance of the proposed method in compressible two-phase flow problems.

[0244] Please refer to Figures 9-10 , Figures 9-10 The method proposed in this invention processes the density schlieren plot of the R22 bubble problem at t=378μs, and compares it with the fifth-order WENO and MUSCL-THINC-BVD methods. [2] The results were compared, among which... Figure 9 and Figure 10 The upper half of the diagrams shows density schlieren images using existing fifth-order WENO and MUSCL-THINC-BVD methods, while the lower half shows density schlieren images using the method proposed in this invention. A significant characteristic of this problem is the instability that develops along the interface, followed by interface roll-up and the generation of small-scale structures. The method proposed in this invention achieves an interface thickness comparable to the fifth-order WENO scheme and better preserves the vortices generated by the baroclinic mechanism, exhibiting a more complex flow morphology than the results of the fifth-order WENO scheme. Simultaneously, the method proposed in this invention captures reflected and transmitted shock waves with clear and compact profiles, demonstrating excellent numerical stability. Although the material interface generated by the method proposed in this invention is slightly thicker than that of the MUSCL-THINC-BVD method, it achieves a comparable level of accuracy in resolving small-scale structures and roll-up modes.

[0245] (4) Interaction between three-dimensional air shock wave and helium bubble: This is used to demonstrate the good performance of the proposed method in compressible two-phase flow problems.

[0246] Please refer to Figure 11 , Figure 11 Numerical simulation results and experimental schlieren plots of the three-dimensional helium bubble problem at different times [3] The comparison shows that, due to the different gas properties of helium and air, the propagation speed of the transmitted shock wave within the helium bubble is faster than that of the incident shock wave in the air, resulting in the complex flow structure shown in the schlieren image. The high-speed air jet penetrates the helium bubble, forming a unique double-ring structure, sometimes referred to as a double-ring vortex. These two vortex rings gradually separate as convection occurs at different speeds, with the downstream vortex ring becoming larger than the upstream one. The method proposed in this invention reproduces the entire process of the formation, evolution, and separation of this double-ring system with sufficient accuracy, matching the experimental results.

[0247] The references are as follows:

[0248] [1]Paula T, Adami S, Adams N A. A robust high-resolution discrete-equations method for compressible multi-phase flow with accurate interfacecapturing[J]. Journal of Computational Physics, 2023, 491: 112371;

[0249] [2] Deng X, Inaba S, Xie B, et al. High fidelity discontinuity-resolving reconstruction for compressible multiphase flows with moving interfaces[J]. Journal of Computational Physics, 2018, 371: 945-966;

[0250] [3]Layes G, Le Métayer O. Quantitative numerical and experimental studies of the shock accelerated heterogeneous bubbles motion[J]. Physics of Fluids, 2007, 19(4).

[0251] Based on the same inventive concept, embodiments of the present invention also propose a readable storage medium storing a computer program that, when executed, can implement the general multi-moment finite volume method for unstructured mesh compressible multiphase flow as described above.

[0252] A readable storage medium can be a tangible device capable of holding and storing instructions for use by an instruction execution device, such as, but not limited to, electrical storage devices, magnetic storage devices, optical storage devices, electromagnetic storage devices, semiconductor storage devices, or any suitable combination thereof. More specific examples of readable storage media (a non-exhaustive list) include: portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), static random access memory (SRAM), portable compact disc read-only memory (CD-ROM), digital multifunction disc (DVD), memory sticks, floppy disks, mechanical encoding devices, such as punch cards or recessed protrusions storing instructions thereon, and any suitable combination thereof. The computer programs described herein can be downloaded from the readable storage medium to various computing / processing devices, or downloaded via a network, such as the Internet, local area network, wide area network, and / or wireless network, to an external computer or external storage device. Networks can include copper transmission cables, fiber optic transmissions, wireless transmissions, routers, firewalls, switches, gateway computers, and / or edge servers. Each computing / processing device's network adapter card or network interface receives and forwards a computer program from the network for storage on a readable storage medium within the respective computing / processing device. The computer program used to perform the operations of this invention can be execution instructions, instruction set architecture (ISA) instructions, machine instructions, machine-dependent instructions, microcode, firmware instructions, status setting data, or source code or object code written in any combination of one or more programming languages, including object-oriented programming languages ​​such as Smalltalk, C++, etc., and conventional procedural programming languages ​​such as "C" or similar languages. The computer program can execute entirely on the user's computer, partially on the user's computer, as a standalone software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In cases involving remote computers, the remote computer can be connected to the user's computer via any type of network, including a local area network (LAN) or a wide area network (WAN), or it can be connected to an external computer (e.g., via the Internet using an Internet service provider). In some embodiments, electronic circuits, such as programmable logic circuits, field-programmable gate arrays (FPGAs), or programmable logic arrays (PLAs), are personalized by utilizing state information from a computer program. These electronic circuits can execute computer-readable program instructions, thereby realizing various aspects of the present invention.

[0253] Various aspects of the present invention are described herein with reference to flowchart illustrations and / or block diagrams of methods, systems, and computer program products according to embodiments of the invention. It should be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by a computer program. These computer programs can be provided to a processor of a general-purpose computer, a special-purpose computer, or other programmable data processing apparatus to produce a machine such that, when executed by the processor of the computer or other programmable data processing apparatus, they create means for implementing the functions / actions specified in one or more blocks of the flowchart illustrations and / or block diagrams. These computer programs can also be stored in a readable storage medium that causes a computer, programmable data processing apparatus, and / or other device to operate in a particular manner; thus, the readable storage medium storing the computer program comprises an article of manufacture including instructions for implementing aspects of the functions / actions specified in one or more blocks of the flowchart illustrations and / or block diagrams.

[0254] A computer program may also be loaded onto a computer, other programmable data processing apparatus, or other device to cause a series of operational steps to be performed on the computer, other programmable data processing apparatus, or other device to produce a computer-implemented process, thereby causing the computer program executing on the computer, other programmable data processing apparatus, or other device to perform the functions / actions specified in one or more boxes of a flowchart and / or block diagram.

[0255] Since the readable storage medium provided by this invention belongs to the same inventive concept as the general multi-moment finite volume method for compressible multiphase flow in unstructured grids described above, the readable storage medium provided by this invention has all the advantages of the general multi-moment finite volume method for compressible multiphase flow in unstructured grids described above. Therefore, the beneficial effects of the readable storage medium provided by this invention will not be elaborated here.

[0256] The above description is merely a description of preferred embodiments of the present invention and is not intended to limit the scope of the invention in any way. Any changes or modifications made by those skilled in the art based on the above disclosure are within the protection scope of the present invention. Obviously, those skilled in the art can make various modifications and variations to the present invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the present invention and its equivalents, the present invention also intends to include these modifications and variations.

Claims

1. A general multi-moment finite volume method for compressible multiphase flow in unstructured meshes, characterized in that, The general multi-moment finite volume method defines the average volume integral as a first-type variable on a computational grid cell, and defines point values ​​at the vertices of the grid cells as a second-type variable, and iteratively executes the following steps: S1. Calculate the point value at the center of each grid cell based on the average value of the volume integral and the point value; S2. Based on the nodes and adjacent mesh cells on the boundary surface, calculate the values ​​and gradients on each boundary surface using the least squares method; S3. Calculate the first and second spatial derivatives within each grid cell and derive the coefficients of each term in the cell reconstruction polynomial; S4. By forcing the average volume integral of the cell reconstruction polynomial on its corresponding grid cell to be equal to the current average volume integral, calculate the constant term coefficients of the cell reconstruction polynomial. S5. Apply a MOGN limiter to the unit reconstruction polynomial; S6. Apply physical correction conditions to the constrained unit reconstruction polynomial; S7. At each boundary surface of each grid cell, the flux is calculated using a Riemann solver; S8. Based on the calculated flux, update the volume integral average value using the integral form of the control equation; S9. For each target vertex, construct a vertex reconstruction polynomial based on its reconstruction template using the least squares method. S10. Apply a MOGN limiter to the vertex reconstruction polynomial; S11. Apply physical correction conditions to the constrained vertex reconstruction polynomial; S12. The flux is calculated at each point in the flux calculation template using a Riemann solver; wherein the flux calculation template is composed of the midpoints of the line segments connecting the target vertex and each point in the reconstruction template; S13. Based on the flux value obtained on the flux calculation template, calculate the flux derivative at the target vertex using the least squares method; S14. Update the point value of the target vertex using the differential control equation based on the flux derivative. S15. Determine if the predetermined calculation time has been reached; if not, return to S1 to update the next time step; if it has been reached, end the simulation.

2. The universal multi-moment finite volume method for compressible multiphase flow in unstructured meshes according to claim 1, characterized in that, The computational grid is an unstructured grid, including at least one grid cell type selected from triangles, quadrilaterals, tetrahedrons, hexahedrons, triangular prisms, or pyramids.

3. The universal multi-moment finite volume method for compressible multiphase flow in unstructured meshes according to claim 1, characterized in that, In S3, the first and second derivative values ​​of the computational variables within each grid cell are obtained according to the Gaussian divergence theorem.

4. The universal multi-moment finite volume method for compressible multiphase flow in unstructured meshes according to claim 1, characterized in that, In S3, the cell reconstruction polynomial of each grid cell is constructed according to the Taylor series, and the coefficients of the cell reconstruction polynomial are calculated.

5. The universal multi-moment finite volume method for compressible multiphase flow in unstructured meshes according to claim 1, characterized in that, The reconstruction template of the target vertex includes the target vertex, the vertices directly adjacent to the target vertex, and the body center of the neighboring mesh cells.

6. The universal multi-moment finite volume method for compressible multiphase flow in unstructured meshes according to claim 1, characterized in that, The Riemann solver is at least one of the Roe, HLL, and AUSM series.

7. The universal multi-moment finite volume method for compressible multiphase flow in unstructured meshes according to claim 1, characterized in that, The physical protection conditions include: ensuring that the density and pressure are positive in single-phase flow simulation, and ensuring that the volume fraction is in the range of [0, 1] and that the component density and internal energy are positive in multiphase flow simulation.

8. The universal multi-moment finite volume method for compressible multiphase flow in unstructured meshes according to claim 1, characterized in that, The flux derivative calculation in S13 specifically includes: solving the flux polynomial coefficients by least squares approximation, and then obtaining the flux divergence at the target vertex.

9. The universal multi-moment finite volume method for compressible multiphase flow in unstructured meshes according to claim 1, characterized in that, The generalized multi-moment finite volume method is time-stepped using the third-order Runge-Kutta method.

10. A readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed, it can implement the general multi-moment finite volume method for compressible multiphase flow in unstructured grids according to any one of claims 1-9.

Images