A horizontal piecewise linear reconstruction method for two-dimensional two-phase interfaces

Abstract

The present invention relates to the field of two-phase flow technology, and particularly relates to a two-dimensional two-phase interface horizontal piecewise linear reconstruction method: the interface direction is preliminarily calculated based on the volume fraction of the main phase, and a local coordinate system is established accordingly to ensure that the main phase is below the interface; the initial intercept of the interface equation is set to obtain the initial interface equation; a series of piecewise horizontal planes are used to approximate the phase interface, and the interface equation is used to calculate the horizontal plane height and calculate the volume fraction of the main phase; the interface equation is iteratively obtained based on the difference between the statistical result and the true result to complete the first interface reconstruction; based on the reconstructed interface, the LS function is calculated using a geometric method, and the interface direction is recalculated using the LS function, and then the above process is repeated to achieve accurate interface reconstruction. This patent obtains a unique interface in the local coordinate system and solves the VOF equation through four arithmetic operations, greatly reducing the difficulty of interface reconstruction and providing technical support for numerical research on two-phase problems.

Description

Technical Field

[0001] The present invention belongs to the technical field of two-phase flow, and particularly relates to a two-dimensional two-phase interface horizontal segmented linear reconstruction method. Background Art

[0002] Two-phase flow refers to a flow process in which two media interact with each other. Two-phase flow is very common in our lives, such as soap bubbles formed during the laundry process; the movement of rainwater in the air, etc. A large number of two-phase flow problems are also involved in industrial application fields, such as oil-water two-phase flow formed by water injection to drive oil during oil extraction; boiling heat transfer widely used in nuclear energy, aerospace, and thermal management of electronic devices, where the liquid coolant undergoes a violent phase change on the heat transfer surface and generates a large number of bubbles to destroy the thermal boundary layer, achieving high intensity. For the above two-phase flow problems, the phase interface usually freely deforms with the flow of the fluid. Accurately capturing the phase interface is the key to studying these two-phase problems and also the difficulty in the numerical study of two-phase flow problems. Currently, the commonly used phase interface tracking methods in numerical calculations are the VOF method, the Level set (LS) method, and the VOSET method. The principles, advantages, and disadvantages of them will be briefly introduced below, and the technical difficulties existing in the current phase interface tracking methods will be analyzed.

[0003] As Figure 2 shown, using the method of the volume of fluid (VOF function), the volume fraction C of the fluid is used to represent the volume fraction of the main phase in a unit grid, and its value ranges from 0 to 1. C = 0 means that the current grid does not contain the main phase; C = 1 means that the current grid is filled with the main phase; 0 < C < 1 means that the current grid contains a vapor-liquid interface. When reconstructing the interface using the VOF function, the mass conservation of the solution object can be guaranteed, but the calculation accuracy of the interface direction and curvature is low, resulting in distortion of the reconstructed interface and large surface tension errors.

[0004] The Level-Set function (abbreviated as LS function) uses a continuous distance function φ to represent the spatial distance from a fluid node to the two-phase interface, and its zero level set is the two-phase interface. Positive and negative signs are used to represent different phases. When the distance function is less than 0, it represents one phase, and when it is greater than 0, it represents the other phase. The LS function has continuous characteristics and can obtain high-precision interface direction and curvature, but it is difficult to ensure the mass conservation of the research object.

[0005] The VOSET method (coupled Volume-of-Fluid and Level Set method) combines the VOF function and the LS function to achieve complementary advantages. The steps are as follows: (1) First, use the initial VOF function to construct the initial phase interface; (2) Then, using the information of the phase interface, use the geometric method to solve the LS function, and use the LS function to reconstruct the phase interface; (3) Repeat step (2) several times to construct an accurate phase interface. The VOSET method inherits the advantages of the VOF method and the LS function method, thus being able to calculate the accurate interface curvature and maintain the conservation of the main phase.

