A two-dimensional gas-liquid interface curve reconstruction method

Abstract

The present invention discloses a method belonging to the technical field of numerical simulation of gas-liquid two-phases, and particularly relates to a two-dimensional gas-liquid interface curve reconstruction method. Under the action of surface tension, the two-dimensional gas-liquid interface appears as a circle, and can be reconstructed using a standard circle: the interface curvature is calculated according to the volume fraction of the main phase or the distance function to obtain the radius of the circle; according to the direction of the interface, a local coordinate system is established to ensure that the main phase is below the interface; a linear reconstruction method is used to obtain the midpoint position of the phase interface; a temporary radius of the circle is set, and the center of the circle is calculated in combination with the midpoint and direction of the phase interface, and the initial circle equation is obtained; a series of segmented horizontal lines are used to approximate the curved phase interface, and the circle equation is used to calculate the height of the horizontal line and the volume fraction of the main phase is counted; the curve equation is iteratively obtained based on the difference between the statistical results and the actual results. High-precision curve reconstruction of the gas-liquid phase interface is achieved through simple iterative steps, providing technical support for high-reduction numerical research on gas-liquid two-phase flow.

Description

Technical Field

[0001] The present invention belongs to the technical field of gas-liquid two-phase numerical simulation, and in particular relates to a two-dimensional gas-liquid interface curve reconstruction method. Background Art

[0002] Two-phase flow is widely present in industrial fields, such as liquid cooling of chips and boiling heat transfer in nuclear pressurized water reactors. Accurate interface capture is crucial for the study of the characteristics and mechanisms of two-phase flow. In the existing technology, different interface capture methods include volume of fluid (VOF), signed distance function (LS) and coupled volume of fluid and level set method (including CVOFLS and VOSET methods). Among them, the VOF method uses a variable to represent the volume fraction of the main phase fluid in the grid unit, which can ensure the conservation of mass of the two-phase flow and is widely used to solve complex two-phase flow problems.

[0003] In the VOF method, the capture of phase interface migration is achieved by solving the pure convection equation of the fluid volume fraction. In the prior art, a high-order format is used to algebraically solve the convection equation. During the iterative solution process, due to numerical pseudo-diffusion, the gas-liquid interface at the nanoscale will become non-physically thickened, occupying multiple grids, resulting in solution errors. To overcome this problem, the prior art solves the VOF equation by reconstructing the two-phase interface using a geometric method: in each time step, a line (two-dimensional problem) or a surface (three-dimensional problem) is first used to approximate the phase interface; then, based on the obtained interface equation, the fluid volume flux at the grid surface is calculated, and the fluid volume fraction is updated to obtain a high-resolution two-phase interface. In the prior art, horizontal or vertical lines parallel to the grid surface are used to reconstruct the two-phase interface. The implementation process of this reconstruction method is relatively simple, but the accuracy is low. To improve accuracy, the classic piecewise linear interface calculation (PLIC) method is used to approximate the real two-phase interface with a straight line, which can greatly improve the accuracy of interface reconstruction. Different algorithms are proposed for this reconstruction method to determine the slope and intercept of the interface linear equation, making this reconstruction method widely used.

[0004] For two-dimensional gas-liquid two-phase flow problems, the interface appears circular due to surface tension. Therefore, the linearly reconstructed interface equation, signed distance function (calculated using geometric methods), interface curvature, surface tension, and volume fraction flux obtained using the classic PLIC method mentioned above all contain errors, which affect the accuracy of the two-phase flow solution. Furthermore, the calculation process involves many special classifications and judgments, and the implementation steps are relatively complex, which is not conducive to the study of the characteristics and mechanisms of two-phase flow. Therefore, it is urgent to construct a two-dimensional gas-liquid interface curve reconstruction method that can achieve high-precision curve reconstruction of the two-dimensional gas-liquid two-phase flow interface through simple iterative steps. Summary of the Invention

[0005] The present invention aims to provide a two-dimensional gas-liquid interface curve reconstruction method, which is characterized in that a standard circle is used to reconstruct the two-dimensional gas-liquid interface, and specifically comprises the following steps:

[0006] S1: Use uniform two-dimensional grid to divide the calculation area;

[0007] S2: For the grid (i, j) where the two-phase interface is located, the volume fraction of the main phase is used Preliminary calculation of the phase interface curvature κ and the circular radius r corresponding to the two-dimensional gas-liquid interface p , the circle radius r p For: r p =1.0 / κ;

[0008] S3: According to the volume fraction of the main phase Preliminary calculation of the phase interface direction n x i,j and n y i,j , where n x i,j is the x-axis component of the interface direction in the grid (i, j) where the phase interface is located, n y i,j is the y-axis component of the interface direction in the grid (i, j) where the phase interface is located; according to the phase interface direction n x i,j and n y i,j Establish a local coordinate system to ensure that the main phase is below the phase interface;

[0009] S4: Use piecewise linear interface calculation method to perform linear reconstruction of the interface and obtain the coordinates of the interface midpoint (x PLIC ,y PLIC );

[0010] S5: Define the circle radius adjustment variable L as: Where: Δx and Δy are the spatial steps of the two-dimensional orthogonal equidivided grid in the x-axis direction and the y-axis direction respectively; the temporary radius of the circle is defined as set up

