Coupling level set and fluid volume high-order interface capturing method based on mixing precision

Abstract

The invention relates to a high-order interface capturing method for a coupling level set and a fluid volume based on mixing precision, which comprises the following steps of: firstly, dividing a calculation area into a grid comprising a plurality of control body units, initializing, then identifying interface units, and reconstructing a THINC interface function by adopting a high-order polynomial based on an LS value; thirdly, introducing a mixed precision iterative algorithm to solve an interface position nonlinear equation, and adaptively switching to multi-precision operation to obtain a robust solution when double precision cannot converge; then, flux is calculated through an accurate interface, a VOF transport equation is solved, and mass conservation is guaranteed; and finally, on the basis of the predicted interface geometry, an LS field which does not need to be reinitialized is directly reconstructed through pure geometric methods such as interface point cloud generation through projection and nearest point searching in a narrow band. According to the method, the capture precision and geometric fidelity of complex interfaces with high curvature, sharp corners and the like on the unstructured grids are remarkably improved, and the method has high precision and strong robustness.

Description

Technical Field

[0001] This invention relates to the field of computational fluid dynamics, and in particular to a high-order interface capture method based on hybrid precision of coupled level sets and fluid volume. Background Technology

[0002] In computational fluid dynamics, accurately simulating multiphase flows involving dynamic interfaces (such as gas-liquid flows, droplet collisions, and spray processes) is a core challenge for numerous engineering and scientific applications. Interface capture techniques, by implicitly representing interfaces within the Eulerian framework, have become a key tool for solving such problems, aiming to accurately track the interface location, geometric features, and their evolution in complex flow fields.

[0003] Currently, existing interface capture methods mainly include the Volume of Fluid (VOF) method, the Level Set (LS) method, and the coupled VOF and LS method (i.e., the CLSVOF method).

[0004] The Volume-of-Flight (VOF) method represents the two-phase distribution using a volume fraction field, and its greatest advantage lies in strictly guaranteeing mass conservation. Among them, the geometric reconstruction-based VOF method explicitly reconstructs the interface through techniques such as piecewise linear interface calculation, reducing numerical dissipation and becoming the most widely used method. However, traditional geometric VOF methods are usually based on the assumption of linear or planar interfaces, which can lead to significant geometric errors when dealing with interfaces with high curvature or sharp features. Although some studies have proposed surface reconstruction algorithms to improve accuracy, these methods are generally computationally complex and mainly limited to structured meshes, making it difficult to extend to unstructured meshes.

[0005] The LS method characterizes the interface by introducing a smooth signed distance function, which possesses favorable mathematical properties, facilitating high-precision calculation of geometric quantities such as interface normals and curvature. However, the properties of the signed distance function are compromised after numerically solving the transport equations, requiring an additional "reinitialization" process to restore them. This step introduces numerical errors, leading to interface offset and mass non-conservation, and is computationally time-consuming and highly parameter-dependent.

[0006] To balance mass conservation and geometric accuracy, the method of coupling Level Set and Volume of Fluid (VOF) has emerged, known as the Coupled Level Set and Volume of Fluid (CLSVOF) method. Existing CLSVOF methods can be divided into two categories based on their VOF components: algebraic VOF and geometric VOF. While algebraic CLSVOF methods are easy to implement on unstructured meshes, their LS fields are constructed using algebraic formulas, resulting in insufficient geometric accuracy and parameter sensitivity. Geometric VOF-based CLSVOF methods can perform explicit interface reconstruction, but the reconstructed interfaces are mostly linear, making it difficult to represent high-order surfaces. They also face the challenge of complex and difficult-to-implement geometric operations on unstructured meshes. Although recent research has attempted to introduce curve or surface interface reconstruction to improve accuracy, its application remains limited to structured meshes, and the fundamental problem of achieving high-order surface representation on general unstructured meshes has not yet been solved.

[0007] Furthermore, the strategy of coupling THINC functions with LS provides a new approach for constructing higher-order interface representations. The THINC method constructs interfaces using functions with analytical expressions, naturally supporting extension to higher-order polynomials. However, existing THINC-LS methods have not yet broken free from the evolutionary framework of traditional LS methods, still relying on complex numerical solutions of partial differential equations and re-initialization processes to maintain LS field properties. The algorithms are complex and difficult to extend to unstructured meshes. Meanwhile, how to efficiently and robustly construct geometric LS fields based on higher-order THINC interfaces, and how to ensure the convergence of the interface position equations under extreme parameters, remain critical technical bottlenecks that need to be addressed.

[0008] In summary, existing interface capture technologies generally face the following challenges in pursuing high precision and high fidelity:

[0009] 1) It is difficult to represent high-order surface interfaces on unstructured meshes;

[0010] 2) It is difficult to achieve both mass conservation and high geometric accuracy simultaneously;

[0011] 3) The re-initialization process of traditional level set methods leads to quality loss and increases computational complexity;

[0012] 4) Existing high-order methods are susceptible to rounding errors when solving for interface positions, resulting in convergence issues.

[0013] Therefore, developing a high-order interface capture method that can accurately capture complex interface geometry, strictly guarantee mass conservation, be applicable to unstructured meshes, and be computationally robust has significant theoretical value and engineering application requirements. Summary of the Invention

[0014] The purpose of this invention is to provide a high-order interface capture method based on hybrid precision coupled level set and fluid volume, so as to improve the numerical convergence of interface position solution, improve the capture accuracy of large curvature interface, and improve the overall accuracy and robustness of interface capture while ensuring computational efficiency.

[0015] To achieve the above objectives, this invention provides a high-order interface capture method based on hybrid precision coupling level set and fluid volume, comprising the following steps:

[0016] S1. Divide the computational domain into a grid containing several control volume elements;

[0017] S2. Initialize the VOF field, LS field and velocity field of each control unit;

[0018] S3. Identify the interface unit, and on the interface unit, based on the LS value of the center of each control body unit on the preset template, reconstruct the THINC interface function using a higher-order polynomial.

[0019] S4. Within the interface unit, based on the constraint relationship that the THINC interface function and the VOF field satisfy the volume fraction conservation, a nonlinear equation about the interface position is established.