[0006] The above VOF and VOSET methods both use piecewise linear interface calculation (PLIC) to reconstruct the phase interface, that is, the grid interface is approximately replaced by the "straight line instead of curved line" method, and the phase interface within the grid uses a straight line to replace the real curved surface phase interface. The PLIC method is more complicated to reconstruct the plane and interface advancement. In terms of interface reconstruction, it is necessary to solve the fourth-order equation of formula (12), and it is necessary to calculate n x =0 and n y = 0 (such as Figure 3 ) for additional discussion; in terms of interface advancement, multiple complex transformations of formula (12) are required to obtain the volume fraction flux on the grid surface, and then update the VOF function. Taking the right interface flux as an example, the special case of the PLIC piecewise linear reconstruction method is as follows Figure 4 As shown, formulas (12) and (13) are transformed to obtain the new interface equation formula (14) and the right interface flux formula (15).

[0007]

[0008] In the above formula, α is the intercept of the line, and H(x) is the Heaviside function, which is defined as follows:

[0009]

[0010]

[0011]

[0012]

[0013] In summary, the PLIC method for interface reconstruction involves tedious and difficult-to-understand operations such as solving quartic equations, classifying special cases, and transforming equations, making it difficult for researchers to apply this method to two-phase problems. Therefore, it is urgent to develop an interface tracking method that balances computational efficiency and operational complexity while ensuring accuracy. This method can accurately describe the interface of two-phase problems, facilitate its implementation in numerical simulations, and provide technical support for the study of two-phase flows. Summary of the Invention

[0014] The object of the present invention is to provide a two-dimensional two-phase interface horizontal piecewise linear reconstruction method, characterized in that it includes the following steps:

[0015] S1: Grid division of the calculation area;

[0016] S2: For the grid (i, j) where the two-phase interface is located, the interface direction n is preliminarily calculated based on the volume fraction of the main phase. x i,j and n y i,j , and establish a local coordinate system according to the phase interface direction to ensure that the main phase is below the phase interface, where i and j are the coordinates of the grid;

[0017] S3: Set the initial interface through the center of the grid (i, j) to obtain the initial interface equation;

[0018] S4: Use R segmented horizontal surfaces to approximate the phase interface and use the interface equation to calculate the height of the segmented horizontal surface;

[0019] S5: Calculate the sum of the areas below the segmented horizontal plane, and obtain the estimated main phase volume fraction C based on the ratio of the sum of the areas to the unit grid area. 0 i,j :

[0020]

[0021] Where R is the number of segmented horizontal planes, S r is the area below the horizontal plane of a single segment, S g is the unit grid area size;

[0022] S6: Define the convergence criterion err for solving the interface equation and compare it with the estimated main phase volume fraction C 0 i,j and the actual main phase volume fraction C i,j , and calculate the difference D between the two:

[0023] D = C 0 i,j - C i,j (2)

[0024] When |D|≤err, execute S9;

[0025] When |D|>err, execute S7;

[0026] S7: Modify the interface equation. When D>0 in S6, modify the interface equation intercept B=BH; when

[0027] When D≤0 in S5, modify B=B+H;

[0028] Where B is the intercept of the interface equation; H is the adjustment parameter, and the initial value of H is the unit grid characteristic size;

[0029] S8: H = H / 2, repeat S4-S8 until |D| ≤ err in S6;

[0030] S9: Obtain the interface equation through S3-S8, and calculate the LS function using the geometric method based on this equation;

[0031] S10: Recalculate the phase interface direction n using the LS function x i,j and n y i,j , then repeat S2-S9 N times to complete multiple interface reconstructions, obtain the final LS function and interface equation, obtain the position and shape of the phase interface, and then describe the effect of the phase interface on flow and heat transfer;

[0032] S11: Solve the VOF function using a geometric method, calculate the flux of the main phase volume fraction at the interface of the current grid (i, j), add them together to get the updated main phase volume fraction C, and use C to update the actual main phase volume fraction C of the next time layer i,j , to achieve the description of phase interface migration and obtain the phase distribution of the entire field;

[0033] S12: Repeat S1-S11 until the preset simulation time is reached.

[0034] The phase interface direction n in S2 x i,j and n y i,j is defined as:

[0035]

[0036]

[0037] Where: C i,j represents the volume fraction of the main phase of the grid (i, j), and Δx and Δy represent the unit grid length in the x and y directions.