[0011] S6: According to the coordinates of the midpoint of the interface (x PLIC ,y PLIC ) position, interface direction n x i,j and n y i,j and the temporary radius of the circle Calculate the temporary center position

[0012] S7: The radius of the circle corresponding to the two-dimensional gas-liquid interface in step S2 is r p and the temporary center position in step S6 Substitute the standard circle equation to obtain the equation of the temporary phase interface curve;

[0013] S8: using M segmented horizontal lines to approximate the curved phase interface, and using the temporary phase interface curve equation in step S7 to calculate the height of the segmented horizontal lines;

[0014] S9: Calculate the sum of the areas below the segmented horizontal lines, and calculate the estimated main phase volume fraction based on the ratio of the sum of the areas below the segmented horizontal lines to the unit grid area.

[0015]

[0016] In formula (1): k is the number of the segmented horizontal line, S k is the area under the kth segment horizontal line, S i,j is the geometric area of ​​the grid where the phase interface is located, and M is the number of segmented horizontal lines;

[0017] S10: Calculate mass conservation error

[0018]

[0019] In formula (2): C i,j is the volume fraction of the main phase fluid in the grid (i, j) where the phase interface is located;

[0020] S11: Determine the mass conservation error like but Otherwise judge like Then repeat the operation of L=L / 2 until the condition is satisfied. like Direct calculation Adjust L = L / 2;

[0021] S12: Repeat steps S6-S11 until the Get the final position of the circle center and get the coordinates of the midpoint of the interface (x PLIC ,y PLIC ), phase interface direction n x i,j and n y i,j And the radius r of the circle corresponding to the two-dimensional gas-liquid interface p Obtain the target grid curve equation;

[0022] S13: Based on the target grid curve equation in step S12, a geometric method is used to calculate the signed distance function φ i,j ;

[0023] S14: Using the signed distance function φ in step S13 i,j , recalculate the phase interface direction n x i,j and n y i,j , the curvature of the phase interface κ and the radius r of the circle corresponding to the two-dimensional gas-liquid interface p ;

[0024] S15: Repeat steps S4-S14 N times to complete N times of interface reconstruction, obtain the final phase interface curve equation, and complete the two-dimensional gas-liquid interface curve reconstruction;

[0025] S16: Using the height of the segmented horizontal line corresponding to the final phase interface curve equation in step S15, a geometric method is used to calculate the fluid volume fraction flux of the interface, update the fluid volume fraction, and complete the two-dimensional gas-liquid interface advancement.

[0026] The steps for calculating the phase interface curvature κ in step S2 are:

[0027]

[0028]

[0029]

[0030]

[0031]

[0032]

[0033]

[0034]

[0035]

[0036] Where: is the component of the interface direction in the x direction at the upper right corner of the grid (i, j), is the component of the interface direction in the x direction at the lower right corner of the grid (i, j), is the component of the interface direction in the x direction at the upper left corner of the grid (i, j), is the component of the interface direction in the x direction at the lower left corner of the grid (i, j), is the component of the interface direction in the y direction at the upper right corner of the grid (i, j), is the component of the interface direction in the y direction at the lower right corner of the grid (i, j), is the component of the interface direction in the y direction at the upper left corner of the grid (i, j), is the component of the interface direction in the y direction at the lower left corner of the grid (i, j), C i,j Represents the volume fraction of the main phase fluid of the grid (i, j); when the signed distance function φ is obtained i,j After that, C in formula (3)-formula (10) i,j Replaced by the signed distance function φ i,j .

[0037] In step S3, the phase interface direction n x i,j and n y i,j is defined as:

[0038]

[0039]

[0040] Where C i,j Represents the volume fraction of the main phase fluid of the grid (i, j). When the signed distance function φ is obtained i,j After that, C in Equation (12) and Equation (13) i,j Replaced by the signed distance function φ i,j .

[0041] The temporary center position in step S6 is defined as:

[0042]

[0043]

[0044] Where: x PLIC is the x-axis coordinate of the midpoint of the interface, y PLIC is the y-axis coordinate of the midpoint of the interface.

[0045] The value of M in step S8 is 30.

[0046] The final phase interface curve equation in step S15 is:

[0047]

[0048] In formula (16): (x p ,y p ) is the center of the final phase interface curve equation.

[0049] The specific steps of repeating steps S4-S14 N times to complete N times of interface reconstruction in step S15 are as follows: set N=3; in the first iteration, use the interface reconstructed for the first time in steps S4-S14 to solve the signed distance function φ i,j ; In subsequent iterations, use equations (17) and (18) to determine the phase interface direction n x i,j and n y i,j , use equations (19) to (26) to calculate the gradient of the corner position, substitute the gradient into equation (11) to obtain the phase interface curvature κ, and obtain the circular radius r corresponding to the two-dimensional gas-liquid interface p , and obtain the final phase interface;

[0050]

[0051]

[0052]

[0053]

[0054]

[0055]

[0056]

[0057]

[0058]

[0059]