[0020] S5. The nonlinear equation is solved using a mixed-precision iterative algorithm to obtain the convergent interface position. The mixed-precision iterative algorithm adaptively switches to multi-precision operation when double-precision operation fails to converge.

[0021] S6. Based on the interface position obtained in S5, determine the THINC interface function at the current time, calculate the interface flux based on the THINC interface function at the current time, and obtain the VOF field at the next time by solving the VOF transport equation.

[0022] S7. Based on the THINC interface function and velocity field at the current moment, predict the interface geometry at the next moment;

[0023] S8. Based on the predicted interface geometry of the next moment, the LS field of the next moment is directly reconstructed by a geometric method within a preset narrow band region surrounding the interface. The geometric method includes generating a point cloud near the interface and searching for the nearest point of the center point of the control volume unit on the continuous interface represented by the higher-order polynomial.

[0024] S9. Repeat S3-S8 to capture the evolution of the interface.

[0025] Optionally, the control unit can be a polygonal or polyhedral unit of any shape.

[0026] Optionally, in S3, the higher-order polynomial is a quadratic or cubic polynomial.

[0027] Optionally, when using a cubic polynomial to refactor the interface, the implementation methods include:

[0028] First, solve for the coefficients of the lower-order terms;

[0029] Then, calculate the coefficients of the cubic term based on the coefficients of the lower-order terms;

[0030] Finally, the coefficients of the lower-order terms are corrected to achieve a higher overall polynomial reconstruction.

[0031] Optionally, in S5, the mixed-precision iterative algorithm includes:

[0032] The nonlinear equations are solved using a second-order homotopy iterative formula in double precision.

[0033] If the nonlinear equation fails to converge under the double precision method, then the nonlinear equation is transformed by variable translation and solved iteratively again under the double precision method.

[0034] If convergence is not achieved after variable shifting, the system adaptively switches to multi-precision arithmetic; wherein the effective number of bits in the multi-precision arithmetic is dynamically determined based on the absolute value of the polynomial value at the integration point in the nonlinear equation.

[0035] Optionally, in the multi-precision calculation, an iterative formula for the dynamic convergence control parameters is used for solving, and the dynamic convergence control parameters are automatically updated in each iteration by solving auxiliary equations.

[0036] Optionally, in S8, the LS field at the next time step is directly reconstructed using geometric methods, specifically including:

[0037] S81. Calculate the LS value at the vertex based on the LS field, obtain the intersection point between the interface and the unit edge according to the LS value at the vertex, and then sample along the intersection point to generate an initial point cloud.

[0038] S82. Project the initial point cloud onto the continuous interface represented by the higher-order polynomial to obtain the interface point cloud located on the continuous interface.

[0039] S83, where the center point of the control unit is, search for the discrete nearest point from the interface point cloud;

[0040] S84. Using the discrete nearest point as the initial value, iteratively search on the continuous interface to obtain the true nearest point of the center point, and calculate its distance as the LS value at the center point.

[0041] Optionally, in S84, the search for the true nearest point on the continuous interface is achieved by establishing and solving a constrained optimization problem constituted by the Lagrange multiplier method.

[0042] Optionally, in S8, the preset narrowband region includes an interface unit, a first-layer neighbor unit that shares vertices with the interface unit, and a second-layer neighbor unit that shares vertices with the first-layer neighbor unit.

[0043] 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 high-order interface capture method for coupling level sets and fluid volumes based on mixed precision as described above.

[0044] The high-order interface capture method based on the coupling level set and fluid volume with mixing accuracy provided by this invention has at least one of the following beneficial effects:

[0045] (1) By introducing a mixed precision strategy, the problem of non-convergence in interface position calculation caused by rounding errors in double-precision arithmetic is solved. Unlike most mixed precision methods that combine single-precision and double-precision arithmetic, this invention combines double-precision and multi-precision arithmetic. Multi-precision arithmetic is specifically used for the derivation and iterative root-finding process of the interface position equation, while the rest of this invention still uses double-precision arithmetic. In this way, rounding errors are effectively eliminated without sacrificing overall computational efficiency;

[0046] (2) A hybrid precision strategy is introduced to obtain the convergent interface location, which provides a reliable basis for the geometric reconstruction of horizontal set field (LS);

[0047] (3) The phase interface is directly represented by quadratic and cubic polynomials, which improves the geometric fidelity of the high curvature interface and is especially suitable for interfaces with sharp corners or high curvature.

[0048] (4) Based on the curved surface interface, the smooth horizontal field is reconstructed through geometric methods, thereby providing a more accurate representation for the THINC formula;

[0049] (5) For the first time, a high-order THINC interface is combined with geometric minimum distance construction to realize a high-order CLSVOF coupling framework that does not require re-initialization and is applicable to unstructured meshes. Attached Figure Description

[0050] 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:

[0051] Figure 1 A flowchart illustrating a high-order interface capture method for coupled level sets and fluid volume based on hybrid precision, provided in an embodiment of the present invention;

[0052] Figure 2 This is a schematic diagram of a template for calculating polynomial coefficients provided in an embodiment of the present invention;

[0053] Figure 3 This is a schematic diagram of the projection point onto a higher-order THINC surface provided in an embodiment of the present invention;

[0054] Figure 4 This is a schematic diagram of finding the nearest point in S8 according to an embodiment of the present invention;

[0055] Figure 5 This is a schematic diagram of the units included in a narrowband according to an embodiment of the present invention;

[0056] Figure 6 This is a graph showing the relationship between the final iterative residual and the precision under double precision and multi-precision conditions, provided as an embodiment of the present invention.

[0057] Figure 7 A comparison diagram of the reconstruction interface and the parsing interface under the double-precision and multi-precision frameworks provided in an embodiment of the present invention;

[0058] Figure 8 A trend diagram of interface position error as a function of steepness factor β is provided for an embodiment of the present invention;

[0059] Figure 9 A graph showing the relationship between the number of iterations and the phase volume fraction under different iterative algorithms provided in an embodiment of the present invention;

[0060] Figure 10 A comparison diagram of the reconstruction interface of the square hole when using quadratic polynomial reconstruction and cubic polynomial reconstruction, provided for an embodiment of the present invention;

