Overlapping grid VOF interpolation method based on THINC and readable storage medium

Abstract

The present invention relates to a THINC-based overlapping grid VOF interpolation method and a readable storage medium. The method comprises: obtaining the body center coordinates and Gaussian integration point coordinates of each grid unit in an overlapping grid system; obtaining grid assembly information; selecting an interpolation boundary unit and a corresponding interpolation template; selecting a Gaussian integration point for the selected interpolation boundary unit; calculating the distance from the Gaussian integration point to the grid body center of each contributing unit in the corresponding interpolation template to find the closest contributing unit; judging whether the closest contributing unit is an interface unit, and if not, the fluid volume fraction value of the Gaussian integration point is the fluid volume fraction value of the closest contributing unit; if so, substituting the Gaussian integration point into the THINC function of the closest contributing unit to calculate the fluid volume fraction value of the Gaussian integration point; repeating the above steps to calculate the fluid volume fraction values ​​of all interpolation boundary units in the overlapping grid system, thereby realizing overlapping grid interpolation calculation of the fluid volume fraction field.

Description

Technical Field

[0001] The present invention relates to the field of CFD technology, and in particular to a THINC-based overlapping grid VOF interpolation method and a readable storage medium. Background Art

[0002] In the field of Computational Fluid Dynamics (CFD), the Volume of Fluid (VOF) method has become an important tool for simulating immiscible multiphase fluid interfaces. The core of this method is to dynamically capture the fluid interface by solving the convection equation of the fluid volume fraction (VOF). When implementing CFD calculations, it is first necessary to discretize the physical calculation domain into grids. These grids are composed of specific types of units. The fluid volume fraction (VOF) represents the volume ratio of a certain phase fluid in the grid unit. In order to facilitate

[0003] For the description, let the area occupied by fluid A and B be Ω A and Ω B And introduce the fluid indicator function To define it, the expression is as follows:

[0004]

[0005] Wherein, the position vector is represented by x = (x, y, z). i In the equation, the volume fraction of fluid A is φ i is defined as:

[0006]

[0007] Where, |Ω i | is used to represent the volume of the grid unit. The fluid volume fraction φ physically follows the following convection equation:

[0008]

[0009] In the formula, t represents the time variable, and u represents the velocity vector field. Therefore, in the process of using the VOF method to capture the fluid motion interface, the core problem becomes how to solve the above convection equation with high accuracy. Specifically, it is necessary to calculate the velocity vector field based on the current time t n The fluid distribution state φ n and the corresponding velocity field u n , to accurately deduce the next moment t n+1 Fluid distribution φ n+1 ,This process is crucial to ensure the accuracy and reliability of multiphase flow numerical simulation results.

[0010] The solution methods for the VOF convection equation are mainly divided into two types: geometric reconstruction and algebraic. Geometric reconstruction methods, such as the piecewise linear interface reconstruction method (PLIC), have high accuracy, but the algorithm complexity is also relatively high. In contrast, algebraic VOF methods are usually simpler and easy to extend to unstructured grids, but their accuracy may be limited, such as the HRIC (High Resolution Interface Capturing) method. However, the THINC (Tangent of Hyperbola Interface Capturing) method, by introducing the hyperbolic tangent function and fully considering the interface geometry information, achieves accuracy comparable to that of the geometric reconstruction method without explicit geometric reconstruction, while maintaining the simplicity of the algebraic method, thus showing good application potential. In order to represent the interface more accurately, the paper "Xie B, Xiao F. Toward Efficient and Accurate Interface Capturing on Arbitrary Hybrid Unstructured Grids: The THINC Method with Quadratic Surface Representation and Gaussian Quadrature [J] Journal of Computational Physics, 2017" mentioned the use of quadratic polynomials to represent the interface, and proposed the THINC / QQ (THINC method with quadratic surface representation and Gauss quadrature) format for any type of unstructured grid. This algorithm greatly simplifies the numerical integration process in multidimensional THINC reconstruction by using Gaussian integrals. Compared with the geometric VOF method, the THINC / QQ format is competitive in terms of numerical accuracy and geometric fidelity, and the algorithm has the advantages of simple algorithm and high computational efficiency on polyhedral unstructured grids.

[0011] The multidimensional hyperbolic tangent function is used to piecewise approximate the VOF index function, and the function is named THINC function. Its specific expression is:

[0012]

[0013] Where β represents the steepness coefficient. When β is large enough, the THINC function will asymptotically approach the step-like VOF index function. Therefore, by adjusting the steepness of the two-phase flow transition region, the width of the gas-liquid two-phase flow interface can be effectively controlled. The two-phase flow interface can be implicitly expressed as a quadratic polynomial:

[0014]

[0015] Where X = xx ic , Y=yy ic , Z=zz ic Represents the displacement relative to the center of the target grid cell. The coefficient of the approximate quadratic polynomial C ijk It is estimated based on the unit normal vector and curvature of the interface using the least squares method. The only unknown coefficient is d i Represents the position of the approximate surface. This coefficient must satisfy the constraint of fluid mass conservation, namely:

[0016]

[0017] In the formula, the volume integral is approximated by the Gaussian integral method. For the solution of equation (6), the Newton iteration method can be used, which can reach convergence after a finite number of iterations. After the THINC function is constructed for each interface grid cell, the finite volume method can be used to update the fluid volume fraction of the target cell. According to the Gaussian divergence theorem, the update process can be expressed as the control volume Ω i Perform volume integration, and its mathematical expression is:

[0018]

[0019] Where, Γ ij represents the shared boundary surface between the target grid cell i and its adjacent grid cell j, and It means that in Γ ij The velocity flux pointing to the outside of the target grid cell. The specific calculation method is as follows:

[0020]

[0021] Where n ij Represents the unit surface normal vector pointing outward from the target grid cell. The surface average of the fluid volume fraction is obtained by the mesh interface Γ ij THINC reconstruction function for upwind grid cells The specific calculation expression is obtained by performing surface integration:

[0022]