[0060] Where: φ i,j represents the signed distance function of the grid (i, j), is the component of the interface direction in the x direction at the upper right corner of the grid (i, j), is the component of the interface direction in the x direction at the lower right corner of the grid (i, j), is the component of the interface direction in the x direction at the upper left corner of the grid (i, j), is the component of the interface direction in the x direction at the lower left corner of the grid (i, j), is the component of the interface direction in the y direction at the upper right corner of the grid (i, j), is the component of the interface direction in the y direction at the lower right corner of the grid (i, j), is the component of the interface direction in the y direction at the upper left corner of the grid (i, j), is the component of the interface direction in the y direction at the lower left corner of the grid (i, j).

[0061] The fluid volume fraction flux of the interface in step S16 includes: the flux of the upper interface volume fraction ΔC t , the flux of the lower interface volume fraction ΔC b , the flux of the left interface volume fraction ΔC l and the flux ΔC of the right interface volume fraction r , the step of updating the fluid volume fraction is:

[0062]

[0063]

[0064]

[0065]

[0066] C t+Δt =C t +ΔC b -ΔC t +ΔC l -ΔC r (31)

[0067] 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 、v l and v r are the velocities of the upper, lower, left, and right interfaces of the current grid respectively; M is the number of rectangular areas divided below the interface; h 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; C t+Δt Indicates the volume fraction of the current layer, C t Indicates the volume fraction of the previous time layer.

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

[0069] The present invention discloses a two-dimensional gas-liquid interface curve reconstruction method, namely the Horizontal Refined Piecewise Curve Interface Reconstruction (HOPCIR) method. The HOPCIR method overcomes the shortcomings of the existing classic PLIC method, which suffers from low accuracy and complex implementation steps when reconstructing curved gas-liquid interfaces. By reconstructing the gas-liquid interface using a portion of a circle, the curve equation can be determined using the circle's center and radius. During the reconstruction of the phase interface curve, the interface direction and curvature are first calculated using the fluid fraction, and the radius of the circle is calculated using the curvature. The PLIC method is then used to calculate the center of the phase interface within the grid. The initial circle center position is then calculated using the interface direction and center point, and the initial curve equation is determined. A series of segmented horizontal lines are used to approximate the true curved phase interface, and the curve equation is used to solve for the height of the horizontal lines and the area below them, thereby calculating the volume fraction of the main phase. The calculated main phase volume fraction is compared with the true main phase volume fraction, and the center of the curve equation is iterated based on the difference between the two, completing the first interface reconstruction. Based on the reconstructed curve, a geometric method is used to calculate the signed distance function of the grid points around the interface. This signed distance function is then used to calculate a more accurate interface direction and curvature. The reconstruction process is then repeated three times. Because the segmented horizontal lines used during the reconstruction process approximate the true curved phase interface, the main phase area within the unit grid is decomposed into numerous small rectangles. The area flowing out of these small rectangles can be calculated based on the interface flow velocity, allowing the fluid volume flux on the grid surface to be easily solved and the fluid volume fraction to be updated. Therefore, the present invention discloses a two-dimensional gas-liquid interface curve reconstruction method, which can achieve high-reduction curve reconstruction of the two-dimensional gas-liquid two-phase flow interface and high-precision solution of the VOF equation through simple steps, providing technical support for the study of two-phase problems.

[0070] In order to verify the effectiveness of the two-dimensional gas-liquid interface curve reconstruction method disclosed in the present invention, the HOPCIR method disclosed in the present invention and the classical PLIC method were compared and tested. The test results show that the accuracy of the gas-liquid interface reconstructed by the HOPCIR method in the interface capture test is higher than that of the classical PLIC method. Compared with the classical PLIC method, the HOPCIR method reduces the average error of the signed distance function by 5.6 times and improves the curvature accuracy by 1.5 times. Since a fine horizontal line is used to approximate the curved interface, the VOF flux obtained by the HOPCIR method is more accurate than that obtained by the classical PLIC method. For horizontally translated bubbles, the average deformation error of the HOPCIR method is 2.86 times lower than that of the classical PLIC method. For the single vortex problem of Zalesak disk rotation and time reversal, the HOPCIR method reduces the deformation error by 28% to 140% compared with the classical PLIC method. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] Figure 1 A flow chart of a two-dimensional gas-liquid interface curve reconstruction method disclosed in the present invention;

[0072] Figure 2 Schematic diagram of the process of establishing local coordinates based on the interface direction of the present invention, where (a) to (h) represent working conditions 1 to 8, respectively, and (i) represents the only case obtained after rotation or mirroring operation;

[0073] Figure 3 This is a schematic diagram of curve reconstruction of the present invention;

[0074] Figure 4 The schematic diagram of the calculation of the signed distance function based on curve reconstruction of the present invention, where (a) is the grid (x mci ,y mcj ) is the marked area, and (b) is an enlarged view of the marked area;

[0075] Figure 5 Schematic diagram of VOF flux calculation at the grid interface of the present invention;

[0076] Figure 6 Comparison of signed distance function and curvature error calculated using PLIC and HOPCIR methods, where (a) is the comparison between signed distance function and maximum curvature error, and (b) is the comparison between signed distance function and mean curvature error;

[0077] Figure 7 Comparison of deformation errors of horizontal translation bubble reconstruction using PLIC and HOPCIR methods;

[0078] Figure 8 The reason why the deformation error of the horizontal translation bubble reconstructed by HOPCIR method is smaller than that of PLIC method is analyzed;