[0061] Figure 11 A comparison diagram of the interface of the notched disk calculated using the present invention and the OpenFOAM geometric VOF algorithm provided in an embodiment of the present invention;

[0062] Figure 12 This is a diagram of the star-shaped rotating flow interface after five cycles, calculated using high-order surface reconstruction, provided in an embodiment of the present invention.

[0063] Figure 13 A comparison chart of numerical error and computational efficiency between the algorithm of the present invention and the THINC-scaling algorithm provided in an embodiment of the present invention;

[0064] Figure 14 A comparison chart showing the reconstruction accuracy of the algorithm of the present invention and two existing advanced algorithms for the same analytical circular interface, provided as an embodiment of the present invention.

[0065] Figure 15 This is a comparison diagram of quadratic and cubic polynomial reconstructions in a frontal flow according to an embodiment of the present invention.

[0066] Figure 16 A comparison diagram of the algorithm interface of this invention with HyMOFLS, MOF and CLSVOF in two-dimensional shear deformation flow provided in an embodiment of this invention;

[0067] Figure 17 The maximum iterative residual and the percentage of cells using mixed precision on three resolution tetrahedral meshes under mixed precision are shown in an embodiment of the present invention.

[0068] Figure 18 The three-dimensional shear deformation flow interface evolution diagram calculated by the algorithm of the present invention is provided in an embodiment of the present invention. Detailed Implementation

[0069] 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.

[0070] 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.

[0071] 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.

[0072] 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.

[0073] Please refer to Figure 1 To address the shortcomings of existing CLSVOF methods and related interface capture methods in terms of numerical stability, interface location accuracy, and high-order geometric representation, this invention provides a high-order interface capture method (named the CLSTOH method) based on hybrid accuracy coupled level sets and fluid volumes. This method improves the numerical convergence of interface location determination, enhances the capture accuracy of interfaces with large curvatures, and improves the overall accuracy and robustness of interface capture while maintaining computational efficiency. The high-order interface capture method includes the following steps:

[0074] S1. Divide the computational domain into a grid containing several control volume elements;

[0075] S2. Initialize the VOF field, LS field and velocity field of each control unit;

[0076] S3. Identify the interface unit, and on the interface unit, based on the LS value of the center of each control body unit on the preset template, reconstruct the THINC interface function using a higher-order polynomial.

[0077] S4. Within the interface unit, based on the constraint relationship that the THINC interface function and the VOF field satisfy the volume fraction conservation, a nonlinear equation about the interface position is established.

[0078] S5. The nonlinear equation is solved using a mixed-precision iterative algorithm to obtain the convergent interface position. The mixed-precision iterative algorithm adaptively switches to multi-precision operation when double-precision operation fails to converge.

[0079] S6. Based on the interface position obtained in S5, determine the THINC interface function at the current time, calculate the interface flux based on the THINC interface function at the current time, and obtain the VOF field at the next time by solving the VOF transport equation.

[0080] S7. Based on the THINC interface function and velocity field at the current moment, predict the interface geometry at the next moment;

[0081] S8. Based on the predicted interface geometry of the next moment, the LS field of the next moment is directly reconstructed by a geometric method within a preset narrow band region surrounding the interface. The geometric method includes generating a point cloud near the interface and searching for the nearest point of the center point of the control volume unit on the continuous interface represented by the higher-order polynomial.

[0082] S9. Repeat S3-S8 to capture the evolution of the interface.

[0083] First, execute S1 to divide the computational domain into a grid containing several control volume elements.

[0084] In this embodiment, the computational region Ω is divided into a series of non-overlapping control volume units of arbitrary shape. (i=1, 2, ..., The volume of its unit is Each control unit contains K vertices. (k=1, 2, ...,K) and J faces (j=1, 2, ..., J), where the unit outward normal vector of each face is... The area is Optionally, the control volume unit can be a polygonal or polyhedral unit of any shape.

[0085] Then, S2 is executed to initialize the VOF field, LS field, and velocity field of each control volume element. In this embodiment, after the mesh generation is completed, the VOF volume average value and the LS point value of the volume center are initialized according to the analytical expression of the interface at the initial time. The VOF initialization tool can be either the setAlphaFields module in OpenFOAM v2206 or the vofTools tool to perform second-order convergent initialization assignment. The initialized VOF field and LS field are denoted as follows: and Initialize the face-center velocity based on the velocity expression at the initial moment. and the speed of the body and mind The volume average velocity is assigned through numerical integration, and the coordinates and weights of the numerical integration points are given by algorithms in existing techniques [FD Witherden, PE Vincent, On the identification of symmetric quadrature rules for finite element methods, Computers & Mathematics with Applications 69 (2015) 1232–1241]. For general polygonal elements, triangulation is used to triangulate the polygon into triangles during numerical integration.

[0086] Next, S3 is executed to identify the interface units, and on the interface units, based on the LS values ​​of the centers of each control unit on the preset template, the THINC interface function is reconstructed using a higher-order polynomial. For those satisfying... The element meeting the condition = 1e-8 is defined as an interface element. A multidimensional hyperbolic tangent function H(x,t) is used on the interface element to approximate the discontinuous phase indicator function piecewise. H can be expressed as:

[0087] (1)

[0088] in, Representation unit The characteristic length is calculated using the following formula:

[0089] (2)

[0090] In the formula, It is the kth vertex Shape function at the location, As vertex The unit normal gradient at the interface is calculated using the following formula:

[0091] (3)

[0092] In actual calculations, the interface normal gradient is given by the unit gradient of the least squares (LS) algorithm. Given the LS values ​​of all units sharing a vertex, the LS gradient at that vertex can be obtained using the least squares method. The interface is reconstructed in the THINC function. It can be represented by a p-th order polynomial, i.e.

[0093] (4)

[0094] In the formula The coefficients are polynomial coefficients, which can be as follows: Figure 2 The calculations were performed on preset templates of different grid cells.

[0095] The known LS value at the center of each control volume element is given on this preset template. The required solution is for the first and second derivatives of the LS values. For convenience, let's denote the set of LS values ​​on the preset template as... .

[0096] In S3, the higher-order polynomial is a quadratic or cubic polynomial.

[0097] Optionally, a second-order polynomial can be used to reconstruct the interface in two dimensions, i.e., p=2.

[0098] (5)

[0099] in, Represents the body-center coordinates of the target element. When the number of templates is sufficient, the polynomial coefficients in the least squares sense... The following conditions must be met.

[0100] (6)

[0101] This therefore generates a series of algebraic equations, namely

[0102] (7)

[0103] in The above system of algebraic equations can be converted into matrix form. ,in

[0104] (8)

[0105] X is the solution vector, and the unknown coefficients can be solved using the Singular Orthogonal Decomposition (SVD) method. get.

[0106] Optionally, the interface can also be represented using a third-order polynomial, i.e.

[0107] (9)

[0108] Despite employing a higher-order reconstructed polynomial, the template for solving the coefficients remains limited to... Figure 2 The process is performed on the compact template shown, thus avoiding parallel and boundary handling issues caused by template expansion. In this embodiment, the cubic polynomial reconstruction mainly consists of the following three steps:

[0109] 1) Determine the coefficients of the lower-order terms. The coefficients of the lower-order terms can be obtained through the steps above. For convenience, the second-order polynomial can be expanded into the following form:

[0110] (10)

[0111] in, , The relationship between polynomial coefficients and the LS derivative can be written as:

[0112] (11)

[0113] 2) Calculation of coefficients for higher-order terms. Adding the coefficients of the cubic terms to the expression for cubic reconstruction yields:

[0114] (12)

[0115] This yields an expression for the cubic coefficient and the coefficients of the lower-order terms, namely

[0116] (13)

[0117] in, , , The derivative can be obtained on the template of the face adjacency through a least-squares linear approximation, expressed by the formula:

[0118] (14)

[0119] in, To be compatible with the target unit The body center position vector of the j-th adjacent control volume unit. Substituting the derivatives of the coefficients in the above formula into the expression for cubic coefficients yields the coefficients of higher-order terms.

[0120] 3) Correcting the coefficients of lower-order terms to ensure numerical accuracy. To achieve overall fourth-order accuracy in the interface function reconstruction, this invention decomposes the cubic polynomial representation into the sum of the third-order principal term and lower-order correction terms. While retaining the coefficients of the third-order polynomial determined in the second stage, the lower-order coefficients are corrected to satisfy predetermined constraints. Specifically, based on the target unit... Within the template region S centered on the template, for any adjacent element... Construct the following relation:

[0121] (15)

[0122] The third-order principal term is represented as follows:

[0123] (16)

[0124] in , The coefficients of the low-order polynomial to be corrected. This represents the volume center coordinates of the m-th control volume element within the template. Furthermore, the above relationship can be expressed in matrix form. ,in

[0125]

[0126] Due to the matrix Since the inverse matrix has been pre-calculated in the previous solution stage, there is no need to reconstruct it in this stage, thus reducing computational overhead and improving algorithm efficiency. Through the above decomposition and correction steps, while maintaining the stability of the third-order terms, the overall fourth-order precision polynomial reconstruction can be achieved by adjusting the lower-order terms, thereby improving the accuracy of the interface geometric representation.

[0127] Next, S4 is executed. Within the interface unit, based on the constraint relationship that the THINC interface function and the VOF field satisfy the volume fraction conservation, a nonlinear equation about the interface position is established. To ensure mass conservation during the interface reconstruction process, this invention provides a method for solving the interface position based on the volume integral constraint of the THINC function, specifically including the following steps:

[0128] S41. Establish the constraint relationship of volume fraction conservation, that is, within the target control volume element, the volume average value of the THINC function integral is equal to the volume fraction within the element, expressed by the formula:

[0129] (17)

[0130] For convenience, remember

[0131] (18)

[0132] in, For interface location, The polynomial obtained from S3 (excluding constant terms).

[0133] S42. Perform numerical integration on the constraint relations.

[0134] (19)

[0135] in, and Let be the coordinates and weights of the numerical integration points.

[0136] (20)

[0137] The equation is further simplified to

[0138] (twenty one)

[0139] The equation is about the interface position. This is a nonlinear algebraic equation that requires numerical iteration to solve. A significant advantage of this invention is that it uses mixed-precision methods to solve this nonlinear algebraic equation to ensure interface position convergence, thereby improving the accuracy of interface reconstruction.

[0140] Then, S5 is executed, and the nonlinear equation is solved using a mixed-precision iterative algorithm to obtain the convergent interface position. The mixed-precision iterative algorithm adaptively switches to multi-precision operation when double-precision operation fails to converge.

[0141] In this embodiment, to solve for the unknowns corresponding to the aforementioned interface position equation, this invention constructs a hybrid precision iterative solution method. This method adaptively switches between double-precision and multi-precision arithmetic based on the computational convergence, thereby ensuring the numerical stability and accuracy of the interface position solution. Its specific implementation includes the following steps:

[0142] S51. First, the interface position equation is defined as the following nonlinear function.

[0143] (twenty two)

[0144] By solving The solution for obtaining the interface position.

[0145] By default, a double-precision algorithm is used to perform iterative calculations. Initial values ​​are determined based on the sign relation. Right now

[0146] (twenty three)

[0147] in It is a very small number, which can be set to 1e-15 under double precision conditions. Given an initial value, the unknown is updated using the second-order homotopy iterative formula, i.e.

[0148] (twenty four)

[0149] In this invention, h = -1 is directly taken in the double-precision case. Iteration starts from k = 0 and continues until the set maximum number of iterations is reached. =30. Furthermore, when

[0150] (25)

[0151] as well as

[0152] (26)

[0153] Where tol is 1e-12, the iteration terminates. This means that the iteration returns when the solution no longer changes within double precision or when the absolute value of the solution is less than the given maximum iteration residual. The value of .

[0154] S52, When the steepness factor in the THINC function When the value is large or the mesh quality deteriorates, the above double-precision calculation may encounter convergence difficulties. Therefore, this invention performs a variable translation transformation on the original objective function, resulting in...

[0155] (27)

[0156] in,

[0157] (28)

[0158] In the above formula Defined as

[0159] (29)

[0160] parameter Based on the adaptive selection of the extreme value of the integration point, the second-order homotopy iterative scheme described above is used again to solve the problem. Under normal circumstances, most interface elements can converge within the double precision range.