[0023] Where, subscript i up To represent the mesh surface Γ ij The index number of the adjacent upwind grid cell.

[0024] The main steps of the overset mesh method include generating a background mesh, generating a body-fitted mesh, and performing mesh assembly. First, the background mesh is generated, which discretizes the entire computational domain without considering the geometry of the specific rigid bodies within it. Next, the body-fitted mesh is generated, which constructs an independent body-fitted mesh for each rigid body in the computational domain. Each body-fitted mesh can move freely with the rigid body or according to a preset motion. Figures 1a-1c A typical two-dimensional overlapping grid system is presented, which includes a background grid and a body-fitting grid. After completing the above steps, the grid assembly step is carried out. The specific approach is to integrate the body-fitting grid into the background grid in an overlapping and embedded manner. In order to solve the flow equations in the overlapping grid system, it is necessary to achieve effective transmission of flow field information between different grids through grid interpolation. In order to establish the interpolation relationship between these grids, the cells of the background grid and the body-fitting grid need to be subdivided into three types: hole cells 1, interpolation boundary cells 2, and calculation cells 3. Referring to Figure 1, the three types of cells are clearly marked in red, green, and blue, respectively. In the overlapping grid system, in order to accurately divide the types of grid cells, the following three core steps must be performed in sequence in each time step:

[0025] (1) Searching for Hole Cell 1: The process of searching for Hole Cell 1 involves identifying and shielding the mesh cells outside the computational domain or inside the rigid body. These cells are often referred to as “Hole Cell 1” in the professional field. Figure 1a As shown in the figure, the mesh elements where the background mesh overlaps with the 2D cylinder are clearly marked as red hole elements 1. When solving the flow equations in the overlapping mesh system, a masking value is assigned to hole element 1, thereby removing it from the calculation. This operation not only effectively decouples the fluid and solid domains but also significantly improves computational efficiency.

[0026] (2) Searching for interpolation boundary cells 2: In the process of searching for interpolation boundary cells 2, two specific types of grid cells must be accurately identified: one is the grid cells adjacent to hole cells 1 on the background grid, and the other is the grid cells bordering the grid boundary on the body-fitting grid. Figure 1a and Figure 1b The green area in the figure represents the location of interpolation boundary cells 2. For each such interpolation boundary cell 2, there is a corresponding interpolation template. This template is composed of multiple contribution cells, and its main function is to effectively exchange flow field information between different grids.

[0027] (3) Contribution unit search: When searching for contribution units, for a specific interpolation boundary unit 2 on a certain grid, one or more adjacent computational units 3 must be accurately selected from another set of overlapping grids (i.e., contribution grids). These selected computational units 3 (i.e., contribution units) will constitute the interpolation template of the interpolation boundary unit 2, and their main function is to transmit the flow field information in the contribution grid to the interpolation boundary unit 2.

[0028] Figure 2 The figure clearly shows an interpolation boundary element 2 (highlighted by a red border) and its corresponding interpolation template 4 (marked by a black border). In the application of overlapping grid technology, a correspondence between interpolation boundary elements 2 and interpolation templates 4 is established to ensure accurate transfer of flow field information between grids. This process is the core step in achieving effective coupling of multiple sets of different grids during flow equation calculation and flow field updates.

[0029] Chan et al. provide a detailed and in-depth discussion of the overset grid method in their paper "Chan WM, Pandya S A. Advances in Distance-Based Hole Cuts on Overset Grids [C] / / 22nd AIAA Computational Fluid Dynamics Conference. Dallas, TX, 2015." When dealing with static overset grids with fixed body-fitting boundaries, the solver only needs to perform a one-time mesh assembly operation before the flow calculation. However, when dealing with dynamic mesh boundaries, the solver must perform multiple overset mesh assembly operations based on the real-time motion of the body-fitting grid.

[0030] Traditional overlapping grid interpolation methods usually use the distance inverse weighted algorithm. In this algorithm framework, the physical quantity φ on the interpolation boundary grid is A It is based on the physical quantity on the contributing grid cell To obtain, the calculation formula is:

[0031]

[0032] Where, subscript D i represents the i-th contributing unit in the interpolation template, and It is used to describe the distance between the center of the contribution unit grid and the center of the target interpolation boundary unit 2 grid. For details, please refer to Figure 2 . Weight coefficient The calculation is based on the inverse of the distance to the mesh center. The inverse distance weighted algorithm is not only computationally efficient but also amenable to parallel processing. However, it should be noted that its accuracy may be affected by changes in the relative position of the mesh center of the contributing element and the interpolation boundary element 2.

[0033] Chandar et al. proposed an innovative polynomial-based overlapping grid interpolation algorithm in the paper "Chandar DD On Overset Interpolation Strategies and Conservation on Unstructured Grids in Open FOAM[J].Computer Physics Communications, 2019, 239:72-83". This algorithm uses linear or high-order polynomials to determine the interface quadratic polynomial coefficients using the contributing unit through the least squares method. However, it is worth noting that this polynomial-based interpolation algorithm is only applicable to overlapping grid interpolation of continuous variables (such as velocity fields and pressure fields). For step functions such as fluid volume fraction, this algorithm is not applicable. Summary of the Invention

[0034] The purpose of the present invention is to provide a THINC-based overlapping grid VOF interpolation method and a readable storage medium, which can accurately interpolate the VOF field on an unstructured overlapping grid system and achieve accurate transmission of VOF field information.

[0035] To achieve the above object, the present invention provides an overlapping grid VOF interpolation method based on THINC, comprising:

[0036] S1. Obtain the body center coordinates and Gaussian integration point coordinates of each grid cell in the overlapping grid system;

[0037] S2. Obtaining grid assembly information in the overlapping grid system, wherein the grid assembly information includes an interpolation boundary unit index list and an index list of contribution units in the interpolation template corresponding to each interpolation boundary unit;

[0038] S3, reading in the interface quadratic polynomial coefficients, interface positions, and fluid volume fraction values ​​of each contributing unit;