[0079] Figure 9 Schematic diagram of Zalesak disk rotation, where (a) is a schematic diagram of circular Zalesak disk rotation, (b) is a schematic diagram of rectangular Zalesak disk rotation, and (c) is a schematic diagram of triangular Zalesak disk rotation;

[0080] Figure 10 Comparison of the Zalesak disk rotation error using the PLIC and HOPCIR methods, where (a) is a comparison of the circular Zalesak disk rotation error, (b) is a comparison of the rectangular Zalesak disk rotation error, and (c) is a comparison of the triangular Zalesak disk rotation error.

[0081] Figure 11 Figure 3 is a diagram of the bubble deformation process under the action of a time-reversed single vortex reconstructed using the PLIC method, where (a) is a schematic diagram of the deformation state at 2 s, (b) is a schematic diagram of the deformation state at 4 s, (c) is a schematic diagram of the deformation state at 6 s, and (d) is a schematic diagram of the deformation state at 8 s.

[0082] Figure 12 Figure 3 is a diagram of the bubble deformation process under the action of a time-reversed single vortex reconstructed using the HOPCIR method, where (a) is a schematic diagram of the deformation state at 2 s, (b) is a schematic diagram of the deformation state at 4 s, (c) is a schematic diagram of the deformation state at 6 s, and (d) is a schematic diagram of the deformation state at 8 s.

[0083] Figure 13 Comparison of the deformation errors of bubbles under the action of a time-reversed single vortex reconstructed using the PLIC and HOPCIR methods. DETAILED DESCRIPTION

[0084] The present invention discloses a two-dimensional gas-liquid interface curve reconstruction method, which uses a standard circle to achieve high-precision reconstruction of the gas-liquid interface; adopts a piecewise approximation curve method to easily solve the curve equation and VOF equation; uses a signed distance function to obtain a high-precision interface curvature (the inverse of the circle radius), and thus combines the VOSET method to perform curve reconstruction. The present invention is further described in detail below with reference to the accompanying drawings. Figure 1 The embodiment of the present invention discloses a method for reconstructing a two-dimensional gas-liquid interface curve, which uses a standard circle to reconstruct the two-dimensional gas-liquid interface, and specifically includes the following steps:

[0085] S1: Use uniform two-dimensional grid to divide the calculation area;

[0086] S2: For the grid (i, j) where the two-phase interface is located, the volume fraction of the main phase is used Preliminary calculation of the phase interface curvature κ and the circular radius r corresponding to the two-dimensional gas-liquid interfacep , the circle radius r p For: r p =1.0 / κ;

[0087] The steps of calculating the phase interface curvature κ are:

[0088]

[0089]

[0090]

[0091]

[0092]

[0093]

[0094]

[0095]

[0096]

[0097] Where: is the component of the interface direction in the x direction at the upper right corner of the grid (i, j), is the component of the interface direction in the x direction at the lower right corner of the grid (i, j), is the component of the interface direction in the x direction at the upper left corner of the grid (i, j), is the component of the interface direction in the x direction at the lower left corner of the grid (i, j), is the component of the interface direction in the y direction at the upper right corner of the grid (i, j), is the component of the interface direction in the y direction at the lower right corner of the grid (i, j), is the component of the interface direction in the y direction at the upper left corner of the grid (i, j), is the component of the interface direction in the y direction at the lower left corner of the grid (i, j), C i,j Represents the volume fraction of the main phase fluid of the grid (i, j); when the signed distance function φ is obtained i,j After that, C in formula (3)-formula (10) i,j Replaced by the signed distance function φ i,j .

[0098] S3: According to the volume fraction of the main phase Preliminary calculation of the phase interface direction n x i,j and ny i,j , where n x i,j is the x-axis component of the interface direction in the grid (i, j) where the phase interface is located, n y i,j is the y-axis component of the interface direction in the grid (i, j) where the phase interface is located; according to the phase interface direction n x i,j and n y i,j Establish a local coordinate system to ensure that the main phase is below the phase interface;