[0161] S53, When the steepness factor in the THINC function When the mesh size is too large or the mesh quality deteriorates sharply, such as in dynamic mesh computation, a small number of cells may fail to converge due to rounding errors or insufficient numerical precision. In such cases, this invention automatically switches to a multi-precision algorithm.

[0162] The significant figures required for multi-precision are adaptively selected according to the following formula.

[0163] (30)

[0164] in, This represents the maximum absolute value of the polynomial at the integration point. As a safety threshold, This is the round-up operator. This adaptive precision strategy avoids the computational overhead of overusing high precision to some extent, maximizing computational efficiency while minimizing the impact of rounding errors.

[0165] Preferably, in the multi-precision calculation, an iterative formula for the dynamic convergence control parameters is used for solving, and the dynamic convergence control parameters are automatically updated in each iteration by solving auxiliary equations.

[0166] In multi-precision applications, an iterative formula for dynamically convergent control parameters is used, namely...

[0167] (31)

[0168] In the formula

[0169] (32)

[0170] After each iteration, the convergence control parameter h is updated by solving the following auxiliary equation, i.e.

[0171] (33)

[0172] in

[0173] (34)

[0174] as well as This dynamic homotopy mechanism requires no manual parameter adjustment and can automatically select a favorable convergence path, improving numerical stability. Through the aforementioned mixed-precision iterative mechanism, this invention can achieve:

[0175] 1. Ensure the interface position equation is stable and convergent;

[0176] 2. Eliminate non-convergence issues caused by double-precision rounding errors;

[0177] 3. Enable multi-precision calculation only in necessary units to avoid increasing overall computational overhead;

[0178] 4. Automatically adjust the convergence path to improve algorithm robustness;

[0179] 5. Provides a reliable basis for interface position reconstruction of high-order curved surfaces.

[0180] Then, S6 is executed. Based on the interface position obtained in S5, the THINC interface function at the current time is determined, and the interface flux is calculated based on the THINC interface function at the current time. The VOF field at the next time is obtained by solving the VOF transport equation.

[0181] Specifically, firstly, in the control body and interval Integrating the transport equation of the phase indicator function yields its semi-discrete scheme.

[0182] (35)

[0183] in, represent The boundary surface, dS represents the direction. The differential area normal vector. According to the definition of the mesh, we have = ,therefore It can be further expressed as

[0184] (36)

[0185] In the formula The surface average value representing the volume fraction. Let represent the velocity flux on the j-th surface. Numerical integration of the above equation using Simpson's integral rule yields...

[0186] (37)

[0187] In the formula, Let the integral weights be the surface average of the volume fraction and the surface average of the velocity fraction at times n+1 / 2 and n+1, respectively, which can be expressed as:

[0188] (38)

[0189] in Let be the volume average velocity at time n.

[0190] Next, S7 is executed. Based on the THINC interface function and velocity field at the current moment, the interface geometry at the next moment is predicted. After completing the interface reconstruction at time n, based on the assumption that the interface shape does not change within one time step, the interface at time n is advanced to obtain the interface at time n+1, expressed by the formula:

[0191] (39)

[0192] based on The LS field at time n+1 is reconstructed for solving the polynomial coefficients in S3. Compared with the traditional THINC method which uses VOF to calculate coefficients, the use of a smoother LS field improves the accuracy of polynomial coefficient calculation.

[0193] Finally, S8 is executed, based on the predicted interface geometry of the next moment, and the LS field of the next moment is directly reconstructed in a preset narrow band region around the interface by a geometric method. The geometric method includes generating a point cloud near the interface and searching for the nearest point of the center point of the control volume unit on the continuous interface represented by the higher-order polynomial.

[0194] In S8, the LS field at the next time step is directly reconstructed using geometric methods, specifically including:

[0195] S81. Calculate the LS value at the vertex based on the LS field, obtain the intersection point between the interface and the unit edge according to the LS value at the vertex, and then sample along the intersection point to generate an initial point cloud.

[0196] S82. Project the initial point cloud onto the continuous interface represented by the higher-order polynomial to obtain the interface point cloud located on the continuous interface.

[0197] S83, where the center point of the control unit is, search for the discrete nearest point from the interface point cloud;

[0198] S84. Using the discrete nearest point as the initial value, iteratively search on the continuous interface to obtain the true nearest point of the center point, and calculate its distance as the LS value at the center point.

[0199] First, S81 is executed, which obtains the LS value at the vertex based on least squares linear interpolation.

[0200] (40)

[0201] in, This represents the target unit Sharing the l-th unit of the k-th vertex, Indicates from vertex Pointing unit The position vector. Once the LS value at the vertex is obtained, each edge of the cell is traversed to obtain the position vector. The coordinates of the intersection points of the contour lines and the boundary. The existence of an intersection point is based on whether the LS values ​​of the two vertices of each edge have opposite signs, such as... Figure 3 As shown.

[0202] Within this unit, vertices i1 and i4 have opposite signs, and vertices i2 and i3 have opposite signs. Therefore, there are two intersection points. The position of the intersection point is calculated using the following formula.

[0203] (41)

[0204] In the formula - Represents the vertex coordinates. Connect. and Scattering points evenly between two points, and taking four sample points for calculation in the two-dimensional case of this invention, the sample points are obtained. - .

[0205] Next, step S82 is executed. After calculating the sample points within each cell, the sample points are projected along the interface normal onto the continuous interface represented by the higher-order polynomial of THINC, resulting in an interface point cloud located on the continuous interface. This process can be described by the formula:

[0206] (42)

[0207] For each point Using the geometric information of the THINC interface, Newton's iteration is employed to obtain the points on the interface. This iterative process only requires a few iterations to converge. For all interface units, repeat steps S81 and S82 above to finally obtain the interface point cloud of the entire interface. ,in , = It represents the collection of all segmented interfaces.

[0208] Then execute S83 to generate point clouds on the interface. Then, given a query point x, it is necessary to start from... Find the point closest to x. Expressed by the formula:

[0209] (43)

[0210] Optionally, a KD-tree can be used to process the point cloud set of the interface. Establish a query structure to find the nearest point of x on the discrete interface. In this invention, the query point is selected as the center point of the unit.