[0039] S4. Select an interpolation boundary unit from the interpolation boundary unit index list for processing, and select a corresponding interpolation template;

[0040] S5. For the selected interpolation boundary unit, select one Gaussian integral point of the interpolation boundary unit for processing;

[0041] S6. Calculate the distance from the Gaussian integral point to the center of each contributing unit grid in the corresponding interpolation template, and find the contributing unit closest to the Gaussian integral point;

[0042] S7. Determine whether the nearest contributing element is an interface element. If the nearest contributing element is not an interface element, select the fluid volume fraction value of the Gaussian integral point as the fluid volume fraction value of the nearest contributing element. If the nearest contributing element is an interface element, substitute the displacement between the Gaussian integral point and the center of the grid of the nearest contributing element into the THINC function of the nearest contributing element to calculate the fluid volume fraction value of the Gaussian integral point.

[0043] S8, repeating S5-S7, calculating the fluid volume fraction values ​​of all Gaussian integral points of the selected interpolation boundary unit, and calculating the fluid volume fraction value of the selected interpolation boundary unit using Gaussian integral method;

[0044] S9. Repeat S4-S8 to calculate the fluid volume fraction values ​​of all interpolation boundary cells in the overlapping grid system, thereby realizing overlapping grid interpolation calculation of the fluid volume fraction field.

[0045] Optionally, in S6, the distance r from the Gaussian integral point to the center of each contributing unit grid in the corresponding interpolation template is calculated. ci The specific formula is as follows:

[0046]

[0047] In the formula, (x gi ,y gi , z gi ) is the coordinate of the Gaussian integral point, (x ci ,y ci , z ci ) is the coordinate of the center of the contributing unit grid.

[0048] Optional, when ∈≤φ i ≤1-∈,∈=10 -8 When , the nearest contributing unit is the interface unit, otherwise it is not; where φ i is the fluid volume fraction value of the nearest contributing unit.

[0049] Optionally, in S7, if the nearest contributing element is the interface element, the specific formula for calculating the fluid volume fraction value of the Gaussian integral point is as follows:

[0050]

[0051] P i (x gi )=C 200 X 2 +C 020 Y 2 +C 002 Z2 +C 110 XY+C 101 XZ+C 011 YZ

[0052] +C 100 X+C 010 Y+C 001 Z

[0053] X=x gi -x ci , Y=y gi -y ci , Z=z gi -z ci

[0054] In the formula, (x gi ,y gi , z gi ) is the coordinate of the Gaussian integral point, (x ci ,y ci , z ci ) is the coordinate of the center of the contributing unit grid; C 200 、C 020 、C 002 、C 110 、C 101 、C 011 、C 100 、C 010 、C 001 are the coefficients of the interface quadratic polynomial.

[0055] Optionally, in S8, the formula for calculating the fluid volume fraction value of the selected interpolation boundary unit using the Gaussian integral method is as follows:

[0056]

[0057] Where, Weights for the fluid volume fraction values ​​at each Gaussian integration point.

[0058] Based on the same technical concept, the present invention further provides a readable storage medium having a computer program stored thereon, which can implement the THINC-based overlapping grid VOF interpolation method described above when the computer program is executed.

[0059] The present invention provides a THINC-based overlapping grid VOF interpolation method and a readable storage medium, which have at least one of the following beneficial effects:

[0060] 1) Traditional interpolation algorithms easily lead to numerical dissipation of VOF values ​​at the interface. The interpolation method provided in this application is based on the THINC function, which accurately records the interface geometry information through quadratic surface polynomials, thereby effectively maintaining the sharpness of the interface.

[0061] 2) Compared with traditional interpolation algorithms, the interpolation method provided by this application demonstrates a stronger interface capture capability. It not only accurately captures the thin film structure of the VOF field, but also effectively restores the VOF field to its original spherical shape. The resulting interface is sharper and the thin film structure of the interface is more refined.

[0062] 3) It is worth mentioning that the interpolation method provided in this application can maintain the thickness of the interface, so that the interface details of the thin film part are more clearly visible compared to the traditional algorithm. In addition, because the THINC function effectively retains the geometric information of the interface during the calculation process, the visual interface on the overlapping grid and the background grid shows a high degree of consistency. In contrast, the traditional algorithm has obvious numerical dissipation when processing thin films, which makes the VOF field unable to accurately restore to its initial spherical shape. In summary, the interpolation method provided in this application has an excellent ability to capture three-dimensional complex moving interfaces on overlapping grid systems and the superiority of high-precision numerical solutions. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Those skilled in the art will appreciate that the accompanying drawings are provided for a better understanding of the present invention and do not constitute any limitation on the scope of the present invention.

[0064] Figure 1a A schematic diagram of a background grid in an existing two-dimensional overlapping grid system;

[0065] Figure 1b Schematic diagram of the body-fitting grid in the existing two-dimensional overlapping grid system;

[0066] Figure 1c is a schematic diagram of an existing two-dimensional overlapping grid system;

[0067] Figure 2 Schematic diagram of the existing traditional interpolation method based on inverse distance weighting;

[0068] Figure 3 A flowchart of a THINC-based overlapping grid VOF interpolation method according to an embodiment of the present invention;

[0069] Figure 4a-4c A schematic diagram of performing VOF field interpolation on overlapping grids according to an embodiment of the present invention;

[0070] Figure 5 A schematic diagram of a two-dimensional overlapping grid system provided by an embodiment of the present invention;

[0071] Figure 6 A schematic diagram of a three-dimensional overlapping grid system provided by an embodiment of the present invention;

[0072] Figure 7 A schematic diagram of a two-dimensional star profile interpolation example provided by an embodiment of the present invention;

[0073] Figure 8 A schematic diagram of a calculation example of three-dimensional star profile interpolation provided by an embodiment of the present invention;

[0074] Figure 9 A visualization interface after overlapping grid interpolation processing provided by an embodiment of the present invention;

[0075] Figure 10 A schematic diagram of a distribution curve of VOF values ​​along the x-axis provided in one embodiment of the present invention;