[0038] The value range of R in S4 is: 40≥R≥10.

[0039] The convergence criterion for solving the interface equation of S6 is err=10e-5.

[0040] The specific steps of repeating S10 N times S2-S9 to complete multiple interface reconstructions are as follows: set N=3; use the interface reconstructed by S2-S9 for the first iteration to solve the LS function, and use formulas (5) and (6) in subsequent iterations to determine the direction of the phase interface, and finally obtain an accurate interface.

[0041]

[0042]

[0043] Where: n x i,j and n y i,j represents the components of the interface direction vector in the x and y directions; φ i,j Represents the size of the LS function of the grid (i, j), Δx and Δy represent the unit grid length in the x and y directions.

[0044] The flux of the main phase volume fraction at the interface of the current grid (i, j) in S11 is CH t 、CH b 、CH l 、CH r , and the updated main phase volume fraction C is defined as:

[0045]

[0046]

[0047]

[0048]

[0049] C=C0+CH b -CH t +CH l -CH r (11)

[0050] Where: subscripts t, b, l and r represent the upper interface, lower interface, left interface and right interface of the current grid respectively; v t 、v b 、u l and u r are the velocities of the upper, lower, left, and right interfaces of the current grid; n is the number of rectangular regions divided below the interface; h is the number of rectangular regions divided below the interface;i is the height of the i-th rectangular area of the current grid; Δt is the time step; Δl is the width of each rectangular area; CH t is the flux of the upper interface volume fraction, CH b is the flux of the lower interface volume fraction, CH l is the flux of the left interface volume fraction, CH r is the flux of the volume fraction of the right interface; C0 represents the volume fraction of the previous time layer.

[0051] The beneficial effects of the present invention are:

[0052] The present invention discloses a two-dimensional two-phase interface horizontal piecewise linear reconstruction method, namely the VOSET-HOPLIRE method. Unlike the traditional VOSET-PLIC and VOF-PLIC methods, it can achieve high-precision phase interface reconstruction and phase interface evolution process solution without a series of complex equation solving, classification, equation transformation and other operations. In the process of determining the interface equation, the initial interface equation is set through the center point of the grid to determine the initial intercept B of the initial interface equation. Then, the idea of "partition reconstruction" is adopted to use a series of segmented horizontal planes to approximate the real phase interface, and the interface equation is used to solve the horizontal plane height and the area below, and then the volume fraction of the main phase is calculated; by comparing and calculating the main phase volume fraction and the real main phase volume fraction, the interface equation intercept is iterated according to the difference between the two to complete the first interface reconstruction; based on the reconstructed interface, the LS function of the points around the interface is calculated by a geometric method, and the interface direction n is recalculated using the LS function. x i,j and n y i,j , then repeat the above process for a second interface reconstruction; after multiple iterations, accurate interface reconstruction is finally achieved; based on the segmented horizontal reconstruction interface, simple arithmetic operations are used to solve the VOF function and achieve phase interface migration. Therefore, the VOSET-HOPLIRE method disclosed in this invention can significantly reduce the difficulty of interface reconstruction while ensuring computational accuracy, allowing more researchers to participate in the study of two-phase problems and providing technical support for the study of two-phase problems.

[0053] Comparative experiments were conducted on the mass and bubble deformation errors when solving shear bubbles using the traditional VOF-PLIC method, the VOSET-PLIC method, and a two-dimensional horizontal piecewise linear reconstruction method for two-phase interfaces (VOSET-HOPLIRE method) disclosed in the present invention. The experimental results show that the VOSET-HOPLIRE method disclosed in the present invention is far superior to the VOF-PLIC reconstruction method in terms of computational accuracy, and is almost the same as the VOSET-PLIC method. In terms of computational efficiency, under the same hardware environment, the VOSET-HOPLIRE method disclosed in the present invention takes 16.08 hours to calculate, while the VOSET-PLIC method takes 17.53 hours, which is a reduction in computational time. The average error of the VOSET-HOPLIRE method disclosed in the present invention remains within 3%, and the computational efficiency is improved by 9.0%, verifying the efficiency and accuracy of the two-dimensional horizontal piecewise linear reconstruction method disclosed in the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 The present invention discloses a flow chart of a horizontal piecewise linear reconstruction method for a two-dimensional two-phase interface;