[0099] like Figure 2 As shown, Figure 2 (a) to (h) represent working conditions 1 to 8 respectively, and (i) represents the only situation obtained after rotation or mirroring. x i,j and n y i,j The size and direction of can be divided into 8 cases, corresponding to Figure 2 Conditions 1 to 8 (corresponding to Figure 8 The grid number (i, j) is omitted in the equation. y i,j |≥|n x i,j |,n x i,j ≤0,n y i,j ≤0, by rotation or mirror operation (condition 1: left and right mirror; condition 2: 90° clockwise rotation, n y i,j and n x i,j Interchange of size; Working condition 3: rotate 90° clockwise, then mirror left and right, n y i,j and n x i,j Interchangeable size; Working condition 4: rotate 180° clockwise; Working condition 5: rotate 180° clockwise and then mirror left and right; Working condition 6: rotate 270° clockwise, n y i,j and n x i,j Interchange of size; Working condition 7: rotate 270° clockwise, then mirror left and right, n y i,j and n x i,j (Size interchangeable; Working condition 8: no transformation required) to ensure that the main phase is below the phase interface in the local coordinate system;

[0100] The phase interface direction n xi,j and n y i,j is defined as:

[0101]

[0102]

[0103] Where C i,j Represents the volume fraction of the main phase fluid of the grid (i, j). When the signed distance function φ is obtained i,j After that, C in Equation (12) and Equation (13) i,j Replaced by the signed distance function φ i,j .

[0104] S4: The interface linear reconstruction is performed using the piecewise linear interface calculation method (PLIC) to obtain the coordinates of the interface midpoint (x PLIC ,y PLIC );

[0105] S5: Define the circle radius adjustment variable L as: Where: Δx and Δy are the spatial steps of the two-dimensional orthogonal equidivided grid in the x-axis direction and the y-axis direction respectively; the temporary radius of the circle is defined as set up

[0106] S6: According to the coordinates of the midpoint of the interface (x PLIC ,y PLIC ) position, interface direction n x i,j and n y i,j and the temporary radius of the circle Calculate the temporary center position

[0107] The temporary center position is defined as:

[0108]

[0109]

[0110] Where: x PLIC is the x-axis coordinate of the midpoint of the interface, y PLIC is the y-axis coordinate of the midpoint of the interface.

[0111] S7: The radius of the circle corresponding to the two-dimensional gas-liquid interface in step S2 is r p and the temporary center position in step S6 Substitute the standard circle equation to obtain the equation of the temporary phase interface curve;

[0112] S8: Figure 3 As shown, M segmented horizontal lines are used to approximate the curved phase interface, and the height of the segmented horizontal lines is calculated using the temporary phase interface curve equation in step S7;

[0113] The temporary phase interface curve equation is:

[0114] Where: is the x-axis coordinate of the temporary circle center, The y-axis coordinate of the temporary circle center.

[0115] The value of M in step S8 is 30. Numerical experiments have shown that using M=30 segmented horizontal lines to approximate the curved phase interface can improve the solution efficiency and reduce memory consumption while ensuring the accuracy of interface reconstruction.

[0116] S9: Calculate the sum of the areas below the segmented horizontal lines, and calculate the estimated main phase volume fraction based on the ratio of the sum of the areas below the segmented horizontal lines to the unit grid area.

[0117]

[0118] In formula (1): k is the number of the segmented horizontal line, S k is the area under the kth segment horizontal line, S i,j is the geometric area of ​​the grid where the phase interface is located, M is the number of segmented horizontal lines; S10: Calculation of mass conservation error

[0119]

[0120] In formula (2): C i,j is the volume fraction of the main phase fluid in the grid (i, j) where the phase interface is located; S11: judgment of mass conservation error like but Otherwise judge like Then repeat the operation of L=L / 2 until the condition is satisfied. like Direct calculation Adjust L = L / 2;

[0121] S12: Repeat steps S6-S11 until the Get the final position of the circle center and get the coordinates of the midpoint of the interface (x PLIC ,y PLIC ), phase interface direction n x i,j and n y i,jAnd the radius r of the circle corresponding to the two-dimensional gas-liquid interface p Obtain a target grid curve equation; the target grid curve equation is a curve equation obtained by reconstructing the gas-liquid interface in the current calculation grid.

[0122] S13: Based on the target grid curve equation in step S12, a geometric method is used to calculate the signed distance function φ i,j ;

[0123] The calculation process of the signed distance function based on curve reconstruction is as follows Figure 4 As shown, Figure 4 In the figure, (a) is the grid (x mci ,y mcj ) is the marked area, and (b) is the enlarged view of the marked area. First, mark the interface grid (x mci ,y mcj ) around the 7×7 grid, and then calculate the signed distance function of each marked grid; secondly, the marked grid (x i ,y j ) as an example, loop through the 7×7 grids around it, record the grids where all interfaces are located, and calculate the center of the circle (x p ,y p ) and (x i ,y j ) The intersection of the extension line and the reconstructed curve is determined, and whether the intersection is located in the grid where these interfaces are located; if so, the distance is recorded as φ0; if not, φ0 is set to a large positive number 10 10 Finally, calculate the intersection of the mesh surface of all interfaces and the reconstructed curve, and calculate the marked mesh (x i ,y j ) to each intersection point, compare these distances with the size of φ0, and find the minimum value, which is the minimum value of the marked grid (x i ,y j ) is the signed distance function of .

[0124] S14: Using the signed distance function φ in step S13 i,j , recalculate the phase interface direction n x i,j and n y i,j , the curvature of the phase interface κ and the radius r of the circle corresponding to the two-dimensional gas-liquid interface p ;

[0125] S15: Repeat steps S4-S14 N times to complete N times of interface reconstruction, obtain the final phase interface curve equation, and complete the two-dimensional gas-liquid interface curve reconstruction;

[0126] The final phase interface curve equation is:

[0127]

[0128] In formula (16): (x p ,y p ) is the center of the final phase interface curve equation.

[0129] The specific steps of repeating steps S4-S14 N times to complete N times of interface reconstruction in step S15 are as follows: set N=3; in the first iteration, use the interface reconstructed for the first time in steps S4-S14 to solve the signed distance function φ i,j ; In subsequent iterations, use equations (17) and (18) to determine the phase interface direction n x i,j and n y i,j , use equations (19) to (26) to calculate the gradient of the corner position, substitute the gradient into equation (11) to obtain the phase interface curvature κ, and obtain the circular radius r corresponding to the two-dimensional gas-liquid interface p , and obtain the final phase interface;

[0130]

[0131]

[0132]

[0133]

[0134]

[0135]

[0136]

[0137]

[0138]

[0139]

[0140] Where: φ i,j represents the signed distance function of the grid (i, j), is the component of the interface direction in the x direction at the upper right corner of the grid (i, j), is the component of the interface direction in the x direction at the lower right corner of the grid (i, j), is the component of the interface direction in the x direction at the upper left corner of the grid (i, j), is the component of the interface direction in the x direction at the lower left corner of the grid (i, j), is the component of the interface direction in the y direction at the upper right corner of the grid (i, j), is the component of the interface direction in the y direction at the lower right corner of the grid (i, j), is the component of the interface direction in the y direction at the upper left corner of the grid (i, j), is the component of the interface direction in the y direction at the lower left corner of the grid (i, j).

[0141] S16: Using the height of the segmented horizontal line corresponding to the final phase interface curve equation in step S15, a geometric method is used to calculate the fluid volume fraction flux of the interface, update the fluid volume fraction, complete the two-dimensional gas-liquid interface advancement, and then describe the influence of the phase interface on flow and heat transfer.

[0142] like Figure 5 As shown, the fluid volume fraction flux of the interface in step S16 includes: the flux of the upper interface volume fraction ΔC t , the flux of the lower interface volume fraction ΔC b , the flux of the left interface volume fraction ΔC l and the flux ΔC of the right interface volume fraction r , the step of updating the fluid volume fraction is:

[0143]

[0144]

[0145]

[0146]

[0147] C t+Δt =C t +ΔC b -ΔC t +ΔC l -ΔC r (31)

[0148] 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 、v l and v r are the velocities of the upper, lower, left, and right interfaces of the current grid respectively; M is the number of rectangular areas divided below the interface; h 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; C t+Δt Indicates the volume fraction of the current layer, C t Indicates the volume fraction of the previous time layer.

[0149] In summary, the HOPCIR method disclosed in this paper provides a simple and accurate method for reconstructing gas-liquid interface curves for numerical studies of two-phase flows. Through simple iterations, high-precision reconstruction of curved interfaces can be achieved, which is beneficial for studying the characteristics and mechanisms of two-phase flows.

[0150] In order to verify the practical effect of the two-dimensional gas-liquid interface curve reconstruction method disclosed in this invention, based on the VOSET method, the classic PLIC method is used as the object, and the advantages of the HOPCIR method are explained by comparing the distance function and curvature accuracy of stationary bubbles, as well as the deformation errors of horizontal bubble movement, Zalesak disk rotation and time-reversed single vortex problems.

[0151] Stationary bubble: Taking a stationary bubble with a radius of R = 0.25m as the object, the signed distance function is calculated based on a grid of R / Δx = 10 to 50. Figure 6 As shown, (a) is a comparison of the maximum signed distance function and curvature error calculated by the PLIC and HOPCIR methods, and (b) is a comparison of the mean signed distance function and curvature error calculated by the PLIC and HOPCIR methods. Compared with the PLIC method, the average error and maximum error of the signed distance function calculated by the HOPCIR method are smaller, with the maximum error reduced by 15.5 times and the average error reduced by 5.6 times. Based on the more accurate signed distance function, the accuracy of the interface curvature calculated by the HOPCIR method is also higher than that of the PLIC method; the errors of the maximum curvature and average curvature calculated by the latter are 1.8 times and 1.5 times the errors calculated by the former, respectively. It can be seen that the curve reconstruction strategy adopted by the HOPCIR method can achieve higher-precision reconstruction of the interface and obtain more accurate signed distance function and curvature.

[0152] Horizontally moving bubble: A bubble with a radius of 0.2m is initially located at (0.3, 0.5). Under the action of a horizontal velocity of 1.0m / s, the bubble moves horizontally to (0.7, 0.5). The test grid range includes R / Δx = 5 to 50. Figure 7 As shown in Figure 2, the deformation error decreases with the increase of grid resolution; when R / Δx≥10, the deformation error of the reconstruction result of the HOPCIR method is significantly lower than that of the reconstruction result of the PLIC method. Quantitatively speaking, the average deformation error of the PLIC method is 2.86 times that of the HOPCIR method. Figure 8 Explain why HOPCIR reconstruction is more accurate. In the PLIC method, using straight lines to reconstruct curved gas-liquid interfaces introduces some errors in the VOF flux calculation on the mesh surface. In contrast, using the HOPCIR method for curved interface reconstruction helps reduce VOF flux errors, and the local errors offset each other, further reducing the error.

[0153] Zalesak disk rotation: This problem involves placing two-dimensional circles, squares, and triangles in a counterclockwise flow field with a certain regularity. These figures will rotate once and then return to their original positions while maintaining their shape. Figure 9 As shown in Figure 1, (a) is a schematic diagram of the rotation of a circular Zalesak disk, (b) is a schematic diagram of the rotation of a rectangular Zalesak disk, and (c) is a schematic diagram of the rotation of a triangular Zalesak disk. The calculation grids include 25×25, 50×50, 100×100, 150×150, and 200×200. The results are shown in Figure 1. Figure 10 As shown in the figure, (a) compares the rotation errors of a circular Zalesak disk, (b) compares the rotation errors of a rectangular Zalesak disk, and (c) compares the rotation errors of a triangular Zalesak disk. For all studied figures under different grid conditions, the deformation error of the PLIC method is consistently greater than that of the HOPCIR method. For the circle, the deformation error is larger when the number of grids is small; when the grid size is larger than 100×100, the deformation error converges to 0.0023. The HOPCIR method, on the other hand, maintains a consistently low deformation error across different grids. When the grid size is 25×25, the deformation error of the HOPCIR method is 0.0043, which is 1.0 / 6.7 of the PLIC method (0.029). When the grid size increases to 100-100, this ratio drops to 1.0 / 1.1. Across different grids, the average deformation error of the PLIC method is 2.40 times that of the HOPCIR method. For rectangular and triangular rotation reconstruction, the advantage of the HOPCIR method over the PLIC method does not change much as the grid resolution increases; under different grid configurations, the average deformation error of the PLIC method is 1.28 and 1.43 times that of the HOPCIR method, respectively.

[0154] Time-reversed single vortex problem: This problem involves placing a two-dimensional bubble in a time-varying vortex. The bubble is sheared by the clockwise vortex in the first half of the time and by the counterclockwise vortex in the second half of the time, eventually returning to its original position while maintaining its shape. The computational grids include 50×50, 100×100, 150×150, and 200×200. Using a 100×100 grid, the deformation process of the PLIC method under the action of the vortex is as follows: Figure 11 As shown in the figure, (a) is a schematic diagram of the deformation state at 2s, (b) is a schematic diagram of the deformation state at 4s, (c) is a schematic diagram of the deformation state at 6s, and (d) is a schematic diagram of the deformation state at 8s. The deformation process of the HOPCIR method under the action of eddy current is as follows Figure 12As shown in the figure, (a) is a schematic diagram of the deformation state at 2s, (b) is a schematic diagram of the deformation state at 4s, (c) is a schematic diagram of the deformation state at 6s, and (d) is a schematic diagram of the deformation state at 8s. Compared with the HOPCIR method, the bubbles reconstructed by the PLIC method have more broken small droplets at the tail during the shearing process. At the end of the simulation, there is a tiny bubble separated from the restored large bubble in the capture result of the PLIC method, and the circle deviation obtained by the PLIC method is more serious than that of the HOPCIR method. As the grid increases to 200-200, the deformation errors of both methods decrease, especially the HOPCIR method. As Figure 13 As shown in Figure 2, the deformation error ratio between the PLIC and HOPCIR methods increases from 1.52 to 3.68 as the mesh size increases from 50-50 to 200-200. On a 100-100 mesh, the HOPCIR method achieves a deformation error of 0.043, which is comparable to the PLIC method's deformation error of 0.040 at a 150-150 mesh resolution. This indicates that the HOPCIR method requires fewer meshes to accurately capture complex deformation problems.

[0155] From the above comparison results, it can be seen that the HOPCIR method proposed in this paper is superior to the classic PLIC method in terms of distance function, curvature calculation and interface capture accuracy, and can provide technical support for high-precision numerical research on gas-liquid two-phase flow.

Claims

1. A two-dimensional gas-liquid interface curve reconstruction method, characterized in that: The two-dimensional gas-liquid interface reconstruction is performed using a standard circle, which specifically includes the following steps: S1: Use uniform two-dimensional grid to divide the calculation area; S2: For the grid (i, j) where the two-phase interface is located, the volume fraction of the main phase is used Preliminary calculation of the phase interface curvature κ and the circular radius r corresponding to the two-dimensional gas-liquid interface p , the circle radius r p For: r p =1.0 / κ; S3: According to the volume fraction of the main phase Preliminary calculation of the phase interface direction n x i,j and n y i,j , where n x i,j is the x-axis component of the interface direction in the grid (i, j) where the phase interface is located, n y i,j is the y-axis component of the interface direction in the grid (i, j) where the phase interface is located; according to the phase interface direction n x i,j and n y i,j Establish a local coordinate system to ensure that the main phase is below the phase interface; S4: Use piecewise linear interface calculation method to perform linear reconstruction of the interface and obtain the coordinates of the interface midpoint (x PLIC ,y PLIC ); S5: Define the circle radius adjustment variable L as: Where: Δx and Δy are the spatial steps of the two-dimensional orthogonal equidivided grid in the x-axis direction and the y-axis direction respectively; the temporary radius of the circle is defined as set up S6: According to the coordinates of the midpoint of the interface (x PLIC ,y PLIC ) position, interface direction n x i,j and n y i,j and the temporary radius r of the circle c Calculate the temporary center position S7: The radius of the circle corresponding to the two-dimensional gas-liquid interface in step S2 is r p and the temporary center position in step S6 Substitute the standard circle equation to obtain the equation of the temporary phase interface curve; S8: using M segmented horizontal lines to approximate the curved phase interface, and using the temporary phase interface curve equation in step S7, calculating the height of the segmented horizontal lines; S9: Calculate the sum of the areas below the segmented horizontal lines, and calculate the estimated main phase volume fraction based on the ratio of the sum of the areas below the segmented horizontal lines to the unit grid area. In formula (1): k is the number of the segmented horizontal line, S k is the area under the kth segment horizontal line, S i,j is the geometric area of ​​the grid where the phase interface is located, and M is the number of segmented horizontal lines; S10: Calculate mass conservation error In formula (2): C i,j is the volume fraction of the main phase fluid in the grid (i, j) where the phase interface is located; S11: Determine the mass conservation error like but Otherwise judge like Then repeat the operation of L=L / 2 until the condition is satisfied. like Direct calculation Adjust L = L / 2; S12: Repeat steps S6-S11 until the Get the final position of the circle center and get the coordinates of the midpoint of the interface (x PLIC ,y PLIC ), phase interface direction n x i,j and n y i,j And the radius r of the circle corresponding to the two-dimensional gas-liquid interface p Obtain the target grid curve equation; S13: Based on the target grid curve equation in step S12, a geometric method is used to calculate the signed distance function φ i,j ; S14: Using the signed distance function φ in step S13 i,j , recalculate the phase interface direction n x i,j and n y i,j , the curvature of the phase interface κ and the radius r of the circle corresponding to the two-dimensional gas-liquid interface p ; S15: Repeat steps S4-S14 N times to complete N times of interface reconstruction, obtain the final phase interface curve equation, and complete the two-dimensional gas-liquid interface curve reconstruction; S16: Using the height of the segmented horizontal line corresponding to the final phase interface curve equation in step S15, a geometric method is used to calculate the fluid volume fraction flux of the interface, update the fluid volume fraction, and complete the two-dimensional gas-liquid interface advancement.

2. The two-dimensional gas-liquid interface curve reconstruction method according to claim 1, characterized in that: The steps for calculating the phase interface curvature κ in step S2 are: Where: is the component of the interface direction in the x direction at the upper right corner of the grid (i, j), is the component of the interface direction in the x direction at the lower right corner of the grid (i, j), is the component of the interface direction in the x direction at the upper left corner of the grid (i, j), is the component of the interface direction in the x direction at the lower left corner of the grid (i, j), is the component of the interface direction in the y direction at the upper right corner of the grid (i, j), is the component of the interface direction in the y direction at the lower right corner of the grid (i, j), is the component of the interface direction in the y direction at the upper left corner of the grid (i, j), is the component of the interface direction in the y direction at the lower left corner of the grid (i, j), C i,j Represents the volume fraction of the main phase fluid in grid (i, j); When the signed distance function φ is obtained i,j After that, C in formula (3)-formula (10) i,j Replaced by the signed distance function φ i,j .

3. The two-dimensional gas-liquid interface curve reconstruction method according to claim 1, characterized in that: In step S3, the phase interface direction n x i,j and n y i,j is defined as: Where C i,j Represents the volume fraction of the main phase fluid of the grid (i, j). When the signed distance function φ is obtained i,j After that, C in Equation (12) and Equation (13) i,j Replaced by the signed distance function φ i,j .

4. The two-dimensional gas-liquid interface curve reconstruction method according to claim 1, characterized in that: The temporary center position in step S6 is defined as: Where: x PLIC is the x-axis coordinate of the midpoint of the interface, y PLIC is the y-axis coordinate of the midpoint of the interface.

5. The two-dimensional gas-liquid interface curve reconstruction method according to claim 1, characterized in that: The value of M in step S8 is 30.

6. The two-dimensional gas-liquid interface curve reconstruction method according to claim 1, characterized in that: The final phase interface curve equation in step S15 is: In formula (16): (x p ,y p ) is the center of the final phase interface curve equation.

7. The two-dimensional gas-liquid interface curve reconstruction method according to claim 1, characterized in that: The specific steps of repeating steps S4-S14 N times to complete N times of interface reconstruction in step S15 are as follows: set N=3; in the first iteration, use the interface reconstructed for the first time in steps S4-S14 to solve the signed distance function φ i,j ; In subsequent iterations, use equations (17) and (18) to determine the phase interface direction n x i,j and n y i,j , use equations (19) to (26) to calculate the gradient of the corner position, substitute the gradient into equation (11) to obtain the phase interface curvature κ, and obtain the circular radius r corresponding to the two-dimensional gas-liquid interface p , and obtain the final phase interface; Where: φ i,j represents the signed distance function of the grid (i, j), is the component of the interface direction in the x direction at the upper right corner of the grid (i, j), is the component of the interface direction in the x direction at the lower right corner of the grid (i, j), is the component of the interface direction in the x direction at the upper left corner of the grid (i, j), is the component of the interface direction in the x direction at the lower left corner of the grid (i, j), is the component of the interface direction in the y direction at the upper right corner of the grid (i, j), is the component of the interface direction in the y direction at the lower right corner of the grid (i, j), is the component of the interface direction in the y direction at the upper left corner of the grid (i, j), is the component of the interface direction in the y direction at the lower left corner of the grid (i, j).

8. The two-dimensional gas-liquid interface curve reconstruction method according to claim 1, characterized in that: The fluid volume fraction flux of the interface in step S16 includes: the flux of the upper interface volume fraction ΔC t , the flux of the lower interface volume fraction ΔC b , the flux of the left interface volume fraction ΔC l and the flux ΔC of the right interface volume fraction r , the step of updating the fluid volume fraction is: C t+Δt =C t +ΔC b -ΔC t +ΔC l -ΔC r (31) 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 、v l and v r are the velocities of the upper, lower, left, and right interfaces of the current grid respectively; M is the number of rectangular areas divided below the interface; h 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; C t+Δt Indicates the volume fraction of the current layer, C t Indicates the volume fraction of the previous time layer.

Images