[0076] Figure 11 A schematic diagram of an overlapping grid system setup for a two-dimensional shear deformation flow example provided by an embodiment of the present invention;

[0077] Figure 12 A schematic diagram of a visualization interface for calculating a two-dimensional shear deformation flow example using a traditional algorithm on a static overlapping grid provided by an embodiment of the present invention;

[0078] Figure 13 A schematic diagram of a visualization interface for calculating a two-dimensional shear deformation flow example using the interpolation method of the present invention on a static overlapping grid provided by one embodiment of the present invention;

[0079] Figure 14 A schematic diagram of an overlapping grid system setup for a three-dimensional spherical deformation flow example provided by an embodiment of the present invention;

[0080] Figure 15 A schematic diagram of a visualization interface for calculating a three-dimensional spherical deformation flow example at time t=T / 2 using the interpolation method of the present invention on a static overlapping grid provided by one embodiment of the present invention;

[0081] Figure 16 A schematic diagram of a visualization interface for calculating a three-dimensional spherical deformation flow example at time t=T using the interpolation method of the present invention on a static overlapping grid provided by one embodiment of the present invention;

[0082] Figure 17 A schematic diagram of a visualization interface for calculating a three-dimensional spherical deformation flow example at time t=T / 2 using a traditional algorithm on a static overlapping grid provided by an embodiment of the present invention;

[0083] Figure 18A schematic diagram of a visualization interface for calculating a three-dimensional spherical deformation flow example at time t=T using a traditional algorithm on a static overlapping grid provided by an embodiment of the present invention.

[0084] in:

[0085] 1-hole unit; 2-interpolation boundary unit; 3-calculation unit; 4-interpolation template;

[0086] 100 - Interpolation boundary element; 110 - Gaussian integration point; 200 - Interface in contributing element; 300 - Interpolation template; 310 - Nearest contributing element; 400 - Background grid; 410 - Overlapping grid; 420 - Visualization interface. DETAILED DESCRIPTION

[0087] In order to make the purpose, advantages and features of the present invention clearer, the present invention is further described in detail below in conjunction with the accompanying drawings and specific embodiments. It should be noted that the drawings are in a very simplified form and use non-precise proportions, which are only used to conveniently and clearly assist in explaining the purpose of the embodiments of the present invention. In order to make the purpose, features and advantages of the present invention more obvious and easy to understand, please refer to the accompanying drawings. It should be noted that the structures, proportions, sizes, etc. illustrated in the drawings of this specification are only used to match the contents disclosed in the specification for people familiar with this technology to understand and read, and are not used to limit the conditions for the implementation of the present invention. Any modification of the structure, change in the proportional relationship or adjustment of the size, under the condition that the effect produced by the present invention and the purpose that can be achieved are the same or similar, should still fall within the scope of the technical content disclosed by the present invention.

[0088] As used in the present invention, the singular forms "a", "an" and "the" include plural objects, unless the content clearly indicates otherwise. As used in the present invention, the term "or" is generally used in a sense including "and / or", unless the content clearly indicates otherwise. As used in the present invention, the term "several" is generally used in a sense including "at least one", unless the content clearly indicates otherwise. As used in the present invention, the term "at least two" is generally used in a sense including "two or more", unless the content clearly indicates otherwise. In addition, the terms "first", "second" and "third" are used for descriptive purposes only and are not to be understood as indicating or implying relative importance or implicitly indicating the number of technical features indicated. Therefore, the features defined as "first", "second" and "third" may explicitly or implicitly include one or at least two of the features.

[0089] In the description of the present invention, unless otherwise expressly specified or limited, the terms "mounted," "connected," "connected," and "fixed" should be understood in a broad sense. For example, they may refer to fixed connection, detachable connection, or integration; mechanical connection or electrical connection; direct connection or indirect connection through an intermediate medium; and internal communication between two components or interaction between two components. Those skilled in the art will understand the specific meanings of the above terms in the present invention in specific circumstances.

[0090] The background art has described in detail how to use the THINC / QQ interface capture algorithm to update the fluid volume fraction field. The core of this algorithm is to perform high-precision segmented reconstruction of the gas-liquid interface by combining quadratic polynomials with THINC functions.

[0091] Based on this, in order to accurately interpolate the VOF field on an unstructured overlapping grid system, an embodiment of the present invention provides an overlapping grid VOF interpolation method based on THINC, which includes:

[0092] S1. Obtain the body center coordinates and Gaussian integration point coordinates of each grid cell in the overlapping grid system;

[0093] S2. Obtaining grid assembly information in the overlapping grid system, wherein the grid assembly information includes an interpolation boundary unit index list and an index list of contribution units in the interpolation template corresponding to each interpolation boundary unit;

[0094] S3, reading in the interface quadratic polynomial coefficients, interface positions, and fluid volume fraction values ​​of each contributing unit;

[0095] S4. Select an interpolation boundary unit from the interpolation boundary unit index list for processing, and select a corresponding interpolation template;

[0096] S5. For the selected interpolation boundary unit, select one Gaussian integral point of the interpolation boundary unit for processing;

[0097] S6. Calculate the distance from the Gaussian integral point to the center of each contributing unit grid in the corresponding interpolation template, and find the contributing unit closest to the Gaussian integral point;

[0098] S7. Determine whether the nearest contributing unit is an interface unit. If the nearest contributing unit is not the interface unit, the fluid volume fraction value of the Gaussian integral point is the fluid volume fraction value of the nearest contributing unit. If the nearest contributing unit is the interface unit, substitute the displacement between the Gaussian integral point and the center of the grid of the nearest contributing unit into the THINC function of the nearest contributing unit to calculate the fluid volume fraction value of the Gaussian integral point.

[0099] S8, repeating S5-S7, calculating the fluid volume fraction values ​​of all Gaussian integral points of the selected interpolation boundary unit, and calculating the fluid volume fraction value of the selected interpolation boundary unit using Gaussian integral method;