[0055] Figure 2 Schematic diagram of the VOF-PLIC piecewise linear reconstruction method in the background technology;

[0056] Figure 3 Schematic diagram of the application of the PLIC piecewise linear reconstruction method in the background technology;

[0057] Figure 4 A schematic diagram of a special case of the PLIC piecewise linear reconstruction method in the background art;

[0058] Figure 5 This is a schematic diagram of the initial state of the interface reconstruction of the present invention;

[0059] Figure 6 This is a schematic diagram of the real interface approximated by a segmented horizontal plane in the present invention;

[0060] Figure 7 A schematic diagram of the geometric method used to advance the interface in the present invention;

[0061] Figure 8 This is a schematic diagram comparing the interface reconstruction results of the present invention;

[0062] Figure 9 Comparison of deformation errors of shear bubbles under three reconstruction methods;

[0063] Figure 10 Comparison of mass conservation errors of shear bubbles under three reconstruction methods;

[0064] Figure 11Comparison of pseudo-velocity distribution of a stationary bubble using three reconstruction methods;

[0065] Figure 12 Comparison of relative quality deviations under three reconstruction methods. DETAILED DESCRIPTION

[0066] The present invention provides a two-dimensional two-phase interface horizontal piecewise linear reconstruction method, and the present invention is further described in detail below with reference to the accompanying drawings. Figure 1 The embodiment of the present invention disclosed herein discloses a horizontal piecewise linear reconstruction method for a two-dimensional two-phase interface, specifically comprising the following steps:

[0067] S1: Grid division of the calculation area;

[0068] Aiming at the two-phase flow problem, a mathematical model is established and the finite volume method is used for differential solution.

[0069] S2: For the grid (i, j) where the two-phase interface is located, the interface direction n is preliminarily calculated based on the volume fraction of the main phase. x i,j and n y i,j , and establish a local coordinate system according to the phase interface direction to ensure that the main phase is below the phase interface, where i and j are the coordinates of the grid;

[0070] The phase interface direction n in S2 x i,j and n y i,j is defined as:

[0071]

[0072]

[0073] Where: C i,j represents the volume fraction of the main phase of the grid (i, j), Δx and Δy represent the unit grid length in the x and y directions;

[0074] S3: Set the initial interface through the center of the grid (i, j) to obtain the initial interface equation;

[0075] In the local coordinate system, using n x i,j and n y i,j And the coordinates of the grid center point are used to solve the initial interface equation, as well as the boundary

[0076] The initial value of the surface equation intercept B, the interface equation is expressed as:

[0077]

[0078] like Figure 5 As shown in the figure, the dotted line represents the initial phase interface. In the local coordinate system, the present invention uses n x i,j and n y i,j And the coordinates of the grid center point are used to solve the initial interface equation and the initial value of the interface equation intercept B. Then the interface equation is determined.

[0079] S4: Use R segmented horizontal surfaces to approximate the phase interface and use the interface equation to calculate the height of the segmented horizontal surface;

[0080] The interface equation is divided into R segmented horizontal planes, and the real phase interface is approximated by R segmented horizontal planes.

[0081] like Figure 6 As shown, in the local coordinate system, the center point coordinate X of the horizontal coordinate of the segmented horizontal plane is i , put it into the interface equation, and solve the height h of the segmented horizontal surface i .

[0082] S5: Calculate the sum of the areas below the segmented horizontal plane, and obtain the estimated main phase volume fraction C based on the ratio of the sum of the areas to the unit grid area. 0 i,j :

[0083]

[0084] Where R is the number of segmented horizontal planes, S r is the area below the horizontal plane of a single segment, S g is the unit grid area size;

[0085] Theoretically, the larger the value of R, the closer the segmented horizontal plane is to the real interface. After a large number of experiments and research on experimental data, the inventor pointed out that when R is too large, the calculation efficiency will decrease, but the accuracy improvement will be limited. On the other hand, when the R value is too small, the calculation accuracy cannot be guaranteed.