[0211] After executing S84 and obtaining the nearest point in the discrete sense, it is necessary to find the true nearest point of x on the THINC piecewise continuous interface. In this invention, the Lagrange multiplier method is used to establish constraint relationships, which satisfy…

[0212] (44)

[0213] It indicates that in Under the constraints, find the shortest distance between y and x. Generally, let x be the center point of the target cell, and y be the nearest point to x on the interface. This optimization problem can be solved using the Newton-Raphson iteration method, letting... Then the iterative formula can be written as

[0214] (45)

[0215] In the formula, for The gradient of H is the Hessian matrix, which can be calculated from the gradient of the polynomial interface in the THINC function and the Hessian matrix.

[0216] (46)

[0217] The initial value for the iteration is given by the following formula.

[0218] (47)

[0219] in, The initial Lagrange multiplier, the initial point Let x be the closest point in the discrete sense of S83. After several iterations converge, let... This yields the nearest point on the continuous x-segment interface. To more intuitively illustrate the specific process of nearest point search in S82 to S84, as follows... Figure 4 As shown.

[0220] Specifically: Figure 4 (a) Demonstrates the generation of interface point clouds using the projection algorithm in step S2. The process Figure 4 (b) demonstrates the process by which S3 finds the nearest point in a discrete sense, regardless of the description. Figure 4 (c) then gives the following: The process of finding the nearest point in the continuous sense, starting from the nearest point in the discrete sense.

[0221] S81-S84 can be applied to all meshes to generate a global LS field. However, in this invention, to maximize computational efficiency while maintaining accuracy, a narrowband local reconstruction mechanism is employed. The elements contained in the narrowband are as follows: Figure 5 As shown, it includes interface units (denoted as...) The first-level neighboring unit that shares vertices with the interface unit. and with Second-level neighbors of shared vertex units .

[0222] For example, for cells within a narrow band Its center point The LS value at that location is calculated using the following formula:

[0223] (48)

[0224] in

[0225] (49)

[0226] Calculated Used for calculating the polynomial coefficients at time n+1.

[0227] The high-order interface capture method based on the coupling level set and fluid volume with hybrid accuracy proposed in this invention has the following significant features:

[0228] 1) A hybrid precision strategy and a smooth LS field are used to improve the accuracy of surface interface reconstruction within the THINC framework, thereby achieving high-fidelity geometric reconstruction of the LS field.

[0229] 2) The coupling technique between algebraic flux calculation and geometric LS reconstruction based on high-order surface interfaces distinguishes this invention from traditional CLSVOF and other THINC-LS methods.

[0230] At the same time, these features bring the following major advantages:

[0231] 1) The mixed precision strategy ensures robust convergence of the interface position. Even when the mesh is highly distorted, the mixed precision iteration strategy in this invention remains effective.

[0232] 2) The converged interface location provides a guarantee for the high-precision geometric reconstruction of the LS field, thus improving the accuracy of the LS geometric reconstruction;

[0233] 3) Relying on the smooth LS function, the interface expression in the THINC function can use a higher-order polynomial, thereby improving the accuracy of capturing interfaces with large curvature.

[0234] 4) The smooth LS field enables more accurate polynomial coefficient calculation, thereby improving the accuracy of interface capture.

[0235] To further demonstrate the technical effects of the present invention, the following examples are used to test the algorithm proposed in this invention.

[0236] Test 1: Reconstructing a circle: Given a circle with radius r and center r... The circle is reconstructed into a surface within the THINC framework, focusing on the convergence of the interface position, the advantages of mixed accuracy, and the advantages of homotopy iteration with dynamic control of convergence factor compared to traditional Newton iteration.

[0237] 1) Advantages of hybrid precision. Given an interface element, its interface equation is expressed as:

[0238] (50)

[0239] in, The volume fraction of the unit is Using 60 integration points, the results were evaluated in both double-precision and multi-precision modes. The value of . Figure 6 The evaluation methods using double precision and multiprecision are presented at different levels of accuracy. The final interface position obtained by solving absolute value of a function The distribution. From Figure 6 It is evident from the image that the rounding error caused... The computational accuracy decreases, and the interface position equation cannot be solved numerically in double precision. Conversely, when and When multi-precision arithmetic is used, rounding errors in the calculation of the hyperbolic tangent function are effectively suppressed. This strongly demonstrates that the mixed-precision strategy in this invention guarantees the numerical solvability of the interface position equation, thereby ensuring the convergence of the interface position.

[0240] 2) Advantages of high precision on twisted meshes. Given a circle with a radius and center at (0.525, 0.464), its phase indicator function is expressed as...

[0241] (51)

[0242] Reconstruct the circular interface on a twisted mesh, and output the polynomial interface for each interface unit. =0 and compare it with the parsing interface. Figure 7 The reconstruction interfaces under double precision and the multi-precision framework of this invention are compared, with the red line representing the reconstruction interface and the black line representing the parsing interface, and steepness factor. .

[0243] Given an expression for a circle, interface position error can be quantitatively calculated. For a given circle, its analytical expression is generally expressed as:

[0244] (52)

[0245] Transform it to a local coordinate system, let , ,as well as

[0246] The converted formula is:

[0247] (53)

[0248] The polynomial coefficients derived in the formula are exact coefficients. The interface position is determined by the analysis. Based on this, a statistical formula for the interface position error can be quantitatively given, namely...

[0249] (54)

[0250] Figure 8 The relationship between interface position error and β after using double-precision and multi-precision arithmetic is presented. When using multi-precision arithmetic, the interface position error decreases as β increases, which is consistent with theoretical expectations, i.e., for larger β values, the reconstructed interface will converge to the analytical interface. Conversely, when only double-precision arithmetic is used, the error increases as β increases, deviating from the expected convergence trend.

[0251] 3) Advantages of the Homotopy Iterative Algorithm with Dynamically Adjusted Convergence Control Factor. This invention employs a dynamic homotopy iterative formula to solve the interface position equation. To demonstrate the advantages of the dynamic homotopy iterative formula in this invention, a comparison of the number of iterations under different iterative algorithms is presented, such as... Figure 9 As shown.