[0100] S9. Repeat S4-S8 to calculate the fluid volume fraction values ​​of all interpolation boundary cells in the overlapping grid system, thereby realizing overlapping grid interpolation calculation of the fluid volume fraction field.

[0101] During the implementation process, the present invention makes use of the cellCellStencil class in the OpenFOAM function library to effectively obtain the interpolation template and contribution unit information corresponding to each interpolation boundary cell after the overlapping grid is moved. This class efficiently stores the topological relationship between the interpolation boundary cell and its contribution unit in the form of a dynamic array, thereby providing support for VOF overlapping grid interpolation based on the THINC function. Specifically, with the help of the member functions of the cellCellStenci class, the interpolation process can be finely controlled. For example, by calling the interpolationCells() function, the index list of the interpolation boundary cell can be obtained; and by using the donorStencil() function, the index list of the contribution unit corresponding to each interpolation boundary cell can be easily accessed.

[0102] Because the THINC function contains relevant geometric information about the gas-liquid interface, such as the normal direction and curvature, the THINC-based overlapping mesh VOF interpolation method proposed in this example can accurately transfer VOF field information between overlapping meshes. This method not only ensures the sharpness of the interface and the conservation of fluid mass, but also eliminates the need for explicit flux corrections.

[0103] Please refer to the Figure 3 , and combined with Figure 4a-4c ,exist Figure 4a In the figure, the red border represents the interpolation boundary unit 100, the red dot represents the Gaussian integration point 110 in the interpolation boundary unit 100, and the blue curve represents the interface 200 in the contribution unit; Figure 4b In FIG, the gray border represents the interpolation template 300 , the red dot represents the selected Gaussian integration point 110 , and the black border represents the nearest contributing unit 310 .

[0104] First, execute S1-S3 to obtain the necessary geometric information for each grid cell in the overlapping grid system, including the body-center coordinates and Gaussian integration point 110 coordinates of each grid cell; obtain the grid assembly information for the overlapping grid system, including a list of interpolation boundary cell indices and a list of contributing cells in the interpolation template 300 corresponding to each interpolation boundary cell 100; and read in the interface quadratic polynomial coefficients, interface position, and fluid volume fraction value for each contributing cell. All of this information can be obtained using known techniques.

[0105] Then, S4 and S5 are executed to select an interpolation boundary cell 100 from the interpolation boundary cell index list for processing, and a corresponding interpolation template 300 is selected. For the selected interpolation boundary cell 100, one of the Gaussian integral points 110 of the interpolation boundary cell 100 is selected for processing. In this process, the method for selecting the interpolation boundary cell 100 and the Gaussian integral point 110 of the interpolation boundary cell 100 is not specified.

[0106] Then, S6 is executed to calculate the distance from the Gaussian integral point 110 to the center of each contributing unit mesh in the corresponding interpolation template 300, and find the contributing unit closest to the Gaussian integral point 110. In this embodiment, in S6, the distance r from the Gaussian integral point 110 to the center of each contributing unit mesh in the corresponding interpolation template 300 is calculated. ci The specific formula is as follows:

[0107] r ci =√(x gi -x ci ) 2 +(y gi -y ci ) 2 +(z gi -z ci ) 2 (11)

[0108] In the formula, (x gi ,y gi , z gi ) is the coordinate of the Gaussian integration point 110, (x ci ,y ci , z ci ) is the coordinate of the center of the contributing unit grid.

[0109] Then, S7 is executed to determine whether the nearest contributing unit 310 is an interface unit. If the nearest contributing unit 310 is not the interface unit, the fluid volume fraction value of the Gaussian integral point 110 is selected as the fluid volume fraction value of the nearest contributing unit 310. If the nearest contributing unit 310 is an interface unit, the displacement between the Gaussian integral point 110 and the grid center of the nearest contributing unit 310 is substituted into the THINC function of the nearest contributing unit 310 to calculate the fluid volume fraction value of the Gaussian integral point 110.

[0110] In this embodiment, when ∈≤φ i ≤1-∈,∈=10 -8 When , the nearest contributing unit 310 is an interface unit, otherwise it is not; where φ i is the fluid volume fraction value of the nearest contributing unit 310 .

[0111] In this embodiment, if the nearest contributing element 310 is not the interface element, the fluid volume fraction value of the Gaussian integral point 110 is selected as the fluid volume fraction value of the nearest contributing element 310 , and the specific formula is as follows:

[0112]

[0113] If the nearest contributing element 310 is the interface element, the specific formula for calculating the fluid volume fraction value at the Gaussian integral point 110 is as follows:

[0114]

[0115] X=x gi -x ci , Y=y gi -y ci , Z=z gi -z ci (15)

[0116] In the formula, (x gi ,y gi , z gi ) is the coordinate of the Gaussian integration point 110, (x ci ,y ci , z ci ) is the coordinate of the center of the contributing unit grid; C 200 、C 020 、C 002 、C 110 、C 101 、C 011 、C 100 、C 010 、C 001The coefficients of the interface quadratic polynomial read from S3.

[0117] Then, S8 is executed to repeat S5-S7 to calculate the fluid volume fraction values ​​of all Gaussian integral points 110 of the selected interpolation boundary unit 100, and the fluid volume fraction value of the selected interpolation boundary unit 100 is calculated using the Gaussian integral method.

[0118] In this embodiment, in S8, the formula for calculating the fluid volume fraction value of the selected interpolation boundary unit 100 using the Gaussian integral method is as follows:

[0119]

[0120] Where, is the weight of the fluid volume fraction value at each Gaussian integration point 110 .

[0121] Finally, S9 is executed to repeat S4-S8 to calculate the fluid volume fraction values ​​of all interpolation boundary cells 100 in the overlapping grid system, thereby realizing the overlapping grid interpolation calculation of the fluid volume fraction field.

[0122] In order to verify the numerical accuracy, the focus is on the actual application effect of the VOF interpolation algorithm on static and dynamic overlapping grids. In this embodiment, by selecting a number of representative examples, key performance indicators such as interpolation accuracy and algorithm convergence are systematically and comprehensively analyzed. In addition, to ensure the accuracy and reliability of the evaluation results, the relative mass error E is defined. m and relative position error E p As the core indicator for measuring numerical error, its calculation formula is:

[0123]

[0124]

[0125] Where, φ ni and |φ ei denote the numerical solution and the exact solution respectively, and N denotes the total number of grid cells.

[0126] In order to verify the performance of the interpolation algorithm proposed in this embodiment in terms of accuracy and conservation, Figure 5-Figure 6 A series of numerical calculations are performed on the static and dynamic overlapping grids shown. Figure 5-Figure 6In the figure, blue represents background mesh 400, and red represents overlay mesh 410. An overlay mesh system consisting of background mesh 400 and overlay mesh 410 was employed. Specifically, the two-dimensional overlay mesh system includes background mesh 400 composed of regular quadrilateral elements and overlay mesh 410 composed of isosceles right triangle elements. The three-dimensional overlay mesh system consists of background mesh 400 composed of regular hexahedral elements and overlay mesh 410 composed of isosceles right prism elements. All meshes used in the calculations were generated using GMSH software. Isosceles right triangle elements and isosceles right prism elements have the greatest non-orthogonality, and the element sizes of background mesh 400 and overlay mesh 410 are consistent to ensure accuracy and comparability of the calculations.

[0127] To verify the advantages of the newly developed algorithm, the numerical accuracy of the VOF interpolation algorithm was verified on two-dimensional and three-dimensional static overlapping grids. The two-dimensional star profile example was first proposed by "Ii S, Sugiyama K, Takeuchi S, et al. An Interface Capturing Method with a Continuous Function: The THINC Method with Multi-Dimensional Reconstruction [J]. Journal of Computational Physics, 2012, 231 (5): 2328-2358". For the two-dimensional star profile interpolation problem, the computational domain and the background grid area are set to the unit square area [0,1] 2 , and the overlapping grid area is defined as [0.25,0.75] 2 On this basis, the VOF field is initialized by the following indicator function:

[0128]

[0129]

[0130] Set the center coordinates of the star outline to (x c ,y c )=(0.5,0.5). In addition, a 3D spherical contour example is included to verify the numerical accuracy of the 3D overlapping grid interpolation method. Specifically, the computational domain and the background grid area are both set to [0,1] 3 , and the overlapping grid area is set to [0.25,0.75] 3 On this basis, the VOF field is initialized according to the following indicator function:

[0131]

[0132]

[0133] The coordinates of the center of the sphere are set as (x c ,y c ,z c )=(0.5,0.5,0.5). Figure 7-Figure 8 The visualization interface of 2D / 3D overlapping grid interpolation examples is shown. Figure 7-Figure 8 In the figure, the visualization interface 420 is presented in blue, while the background grid 400 and the overlapping grid are marked in gray and green respectively. Figure 5-Figure 6 The unit types shown remain consistent.

[0134] In order to verify the effectiveness and advantages of the algorithm proposed in this embodiment, Figure 9 In 128 2 At the grid resolution, the visualization interface 420 after overlapping grid interpolation (VOF = 0.5) is shown. Observation shows that the visualization interface 420 in the overlapping grid is highly consistent with the interface in the background grid 400. For further analysis, the present invention selects the edge cell area, specifically at y = 0.25 and with the x-axis range between 0.25 and 0.75, and plots the distribution curve of the VOF value along the x-axis, as shown in the figure. Figure 10 As shown. Figure 10 In the figure, the black dots represent the exact solution, the blue triangles and the red squares represent the numerical solutions obtained by the traditional interpolation method and the interpolation method provided by the present invention, respectively. By comparing the numerical results of the interpolation method of the present invention and the traditional overlapping grid interpolation method, it can be found that the traditional interpolation algorithm easily causes numerical dissipation of the VOF value at the interface. The interpolation method provided by the present invention, based on the THINC function, accurately records the interface geometry information through the quadratic surface polynomial, thereby effectively maintaining the sharpness of the interface. Therefore, the VOF field obtained after overlapping grid interpolation using the interpolation method provided by the present invention is closer to the exact solution, further verifying the advantages of the algorithm.

[0135] The THINC / QQ interface capture algorithm can be combined with the overlapping grid interpolation algorithm based on the THINC function to develop an interface capture algorithm suitable for overlapping grids. In order to verify the ability of the interpolation algorithm provided by the present invention to capture interface distortion, stretching and large deformation on static overlapping grids, two-dimensional shear deformation flow and three-dimensional spherical deformation flow were designed as numerical verification examples. At the same time, in order to conduct an in-depth comparative analysis, the numerical results obtained by using the MULES interface capture algorithm combined with the inverse distance weighted interpolation algorithm were also statistically analyzed. It is worth noting that the overInterDyMFoam solver in OpenFOAM also uses MULES and inverse distance weighted interpolation as standard interface capture and overlapping grid interpolation algorithms.

[0136] The paper "Rider WJ, Kothe DB. Reconstructing Volume Tracking [J]. Journal of Computational Physics, 1998, 141 (2): 112-152" first introduced a two-dimensional shear deformation flow example, which covers the large deformation and stretching phenomena of the interface. The two-dimensional shear deformation flow example is set up as follows Figure 11 In the figure, the visualization interface 420 (VOF = 0.5) at the initial time (t = 0) is marked in blue, the background grid 400 and the overlapping grid 410 are distinguished by gray and green respectively, and the grid unit type is the same as Figure 5 Keep it consistent. The calculation domain and the background grid 400 area are both set to [0,1] 2 , and the overlapping grid 410 area is set to [0.25, 0.75] × [0.0625, 0.75]. At t = 0, the liquid phase occupies a circle with a radius of r = 0.15 and a center coordinate (x c ,y c )=(0.5,0.75). In addition, the velocity field of this example is accurately described by a stream function:

[0137]

[0138] The initial VOF field deforms under the influence of the velocity field and is gradually stretched into a spiral shape with a slender tail. This process lasts for half a cycle, t = T / 2. Subsequently, the direction of the velocity field reverses, allowing the VOF field to return to its initial state at t = T.