[0086] Therefore, in an optional embodiment, the inventor proposes that the value range of R in S4 is: 40≥R≥10,

[0087] It can strike a good balance between computational efficiency and accuracy.

[0088] like Figure 6 As shown, the sum of the areas of the R segmented horizontal planes can be solved by simple arithmetic operations, and the estimated main phase volume fraction C can be obtained based on the ratio of the sum of the areas to the unit grid area. 0 i,jThis method has low computational difficulty, high computational efficiency, and good interface reconstruction effect.

[0089] S6: Define the convergence criterion err for solving the interface equation and compare it with the estimated main phase volume fraction C 0 i,j and the actual main phase volume fraction C i,j , and calculate the difference D between the two:

[0090] D = C 0 i,j - C i,j (2)

[0091] When |D|≤err, execute S9;

[0092] When |D|>err, execute S7;

[0093] Considering the estimated main phase volume fraction C 0 i,j and the actual main phase volume fraction C i,j There is an error, so step S6 needs to be executed to correct the initial interface reconstruction result, and directly use the estimated main phase volume fraction C of the current time layer grid (i, j) calculated in S5 0 i,j The actual main phase volume fraction C of the current layer is preset i,j The comparison is performed to calculate the difference in volume fraction of the main phase, D. The convergence criterion, err, is set according to the specific application scenario. The absolute value of the difference D is compared with err to complete the correction process. The correction method disclosed in this embodiment is simple, rapid, and can be adjusted according to the specific application scenario, with good computational efficiency and compatibility.

[0094] S7: Modify the interface equation. When D>0 in S6, modify the interface equation intercept B=BH; when

[0095] When D≤0 in S5, modify B=B+H;

[0096] Where B is the intercept of the interface equation; H is the adjustment parameter, and the initial value of H is the unit grid characteristic size;

[0097] like Figure 5 As shown, when D>0, it indicates that the estimated main phase volume fraction C 0 i,j Than the preset actual main phase volume fraction C i,j If D is large, it means that the initial value of intercept B is too large. Set B=BH to reduce the value of intercept B and recalculate. When D≤0, it means that the estimated main phase volume fraction C 0 i,jThan the preset actual main phase volume fraction C i,j If the value is small, it means that the initial value of intercept B is too small. Set B=B+H, increase the value of intercept B and recalculate.

[0098] S8: H = H / 2, repeat S4-S8 until |D| ≤ err in S6;

[0099] When |D|≤err in S6, the estimated main phase volume fraction C is obtained. 0 i,j The actual main phase volume fraction C of the current layer is preset i,j are consistent and reach the convergence condition.

[0100] In an optional embodiment, through verification by a large number of experiments, the present invention proposes that the convergence criterion of S6 is err=10e-5; the value of err can better balance the calculation accuracy with the calculation efficiency.

[0101] S9: Obtain the interface equation through S3-S8, and calculate the LS function using the geometric method based on this equation;

[0102] The above calculations yield an accurate intercept for the interface equation, completing the first reconstruction. Since the volume function is essentially a step function, the resulting interface is not precise after the first reconstruction. Given the inherent fineness and smoothness of the LS function, a geometric method is employed to calculate the LS function around the interface, building upon the initial reconstruction and completing the second reconstruction.

[0103] S10: Recalculate the phase interface direction n using the LS function x i,j and n y i,j , then repeat S2-S9 N times to complete multiple interface reconstructions, obtain the final LS function and interface equation, obtain the position and shape of the phase interface, and then describe the effect of the phase interface on flow and heat transfer;

[0104] On the basis of the second interface reconstruction, the obtained LS function is used to recalculate the phase interface direction n x i,j and n y i,j , and then solve the LS function through N iterations (generally N = 3), and so on to complete multiple interface reconstructions, update the LS function and interface equations, and achieve accurate reconstruction of the interface.

[0105] In an optional embodiment, the specific steps of repeating S10 N times S2-S9 to complete multiple interface reconstructions are as follows: set N=3; in the first iteration, use the interface reconstructed by S2-S9 for the first time to solve the LS function, and use the LS function to calculate the phase interface direction n x i,j and n y i,j When , the subsequent iterations use formulas (5) and (6) to determine the phase interface direction. After three iterations, the accurate interface is finally obtained.