[0252] Compared to the traditional Newton's iteration algorithm and the second-order homotopy iteration (HAM) formula with h=-1, the dynamic h homotopy iteration algorithm used in this invention significantly reduces the number of iterations, which can improve computational efficiency in mixed-precision mechanisms. Specifically, compared to the traditional Newton method, the dynamic h homotopy iteration algorithm can reduce the number of iterations to less than half. Please refer to Table 1 below for a comparison of the average number of iterations for different algorithms.

[0253] Table 1 Comparison of average number of iterations for different iterative algorithms

[0254]

[0255] Test 2: 1) This invention uses cubic polynomials for interface reconstruction, which shows a significant advantage over traditional linear and quadratic polynomial reconstructions in capturing interfaces with large curvatures, thus improving interface capture accuracy. This invention uses classic VOF examples for verification, including the translational motion of a square with a hole along its diagonal, the periodic rotation of a notched disk, and the periodic rotation of a star-shaped object.

[0256] Figures 10-12 A comparison between reconstructed and analyzed interfaces under different interface shapes is presented. The cubic polynomial method used in this invention maintains higher geometric fidelity when capturing sharp corners and interfaces with high curvature, reducing interface capture errors by approximately 50%. This significantly demonstrates the advantages of using high-order surface reconstruction in this invention. Furthermore, compared to the mainstream high-precision geometric VOF class in the current open-source software OpenFOAM, high-order surface reconstruction, compared to linear reconstruction, ensures the symmetry of the interface, thus improving the overall accuracy of interface capture.

[0257] 2) The high-order surface reconstruction and interface capture method of the present invention has higher computational efficiency and accuracy, referring to... Figure 13 Compared to the THINC-scaling algorithm, the present invention uses a quadratic polynomial to reduce the computation time to 1 / 5 to 1 / 12, and uses a cubic polynomial to improve the computation efficiency by 300% to 500%, while the numerical error is only half that of THINC-scaling. Table 2 below shows the quantitative numerical error.

[0258] Table 2. Errors between the proposed algorithm and other algorithms in rotating flow on a notched disk.

[0259]

[0260] The references mentioned in the table are as follows:

[0261] [1] D. Chen, B. Xie, F. Xiao, Revisit to the thinc / qq scheme: Recent progress to improve accuracy and robustness, International Journal for Numerical Methods in Fluids 94 (2022) 719–755

[0262] [2] Y. Chen, W. Lu, D. Sun, B. Yu, J. Gong, W. Zhang, W. Tao, Ahorizontal refined piecewise curve interface reconstruction (hopcir) algorithm for reconstructing the vapor-liquid interface, International Journal of Multiphase Flow 178 (2024) 104905.

[0263] Test 3: In this invention, the geometrically reconstructed LS field is used to calculate the coefficients of the polynomial in the THINC function. Compared with the traditional VOF algorithm, the smooth LS field improves the accuracy of coefficient calculation, thereby ensuring the accuracy of higher-order geometric information of the interface, such as normal vectors and curvature. The translation of a two-dimensional circle along its diagonal is used as an example to illustrate the advantages of using LS to calculate the polynomial coefficients of the THINC function in this invention.

[0264] from Figure 14 It is evident that, on two different grids, after one cycle, the numerical interface of the THINC method based on VOF for calculating coefficients deviates significantly from the analytical interface, while the numerical interface of the algorithm in this invention is almost indistinguishable from the analytical interface to the naked eye after one cycle. This further demonstrates the significant advantage of the VOF and LS coupling strategy adopted in this invention in maintaining the interface geometry.

[0265] Test 4: Two-dimensional shear deformation flow and frontal flow. This test demonstrates the algorithm's ability to analyze interfaces.

[0266] from Figure 15It can be seen that the cubic polynomial reconstruction in this invention can resolve extremely thin interfaces. Compared with the quadratic polynomial reconstruction in the traditional THINC algorithm, the cubic polynomial reconstruction provides a better match between the interface and the analytical interface, further demonstrating the advantages of high-order surface reconstruction. Furthermore, in two-dimensional shear deformation flow, the algorithm in this invention significantly outperforms the traditional CLSVOF algorithm in resolving interfaces on structured meshes, and is comparable to MOF and hybrid algorithms combining MOF and LS. However, the algorithm in this invention is applicable to arbitrary unstructured meshes, exhibiting broader mesh universality, and is simpler, making it easier for researchers in this field to use. Please refer to [reference needed]. Figure 16 , Figure 16 For two-dimensional shear deformation flow, the algorithm interface of this invention ( Figure 16 d- Figure 16 f) and HyMOFLS ( Figure 16 a), MOF ( Figure 16 b) and CLSVOF ( Figure 16 c) Comparison.

[0267] Test 5: Three-dimensional shear deformation flow, designed to demonstrate the effectiveness of the hybrid precision strategy of this invention in capturing large deformation interfaces. Three-dimensional shear deformation flow is a classic convection example. The initial interface is a sphere, which reaches its maximum deformation in half a cycle under the action of the shear velocity field. Then the velocity field reverses, and ideally, it can recover to its original shape.

[0268] Figure 17 The performance of the proposed mixed-precision iterative strategy in solving for interface locations is demonstrated on low-quality tetrahedral meshes. It can be seen that all interface elements converge on all meshes due to the mixed-precision strategy; however, the mixed-precision elements account for less than 1% of all interface elements. This indicates that the proposed mixed-precision iterative strategy overcomes the influence of rounding errors while maintaining the high efficiency of a full double-precision strategy, making it suitable for large-scale parallel computing. Figure 18 The algorithm demonstrates its excellent performance in capturing large-scale 3D deformation interfaces. Overall, the interface shape at time t=0 and time t=T is almost unchanged, which further demonstrates the advantages of this invention in interface capture.

[0269] Based on the same inventive concept, embodiments of the present invention also propose a readable storage medium storing a computer program thereon, which, when executed, can realize the high-order interface capture method for coupling level sets and fluid volumes based on mixed precision as described above.

[0270] 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.

[0271] 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.

[0272] 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.

[0273] Since the readable storage medium provided by this invention belongs to the same inventive concept as the high-order interface capture method based on the coupling level set and fluid volume of the above-described mixed precision, the readable storage medium provided by this invention has all the advantages of the high-order interface capture method based on the coupling level set and fluid volume of the above-described mixed precision. Therefore, the beneficial effects of the readable storage medium provided by this invention will not be described in detail here.