[0139] In order to verify the numerical accuracy of the interpolation method provided by the present invention, the numerical results obtained by the VOF equation solver combining the traditional inverse distance weighted interpolation and the MULES interface capture algorithm are compared. Figure 12-13 The comparison of the visualization interface 420 (VOF=0.5) of a two-dimensional shear deformation flow case calculated using different algorithms on a static overlapping grid is shown. The red and blue lines in the figure represent the numerical solutions at times t=T / 2 and t=T, respectively, while the black line represents the exact solution at time t=T. At t=T / 2, the interface undergoes significant deformation and is stretched into a spiral with an extremely long and thin tail. It is worth noting that because the distorted interface becomes extremely thin, it exceeds the resolvability of the finite resolution grid, and the tail of the spiral VOF field is easily broken into small droplets. Figure 13In the numerical simulation, the interpolation method provided by the present invention was applied. The results show that the interpolation method provided by the present invention can fully resolve the slender tail at t = T / 2 and show a more refined interface structure. Finally, the interface returns to a circular shape, which is highly consistent with the initial contour. In contrast, Figure 12 The results of numerical simulations using traditional algorithms are shown. Due to the numerical dissipation introduced by the MULES algorithm and inverse distance weighted interpolation, the distorted interface is smoothed. However, the elongated tail is severely damaged, resulting in inaccurate alignment between the visualized interface 420 on the overlay mesh and the interface on the background mesh. In this case, the interface loses its geometric characteristics, making it difficult to restore it to its initial state at time t = T.

[0140] In order to verify the numerical accuracy of the interpolation method provided by the present invention in a three-dimensional overlapping grid system, a three-dimensional spherical deformation flow verification example is designed. This example was originally proposed by Enright et al.

[12] . At the initial time t = 0, the liquid phase is confined to a radius r = 0.15 and the coordinates of the sphere center are located at (x c ,y c ,z c )=(0.35,0.35,0.35). The velocity field expression of this example is:

[0141]

[0142] In the formula, G(t) = cos(πt / T) is defined as a cosine function that varies with time. In the three-dimensional spherical deformation flow example, the overlapping grid system is set as follows Figure 14 As shown. Figure 14 In the figure, blue is used to represent the visualization interface 420 (i.e., the isosurface of VOF=0.5), the background grid 400 and the overlapping grid 410 are distinguished by gray and green respectively, and the cell types of these two types of grids are the same. Figure 5 In addition, the computational domain and background grid 400 area of ​​this example are set to [0,1] 3 , and the area of ​​the overlapping grid 410 is specifically defined as [0.25,0.75] 3 . In the example of three-dimensional spherical deformation flow, the liquid phase sphere is strongly stretched by two vortices. At half a cycle t=T / 2, the shear rate of the sphere reaches its maximum value, the interface is stretched to the limit, and a thin film is formed. If the grid resolution is insufficient, the details of the film cannot be accurately captured. Subsequently, affected by the reverse velocity field, the deformed sphere will gradually return to its original shape. The present invention uses an interpolation algorithm based on the THINC / QQ interface capture format and the THINC function for numerical calculations.

[0143] Figures 15-18The visualization interface 420 (VOF=0.5) of the three-dimensional spherical deformation flow at two times t=T / 2 and t=T on the static overlapping grid 410 is depicted. Figures 15-18 , the visual interface 420 of the background grid 400 and the overlapping grid 410 are distinguished by blue and red respectively. It can be observed that compared with the traditional algorithm, the interpolation method provided by the present invention exhibits a stronger interface capture capability, which can not only accurately capture the thin film structure of the VOF field, but also effectively restore the VOF field to its original spherical shape. From the numerical results, the performance of the interpolation method provided by the present invention is significantly better than the traditional method. The interface obtained by the interpolation method provided by the present invention is sharper, and the thin film structure of the interface is finer. It is worth mentioning that the interpolation method provided by the present invention can maintain the thickness of the interface, so that the interface details of the thin film part are clearer and more visible than the traditional algorithm. In addition, since the THINC function effectively retains the geometric information of the interface during the calculation process, the visual interface 420 on the overlapping grid 410 and the background grid 400 shows a high degree of consistency. In contrast, the traditional algorithm has obvious numerical dissipation when processing thin films, which makes the VOF field unable to accurately restore to its initial spherical shape. In summary, the numerical results fully verify the excellent ability of the interpolation method provided by the present invention in capturing three-dimensional complex motion interfaces on overlapping grid systems and the superiority of high-precision numerical solutions.

[0144] Based on the same inventive concept, an embodiment of the present invention further provides a readable storage medium having a computer program stored thereon. When the computer program is executed, the THINC-based overlapping grid 410 VOF interpolation method described above can be implemented.

[0145] The readable storage medium can be a tangible device that can keep and store the instruction used by the instruction execution device, such as, but not limited to, an electrical storage device, a magnetic storage device, an optical storage device, an electromagnetic storage device, a semiconductor storage device or any suitable combination thereof. The more specific example (non-exhaustive list) of readable storage medium includes: a portable computer disk, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), a static random access memory (SRAM), a portable compact disc read-only memory (CD-ROM), a digital versatile disk (DVD), a memory stick, a floppy disk, a mechanical encoding device, such as a punch card or a convex structure in a groove having instructions stored thereon, and any suitable combination thereof. The computer program described herein can be downloaded to each computing / processing device from the readable storage medium, or downloaded to an external computer or external storage device by a network, such as the Internet, a local area network, a wide area network and / or a wireless network. The network can include copper transmission cables, optical fiber transmission, wireless transmission, routers, firewalls, switches, gateway computers and / or edge servers. The network adapter card or network interface in each computing / processing device receives the computer program from the network and forwards the computer program for storage in a readable storage medium in each computing / processing device. The computer program for performing the operations of the present invention can be assembly instructions, instruction set architecture (ISA) instructions, machine instructions, machine-dependent instructions, microcode, firmware instructions, state 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 programming languages. The computer program can be executed entirely on the user's computer, partially on the user's computer, as a separate software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In the case of a remote computer, the remote computer can be connected to the user's computer through any type of network, including a local area network (LAN) or a wide area network (WAN), or can be connected to an external computer (for example, through the Internet using an Internet service provider). In some embodiments, various aspects of the present invention are implemented by utilizing state information of a computer program to personalize an electronic circuit, such as a programmable logic circuit, a field programmable gate array (FPGA), or a programmable logic array (PLA), which can execute computer-readable program instructions.