[0106]

[0107]

[0108] Where: n x i,j and n y i,j represents the components of the interface direction vector in the x and y directions; φ i,j Represents the size of the LS function of the grid (i, j), Δx and Δy represent the unit grid length in the x and y directions.

[0109] S11: Solve the VOF function using a geometric method, calculate the flux of the main phase volume fraction at the interface of the current grid (i, j), add them together to get the updated main phase volume fraction C, and use C to update the actual main phase volume fraction C of the next time layer i,j , to achieve the description of phase interface migration and obtain the phase distribution of the entire field;

[0110] After obtaining an accurate and smooth reconstruction interface through N iterations of S9, the geometric method is used to advance the interface at the next moment. The calculation diagram is as follows: Figure 7 As shown, Figure 7 (a) represents the upper interface, (b) the lower interface, (c) the left interface, and (d) the right interface. The volume function of the four surfaces, top, bottom, left, and right, is solved for each unit grid, and the true volume fraction of the primary phase in the next time step is calculated.

[0111] The flux of the main phase volume fraction at the interface of the current grid (i, j) in S11 is CH t 、CH b 、CH l 、CH r , and the updated main phase volume fraction C is defined as:

[0112]

[0113]

[0114]

[0115]

[0116] C=C0+CH b -CH t +CH l -CH r (11)

[0117] Where: subscripts t, b, l and r represent the upper interface, lower interface, left interface and right interface of the current grid respectively; v t 、v b 、u l and u r are the velocities of the upper, lower, left, and right interfaces of the current grid; n is the number of rectangular regions divided below the interface; h is the number of rectangular regions divided below the interface; i is the height of the i-th rectangular area of the current grid; Δt is the time step; Δl is the width of each rectangular area; CH t is the flux of the upper interface volume fraction, CH b is the flux of the lower interface volume fraction, CH l is the flux of the left interface volume fraction, CH r is the flux of the volume fraction of the right interface; C0 represents the volume fraction of the previous time layer.

[0118] S12: Repeat S1-S11 until the preset simulation time is reached.

[0119] like Figure 8 As shown, completing steps S1-S11 means completing the solution of the volume function in the current time layer, and then assigning the volume function of the current time layer to the previous time layer. Repeating the above steps until the preset simulation time is reached completes the horizontal piecewise linear reconstruction and tracking of the two-dimensional two-phase interface.

[0120] In summary, the VOSET-HOPLIRE method disclosed in this paper provides a new interface capture method for two-phase flow numerical research. While ensuring computational accuracy and solution efficiency, this method significantly reduces the difficulty of two-phase flow numerical simulation, lowering the threshold for two-phase flow numerical simulation research. This will allow more outstanding scholars to participate in two-phase flow research and promote the application of two-phase flow-related technologies and engineering.

[0121] In order to verify the practical effect of the horizontal piecewise linear reconstruction method of a two-dimensional two-phase interface disclosed in the present invention, a comparative experiment was conducted on the bubble deformation error (Serr) and mass error (Merr) when solving shear bubbles. The comparative experiment used the traditional VOF-PLIC method, the VOSET-PLIC method and the horizontal piecewise linear reconstruction method of a two-dimensional two-phase interface (VOSET-HOPLIRE method) disclosed in the present invention, respectively for the two cases of 100×100 and 200×200 grid sizes. In terms of bubble deformation error, the experimental results are as follows Figure 9 As shown in Figure 2, the VOSET-HOPLIRE method disclosed in the present invention is superior to the VOF-PLIC method and has similar results to the VOSET-PLIC method. Figure 10 As shown, the VOSET-HOPLIRE method disclosed in the present invention can keep the mass error within a reasonable range while adopting a simpler calculation method compared with the VOF-PLIC and VOSET-PLIC methods.