[0274] In summary, to address the problems of the traditional CLSVOF method and the THINC-LS coupled method, this invention proposes a high-order interface reconstruction and capture method based on hybrid precision for accurately capturing moving interfaces. Firstly, this invention solves the problem of non-convergence in interface position calculation caused by rounding errors in double-precision arithmetic by introducing a hybrid precision strategy. Unlike most hybrid precision methods that combine single-precision and double-precision, this invention combines double-precision and multi-precision arithmetic. Multi-precision arithmetic is specifically used for the derivation and iterative root-finding process of the interface position equation, while the rest of the invention still uses double-precision arithmetic. This effectively eliminates rounding errors without sacrificing overall computational efficiency. Furthermore, this invention implements adaptive precision adjustment to avoid excessive computational overhead caused by unnecessarily high precision. The resulting hybrid precision algorithm guarantees the numerical validity of the interface position equation and obtains converged interface positions. Based on these converged positions, a smoother horizontal field is directly geometrically reconstructed using the nearest-point algorithm on the high-order surface interface reconstructed by THINC. Then, the reconstructed horizontal field is used to calculate the coefficients of the high-order polynomial in subsequent time steps, thereby improving the overall accuracy of interface capture. Numerous numerical experiments on interface advection benchmarks demonstrate that the algorithm proposed in this invention has stronger numerical robustness to distorted meshes, employs quadratic and cubic polynomial representations to improve geometric fidelity, achieves higher geometric fidelity for interfaces with sharp corners and high curvature, and is competitive in accuracy with state-of-the-art VOF and CLSVOF methods.

[0275] 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 high-order interface capture method based on hybrid precision coupling level set and fluid volume, characterized in that, Includes the following steps: S1. Divide the computational domain into a grid containing several control volume elements; S2. Initialize the VOF field, LS field and velocity field of each control unit; S3. Identify the interface unit, and on the interface unit, based on the LS value of the center of each control body unit on the preset template, reconstruct the THINC interface function using a higher-order polynomial. S4. Within the interface unit, based on the constraint relationship that the THINC interface function and the VOF field satisfy the volume fraction conservation, a nonlinear equation about the interface position is established. S5. The nonlinear equation is solved using a mixed-precision iterative algorithm to obtain the convergent interface position. The mixed-precision iterative algorithm adaptively switches to multi-precision operation when double-precision operation fails to converge. S6. Based on the interface position obtained in S5, determine the THINC interface function at the current time, calculate the interface flux based on the THINC interface function at the current time, and obtain the VOF field at the next time by solving the VOF transport equation. S7. Based on the THINC interface function and velocity field at the current moment, predict the interface geometry at the next moment; S8. Based on the predicted interface geometry of the next moment, the LS field of the next moment is directly reconstructed by a geometric method within a preset narrow band region surrounding the interface. The geometric method includes generating a point cloud near the interface and searching for the nearest point of the center point of the control volume unit on the continuous interface represented by the higher-order polynomial. S9. Repeat S3-S8 to capture the evolution of the interface.

2. The high-order interface capture method for coupled level sets and fluid volume based on hybrid precision according to claim 1, characterized in that, The control unit is a polygonal or polyhedral unit of arbitrary shape.

3. The high-order interface capture method for coupled level sets and fluid volume based on hybrid precision according to claim 1, characterized in that, In S3, the higher-order polynomial is a quadratic or cubic polynomial.

4. The high-order interface capture method for coupled level sets and fluid volume based on hybrid precision according to claim 3, characterized in that, When using a cubic polynomial to reconstruct the interface, the implementation methods include: First, solve for the coefficients of the lower-order terms; Then, calculate the coefficients of the cubic term based on the coefficients of the lower-order terms; Finally, the coefficients of the lower-order terms are corrected to achieve a higher overall polynomial reconstruction.

5. The high-order interface capture method for coupled level sets and fluid volume based on hybrid precision according to claim 1, characterized in that, In S5, the mixed-precision iterative algorithm includes: The nonlinear equations are solved using a second-order homotopy iterative formula in double precision. If the nonlinear equation fails to converge under the double precision method, then the nonlinear equation is transformed by variable translation and solved iteratively again under the double precision method. If convergence is not achieved after variable shifting, the system adaptively switches to multi-precision arithmetic; wherein the effective number of bits in the multi-precision arithmetic is dynamically determined based on the absolute value of the polynomial value at the integration point in the nonlinear equation.

6. The high-order interface capture method for coupled level sets and fluid volume based on hybrid precision according to claim 5, characterized in that, In the multi-precision calculation, an iterative formula for the dynamic convergence control parameters is used for solving. The dynamic convergence control parameters are automatically updated in each iteration by solving auxiliary equations.

7. The high-order interface capture method for coupled level sets and fluid volume based on hybrid precision according to claim 1, characterized in that, In S8, the LS field at the next time step is directly reconstructed using geometric methods, specifically including: S81. Calculate the LS value at the vertex based on the LS field, obtain the intersection point between the interface and the unit edge according to the LS value at the vertex, and then sample along the intersection point to generate an initial point cloud. S82. Project the initial point cloud onto the continuous interface represented by the higher-order polynomial to obtain the interface point cloud located on the continuous interface. S83, where the center point of the control unit is, search for the discrete nearest point from the interface point cloud; S84. Using the discrete nearest point as the initial value, iteratively search on the continuous interface to obtain the true nearest point of the center point, and calculate its distance as the LS value at the center point.

8. The high-order interface capture method for coupled level sets and fluid volume based on hybrid precision according to claim 7, characterized in that, In S84, the search for the true nearest point on the continuous interface is achieved by establishing and solving a constrained optimization problem constructed by the Lagrange multiplier method.

9. The high-order interface capture method for coupled level sets and fluid volume based on hybrid precision according to claim 1, characterized in that, In S8, the preset narrowband region includes an interface unit, a first-layer neighbor unit that shares vertices with the interface unit, and a second-layer neighbor unit that shares vertices with the first-layer neighbor unit.

10. A readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed, it can implement the high-order interface capture method based on hybrid precision coupling level set and fluid volume according to any one of claims 1-9.

Images