[0146] Various aspects of the present invention are described herein with reference to the flowcharts and / or block diagrams of the methods, systems, and computer program products according to embodiments of the present invention. It should be understood that each block of the flowcharts and / or block diagrams, as well as the combination of blocks in the flowcharts and / or block diagrams, can be implemented by computer programs. These computer programs can be provided to a processor of a general-purpose computer, a special-purpose computer, or other programmable data processing device, thereby producing a machine such that when these programs are executed by the processor of the computer or other programmable data processing device, a device is generated that implements the functions / actions specified in one or more blocks in the flowcharts and / or block diagrams. These computer programs can also be stored in a readable storage medium, which causes the computer, programmable data processing device, and / or other device to operate in a specific manner, so that the readable storage medium storing the computer program comprises an article of manufacture comprising instructions for implementing various aspects of the functions / actions specified in one or more blocks in the flowcharts and / or block diagrams.

[0147] The computer program may also be loaded onto a computer, other programmable data processing apparatus, or other device so that a series of operational steps are performed on the computer, other programmable data processing apparatus, or other device to produce a computer-implemented process, thereby causing the computer program executed on the computer, other programmable data processing apparatus, or other device to implement the functions / actions specified in one or more blocks in the flowchart and / or block diagram.

[0148] The above description is merely a description of preferred embodiments of the present invention and does not limit the scope of the present invention. Any changes or modifications made by persons skilled in the art based on the above disclosure are within the scope of protection of the present invention. Obviously, various modifications and variations may be made by persons skilled in the art without departing from the spirit and scope of the present invention. Thus, if such modifications and variations fall within the scope of the present invention and its equivalents, the present invention is intended to include such modifications and variations.

Claims

1. A THINC-based overlapping grid VOF interpolation method, characterized in that: include: S1. Obtain the body center coordinates and Gaussian integration point coordinates of each grid cell in the overlapping grid system; S2. Obtaining grid assembly information in the overlapping grid system, wherein the grid assembly information includes an interpolation boundary unit index list and an index list of contribution units in the interpolation template corresponding to each interpolation boundary unit; S3, reading in the interface quadratic polynomial coefficients, interface positions, and fluid volume fraction values ​​of each contributing unit; S4. Select an interpolation boundary unit from the interpolation boundary unit index list for processing, and select a corresponding interpolation template; S5. For the selected interpolation boundary unit, select one Gaussian integral point of the interpolation boundary unit for processing; S6. Calculate the distance from the Gaussian integral point to the center of each contributing unit grid in the corresponding interpolation template, and find the contributing unit closest to the Gaussian integral point; S7. Determine whether the nearest contributing unit is an interface unit. If the nearest contributing unit is not the interface unit, the fluid volume fraction value of the Gaussian integral point is the fluid volume fraction value of the nearest contributing unit. If the nearest contributing unit is the interface unit, substitute the displacement between the Gaussian integral point and the center of the grid of the nearest contributing unit into the THINC function of the nearest contributing unit to calculate the fluid volume fraction value of the Gaussian integral point. S8, repeating S5-S7, calculating the fluid volume fraction values ​​of all Gaussian integral points of the selected interpolation boundary unit, and calculating the fluid volume fraction value of the selected interpolation boundary unit using Gaussian integral method; S9. Repeat S4-S8 to calculate the fluid volume fraction values ​​of all interpolation boundary cells in the overlapping grid system, thereby realizing overlapping grid interpolation calculation of the fluid volume fraction field.

2. The THINC-based overlapping grid VOF interpolation method according to claim 1, characterized in that: In S6, the distance r from the Gaussian integral point to the center of each contributing unit grid in the corresponding interpolation template is calculated. ci The specific formula is as follows: In the formula, (x gi ,y gi , z gi ) is the coordinate of the Gaussian integral point, (x ci ,y ci , z ci ) is the coordinate of the center of the contributing unit grid.

3. The THINC-based overlapping grid VOF interpolation method according to claim 1, characterized in that: When ∈≤φ i ≤1-∈,∈=10 -8 When , the nearest contributing unit is the interface unit, otherwise it is not; where φ i is the fluid volume fraction value of the nearest contributing unit.

4. The THINC-based overlapping grid VOF interpolation method according to claim 1, characterized in that: In S7, if the nearest contributing element is the interface element, the specific formula for calculating the fluid volume fraction value of the Gaussian integral point is as follows: P i (x gi )=C 200 X 2 +C 020 Y 2 +C 002 Z 2 +C 110 XY+C 101 XZ+C 011 YZ +C 100 X+C 010 Y+C 001 Z X=x gi -x ci ,Y=y gi -y ci ,Z=z gi -z ci In the formula, (x gi ,y gi , z gi ) is the coordinate of the Gaussian integral point, (x ci ,y ci , z ci ) is the coordinate of the center of the contributing unit grid; C 200 、C 020 、C 002 、C 110 、C 101 、C 011 、C 100 、C 010 、C 001 are the coefficients of the interface quadratic polynomial.

5. The THINC-based overlapping grid VOF interpolation method according to claim 4, characterized in that: In S8, the formula for calculating the fluid volume fraction value of the selected interpolation boundary unit using the Gaussian integral method is as follows: Where, Weights for the fluid volume fraction values ​​at each Gaussian integration point.

6. A readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed, the THINC-based overlapping grid VOF interpolation method according to any one of claims 1 to 5 can be implemented.

Images