[0122] In terms of accuracy: those skilled in the art usually use pseudo-velocity to evaluate accuracy, and pseudo-velocity is inversely proportional to accuracy. In this example, a comparative experiment was conducted on a stationary bubble using the traditional VOF-PLIC method, the VOSET-PLIC method, and the VOSET-HOPLIRE method disclosed in the present invention. The pseudo-velocity distribution experimental results of the three methods are shown in Figure 2. Figure 11 As shown in the figure, the pseudo velocity around the bubble simulated by the VOF-PLIC method is relatively large, with a maximum value of 4.731×10 -3 m / s; the maximum pseudo-velocity of bubbles simulated using the VOSET-PLIC and VOSET-HOPLIRE methods is 3.233×10 -4 m / s and 3.237×10 -4 m / s, with a difference of only 0.12%. Therefore, the VOSET-HOPLIRE method is more accurate than the VOF-PLIC method, and its accuracy is almost the same as that of the VOSET-PLIC method, which also has high accuracy.

[0123] In terms of conservation: the above three methods are used to analyze the representative moments of the floating bubbles. Figure 12 The bubble mass simulated by the three methods during the bubble rise process does not change significantly. In quantitative analysis, the relative mass deviations of the VOF-PLIC method, the VOSET-PLIC method, and the VOSET-HOPLIRE method compared with the initial bubble mass are 3.82×10 -6 , 1.02×10 -4 and 6.62×10-5 It can be seen that the VOSET-HOPLIRE method disclosed in the present invention is slightly better than the VOSET-PLIC method in terms of conservation, and achieves an effect close to that of the VOF-PLIC method while ensuring accuracy.

[0124] In terms of implementation difficulty, the PLIC method commonly used in the existing technology for interface reconstruction and interface advancement requires a series of tedious coordinate system transformation, rotation, iteration and other operations to determine the intercept α of the line. This method has problems with low computational efficiency and high computational difficulty. It is also difficult to expand the two-dimensional model to a three-dimensional model, which is not conducive to the promotion of numerical simulation methods.

[0125] In the VOSET-HOPLIRE method disclosed in the present invention, the initial interface equation initially passes through the center point of the grid to obtain the initial value of the intercept B, and then the main phase area surrounded by the interface equation and the grid interface of the grid (i, j) is divided into R horizontal areas, and the areas of the R areas are summed to obtain the estimated main phase volume fraction C of the current time layer grid (i, j) 0 i,j , C 0 i,j The actual volume function value C of the current layer is preset i,j Compare and adjust the intercept B until the convergence criterion err is met, and finally determine the intercept of the straight line L, and solve the volume function value C by a simple summation method. i,j Therefore, the VOSET-HOPLIRE method disclosed in this invention has low computational difficulty and high computational efficiency. It can expand two-dimensional models to three-dimensional models, which is conducive to the promotion of numerical simulation methods. While ensuring computational accuracy, it can also take into account computational efficiency and reduce operational difficulty to achieve interface tracking. This method can achieve an accurate description of the phase interface of two-phase flow, significantly reducing the difficulty of two-phase flow numerical research, and providing technical support for two-phase flow numerical research.

[0126] Under the same hardware environment, the VOSET-HOPLIRE method disclosed in this invention took 16.08 hours to compute, while the VOSET-PLIC method took 17.53 hours, significantly reducing computational time. The average error of the VOSET-HOPLIRE method disclosed in this invention remained within 3%, improving computational efficiency by 9.0%, demonstrating the high efficiency and accuracy of the disclosed method for horizontal piecewise linear reconstruction of a two-dimensional two-phase interface.

[0127] The present invention is superior to VOF-PLIC in terms of accuracy, can strictly guarantee the conservation of the research object in terms of conservation, has slightly higher computational efficiency than VOSET-PLIC, is much less difficult to operate than VOSET-PLIC, and can be more conveniently expanded to three dimensions.

Claims

1. A two-dimensional two-phase interface horizontal piecewise linear reconstruction method, characterized in that: The steps include: S1: Grid division of the calculation area; S2: For the grid (i, j) where the two-phase interface is located, the interface direction n is preliminarily calculated based on the volume fraction of the main phase. x i,j and n y i,j , and establish a local coordinate system according to the phase interface direction to ensure that the main phase is below the phase interface, where i and j are the coordinates of the grid; S3: Set the initial interface through the center of the grid (i, j) to obtain the initial interface equation; S4: Use R segmented horizontal surfaces to approximate the phase interface and use the interface equation to calculate the height of the segmented horizontal surface; S5: Calculate the sum of the areas below the segmented horizontal plane, and obtain the estimated main phase volume fraction C based on the ratio of the sum of the areas to the unit grid area. 0 i,j : Where R is the number of segmented horizontal planes, S r is the area below the horizontal plane of a single segment, S g is the unit grid area size; S6: Define the convergence criterion err for solving the interface equation and compare it with the estimated main phase volume fraction C 0 i,j and the actual main phase volume fraction C i,j , and calculate the difference D between the two: D=C 0 i,j -C i,j (2) When |D|≤err, execute S9; When |D|>err, execute S7; S7: Modify the interface equation. When D>0 in S6, modify the interface equation intercept B=BH; when When D≤0 in S5, modify B=B+H; Where B is the intercept of the interface equation; H is the adjustment parameter, and the initial value of H is the unit grid characteristic size; S8: H = H / 2, repeat S4-S8 until |D| ≤ err in S6; S9: Obtain the interface equation through S3-S8, and calculate the LS function using the geometric method based on this equation; S10: Recalculate the phase interface direction n using the LS function x i,j and n y i,j , then repeat S2-S9 N times to complete multiple interface reconstructions, obtain the final LS function and interface equation, obtain the position and shape of the phase interface, and then describe the effect of the phase interface on flow and heat transfer; S11: Solve the VOF function using a geometric method, calculate the flux of the main phase volume fraction at the interface of the current grid (i, j), add them together to get the updated main phase volume fraction C, and use C to update the actual main phase volume fraction C of the next time layer i,j , to achieve the description of phase interface migration and obtain the phase distribution of the entire field; S12: Repeat S1-S11 until the preset simulation time is reached; The specific steps of repeating S10 N times S2-S9 to complete multiple interface reconstructions are as follows: set N=3; use the interface reconstructed by S2-S9 for the first iteration to solve the LS function, and use formulas (5) and (6) in subsequent iterations to determine the phase interface direction, and finally obtain an accurate interface; Where: n x i,j and n y i,j represents the components of the interface direction vector in the x and y directions; φ i,j Represents the size of the LS function of the grid (i, j), Δx and Δy represent the unit grid length in the x and y directions.

2. A two-dimensional two-phase interface horizontal piecewise linear reconstruction method according to claim 1, characterized in that: The direction of the phase interface in S2 and is defined as: Where: C i,j represents the volume fraction of the main phase of the grid (i, j), and Δx and Δy represent the unit grid length in the x and y directions.

3. The horizontal piecewise linear reconstruction method of a two-dimensional two-phase interface according to claim 1, characterized in that: The value range of R in S4 is: 40≥R≥10.

4. The horizontal piecewise linear reconstruction method of a two-dimensional two-phase interface according to claim 1, characterized in that: The convergence criterion for solving the interface equation of S6 is err=10e-5.

5. The horizontal piecewise linear reconstruction method of a two-dimensional two-phase interface according to claim 1, characterized in that: The flux of the main phase volume fraction at the interface of the current grid (i, j) in S11 is CH t 、CH b 、CH l 、CH r , and the updated main phase volume fraction C is defined as: C=C0+CH b -CH t +CH l -CH r (11) Where: subscripts t, b, l and r represent the upper interface, lower interface, left interface and right interface of the current grid respectively; v t 、v b 、u l and u r are the velocities of the upper, lower, left, and right interfaces of the current grid; n is the number of rectangular regions divided below the interface; h is the number of rectangular regions divided below the interface; i is the height of the i-th rectangular area of the current grid; Δt is the time step; Δl is the width of each rectangular area; CH t is the flux of the upper interface volume fraction, CH b is the flux of the lower interface volume fraction, CH l is the flux of the left interface volume fraction, CH r is the flux of the volume fraction of the right interface; C0 represents the volume fraction of the previous time layer.

Images