Two-phase flow interface capturing method for underwater gas jet flow field

Abstract

The invention discloses a two-phase flow interface capturing method of an underwater gas jet flow field, which belongs to the technical field of fluid mechanics, and specifically comprises the following steps: reconstructing an interface in each unit, representing a multiphase flow interface in each grid by using an inclined plane, the normal direction of a unit interface can be calculated through the volume fraction of each unit and the volume fraction of the adjacent unit, the position of the interface can be determined by combining the volume fraction of the unit, and the fluid volume transportation volume of each control surface of the unit can be calculated according to a volume cutting method after the position of the interface is determined. And subtracting the volume net flux of the outflow unit from the main phase fluid volume of the grid unit at the previous moment, and dividing by the total volume of the grid unit to obtain the volume fraction of the unit at the next moment, thereby continuously propelling to obtain the two-phase flow interface of each grid unit.

Description

Technical Field

[0001] The invention relates to the technical field of fluid mechanics, and in particular to a two-phase flow interface capture method for an underwater gas jet flow field. Background Art

[0002] During underwater ignition, when high-pressure gas is ejected from the nozzle, due to the presence of the surrounding water medium, the gas will expand rapidly and form bubbles. This process involves the interaction between the gas phase (gas) and the liquid phase (water), thus constituting a typical two-phase flow problem.

[0003] In this complex flow phenomenon, accurately capturing and describing the behavior of the gas-liquid interface is crucial, but also extremely challenging. This is because with the rapid diffusion of gas and the propagation of pressure waves, the interface position, shape and even the entire flow field structure will change dramatically; in addition, the influence of factors such as surface tension and viscosity effects must be considered, which increases the difficulty of accurately modeling such problems. Therefore, developing technologies or methods that can effectively handle such dynamic interface changes is of great significance for understanding the relevant physical mechanisms and optimizing engineering design. Summary of the invention

[0004] The purpose of the present invention is to provide a two-phase flow interface capture method for an underwater gas jet flow field to solve the following technical problems: How to simulate and capture the change process of two-phase flow interface caused by underwater gas ignition.

[0005] The purpose of the present invention can be achieved through the following technical solutions: A method for capturing a two-phase flow interface in an underwater gas jet flow field comprises the following steps: S1: Embed the fluid interface into the static grid through the scalar value, and define the volume fraction scalar field function f(x, t), where x represents the position variable, t represents the time variable, f(x, t) = 1 represents the primary phase fluid, and f(x, t) = 0 represents the secondary phase fluid; S2: Calculate the volume fraction F of the main phase fluid in each grid, that is, the ratio of the integral of the labeling function f(x, t) in the entire grid to the grid volume; S3: Reconstruct the interface through the PLIC method, determine the interface normal vector, and solve the intersection points of the interface and each edge of the mesh through the geometric method to obtain the interface position; S4: Calculate the volume flux flowing through each control surface of the grid, and update the volume fraction F of each grid at the next moment to advance the flow field in time; S5: Repeat S3 and S4 until the set end time.

[0006] As a further solution of the present invention: in S2, the process of calculating the volume fraction F of the primary phase fluid in each grid, that is, the ratio of the integral of the marker function f(x,t) over the entire grid to the grid volume, is as follows: The calculation formula for the volume fraction F is: ; In the formula, V c represents the volume of a single grid. If F = 0, it means an empty grid and there is no primary phase fluid in the current grid. If F = 1, it means a full grid and the current grid is filled with the primary phase fluid. If 0 < F < 1, it means a half grid and the fluid interface is located in the current grid.

[0007] As a further solution of the present invention: in S3, the process of reconstructing the interface by the PLIC method and determining the interface normal vector is as follows: Reconstruct the interface within each grid, and reconstruct the fluid interface in the grid as a plane. The expression of the intercept form of the interface is: ; In the formula, is the normal vector of the two-phase flow interface, represents the position coordinate vector of the intersection point of the interface and the grid. a is a constant, and the interface normal is determined by calculating the gradient of the volume fraction. The calculation formula is: ; In the formula, represents the gradient of the volume fraction F, and it is set that the interface normal is positive from the inside of the primary phase fluid to the outside of the primary phase fluid.

[0008] As a further solution of the present invention: in S3, the gradient of the volume fraction F is calculated by gradient interpolation based on the MPS method.

[0009] As a further solution of the present invention: in S3, the process of obtaining the interface position by solving the intersection points of the interface and each side of the grid by a geometric method is as follows: The grid is a tetrahedron. Mark the 4 vertices of the tetrahedron as A0, A1, A2, A3. Set a plane P. The normal of the plane P is the same as the interface normal, and the plane P passes through any vertex of the tetrahedron, which is expressed as: ; In the formula, represents the direction vector from any vertex of the grid to another vertex. Obtain the distances from all vertices to the plane P, sort all vertices from near to far according to the distances, mark them as A, B, C, D, and obtain the distances d A d B d C dD ; Set planes B and C with the same normal as the interface through point B and point C respectively, mark the volume ratio of the tetrahedron separated from the mesh by plane B to the whole mesh as fB, and mark the volume ratio of the tetrahedron separated from the mesh by plane C to the whole mesh as fC; For any grid, when the fluid volume fraction F≤fB, the interface is a triangle, and the interface in the grid is marked as the first type of interface unit; when the fluid volume fraction F≥1-fC, the interface is a triangle, and the interface in the grid is marked as the third type of interface unit; when fB<F<1-fC, the interface is a quadrilateral, and the interface in the grid is marked as the second type of interface unit; Get the volume fraction F and distance d for each interface type A d B d C d D , select the corresponding geometric method according to the interface type, and calculate the position coordinate vector of the intersection of the interface and the grid.

[0010] As a further solution of the present invention: In S4, the process of calculating the volume flux flowing through each grid internal interface is: The volume flux through each control surface of the mesh is calculated based on the geometric cutting method. The control surface is the boundary plane of the tetrahedral mesh. For the control surface whose volume flux is to be calculated, a plane μ parallel to the control surface is set. The intersection points of the plane μ and the mesh edge are marked as P1, P2 and P3. The distance between the plane μ and the corresponding control surface is V n Δt, where V n is the normal velocity of the fluid in the direction of the control surface, Δt is the time step at each moment; When the interface is a triangular interface, the control surface is divided into a main control surface that intersects with the fluid interface and a secondary control surface that does not intersect with the fluid interface. For the main control surface, obtain the interface grid intersection point that is not on the main control surface, mark it as point X, and mark the distance from point X to the main control surface as h. X , if h X ≤V n ·Δt, then the volume flux outflow in the time step is all the main phase fluid in the grid, if h X >V n ·Δt, then the volume flux flowing out in the time step is the main phase fluid between the main control surface and plane μ; For the secondary control surface, obtain the interface grid intersection point closest to the secondary control surface, marked as point I; mark the distance from point I to the main control surface as h Y , if h Y ≤V n·Δt, then the volume flow out of the time step is the main phase fluid between the secondary control surface and plane μ. If h Y >V n ·Δt, then the outflow volume flux in the time step is 0; When the interface is a quadrilateral interface, for any control surface, obtain the two interface mesh intersection points that are not on the control surface, mark them as point Z and point W, and mark the distances from point Z and point W to the main control surface as h. Z and h W ; The control surface with a quadrilateral intersection between the interface and the control surface is marked as the main control surface, and the control surface with a triangle intersection between the interface and the control surface is marked as the secondary control surface. For the main control surface, when h Z and h W Both are greater than V n ·Δt, the flux outflowing from the main control surface is the main phase fluid between the main control surface and plane μ; when h Z and h W All are less than V n When Δt, the volume flux flowing out in the time step is all the main phase fluid in the grid; For the secondary control surface, when h Z and h W Both are greater than V n ·Δt, the flux outflowing from the secondary control surface is the main phase fluid between the secondary control surface and plane μ. Z and h W All are less than V n When Δt is equal to 0, the outflow volume flux in the time step is 0.

[0011] As a further solution of the present invention: in S5, the volume fraction F of each unit at the next moment is updated, and the process of time advancement of the flow field is as follows: The volume fraction F of each grid at the kth moment is known k and normal , subtract the net outflow volume flow from the fluid volume at the kth moment to get the volume difference, and divide the volume difference by the total volume of the grid to get the volume fraction F at the k+1th moment k+1 , according to the volume fraction F k+1 Change the flow field properties, advance the flow field in time, and obtain the true conservation variable values ​​and original variable values ​​of the flow field at time k+1.

[0012] Beneficial effects of the present invention: The present invention reconstructs the interface in each unit, and the multiphase flow interface in each grid is represented by an inclined plane. The normal of the unit interface can be calculated through the volume fraction of each unit and its adjacent units, and the position of the interface can be determined in combination with the volume fraction of the unit. After determining the interface position, the fluid volume transport of each control surface of the unit can be calculated according to the volume cutting method, and then the volume of the main phase fluid of the grid unit at the previous moment is subtracted from the volume net flux outflowing the unit and then divided by the total volume of the unit to obtain the volume fraction of the unit at the next moment, thereby continuously advancing the acquisition of the two-phase flow interface of each grid unit, improving the interface capture capability in complex environments such as translational flow fields and shear flow fields, and having good conservation properties. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] The present invention will be further described below in conjunction with the accompanying drawings.

[0014] Figure 1 It is a schematic diagram of the process of the present invention; Figure 2 is a schematic diagram of the unstructured grid gradient model of the present invention; Figure 3 It is a schematic diagram of interpolation of the MPS method of the present invention; Figure 4 It is a schematic diagram of a unit in the original state of the present invention and a unit after the vertex order is changed; Figure 5 is a classification diagram of the half-grid unit of the present invention; Figure 6 It is a schematic diagram of the first type of half-grid flux solution of the present invention; Figure 7 It is a schematic diagram of the second type of half-grid flux solution of the present invention; Figure 8 It is a schematic diagram of the third type of half-grid flux solution of the present invention; Fig. 9 It is the phase diagram at each moment on the diagonal cross section of the hollow sphere of the present invention; Fig.10 is the isosurface diagram of the interface vof at each moment of rotation of the notched disk of the present invention; Fig.11 This is a phase diagram of the initial state of the ball moving in the vortex flow field of the present invention; Fig.12 This is the evolution process of the sphere of the present invention in the vortex flow field. DETAILED DESCRIPTION

[0015] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0016] Please refer to Figure 1-Figure 12 As shown, the present invention is a two-phase flow interface capturing method for an underwater gas jet flow field:

[0017] In the first embodiment, the present invention uses the VOF method as the basis for interface capturing. The VOF method is an interface capturing method based on Eulerian grids, which captures the fluid interface by embedding the fluid interface into a static grid through a scalar value. A volume fraction scalar field function f( x , t) is defined for each fluid, as shown in the following formula: ; If f( x , t) = 1, it represents the primary phase fluid. If f( x , t) = 0, it represents the secondary phase fluid, and the volume fraction F of the primary phase fluid in each cell is calculated, which is the ratio of the integral of the marker function f( x , t) over the entire grid to the grid volume, as shown in the following formula: ; If F = 0, it represents an empty grid, and the grid does not contain the primary phase fluid. If F = 1, it represents a full grid, and the grid is filled with the primary phase fluid. If 0 < F < 1, it is a half grid, and the interface is within this grid.

[0018] In the formula represents the volume of the cell, and the transport equation of the VOF scalar function f( x , t) is as shown in the following formula.

[0019] ; is the velocity vector. By integrating over the grid control volume and using Gauss's theorem, the above formula can be transformed into its integral form, as shown in the following formula.

[0020] ; In the formula is the area of the k-th face of the cell, is the outward normal direction of the k-th face of the cell. By performing time discretization, the following formula can be obtained.

[0021] ; ; Where the normal velocity of the kth surface is . The above formula can be used to transform the problem into solving the net volume flux of the VOF scalar function through each unit control surface. To calculate the volume flux flowing through each unit control surface, it is necessary to first reconstruct the interface position based on the volume fraction and the interface normal vector. After solving the volume flux, the volume fraction of each unit at the next moment can be obtained, and then the flow field properties can be changed according to the volume fraction, and the time of the flow field can be advanced to obtain the actual conservation variable value and the value of the original variable of the flow field at the next moment.

[0022] Determination of interface normal: The present invention adopts the PLIC method to reconstruct the interface in each unit. In the two-dimensional case, a line segment with a slope is used to replace the interface in this unit. In the three-dimensional case, an inclined plane is used to replace the interface in this unit.

[80] The intercept form of the interface can be expressed as follows.

[0023] ; In the formula, is the normal vector of the two-phase flow interface, is the position coordinate vector on the interface, and a is a constant. In the algebraic VOF method, the solution is obtained by iteration, but the number of iterations is large and the amount of calculation is large. The PLIC method adopted in the present invention reconstructs the interface by geometric method, directly solves the equation to obtain the intercept between the interface and each side of the unit, and can save computing resources.

[0024] There are many ways to calculate the normal direction of the interface in the formula. Since the volume fraction value in the main phase fluid is 1 and the volume fraction value in the secondary phase fluid is 0, the normal direction should point to the direction where the volume fraction value changes most dramatically. The present invention directly uses the gradient of the volume fraction to replace the normal direction, as shown in the following formula, and stipulates that the direction from the inside of the main phase fluid to the outside of the main phase fluid is positive.

[0025] ; The unit normal vector can be expressed as follows.

[0026] ; Taking the two-dimensional case as an example, the specific method of calculating the interface normal is introduced. The three-dimensional case is similar, except that there is an additional coordinate in the z direction. Figure 2 A gradient model of an unstructured grid system is given. For the cells in the model, assuming that the volume fraction F can be expanded by Taylor series at cell i, the volume fraction of cell j is It can be written as follows: ; The gradient vector of unit i can be derived from the gradient vector of unit j, which can be written as follows: ; ; The most important thing is how to use the volume fraction interpolation at the adjacent units to obtain the volume fraction gradient of this unit. For structured grids, the nine-point difference method can be used to determine it. For unstructured grids, the more commonly used methods are direct averaging, moving particle semi-implicit method (MPS) and least squares method. The direct averaging method directly averages the volume fractions of the units coplanar with the current unit. Considering the accuracy problem, it is necessary to let as many adjacent units as possible participate in the solution, but this method only considers coplanar units, so the accuracy is low. The MPS method and the least squares method not only take into account the coplanar units but also the co-point units, so the accuracy is higher. Considering that the least squares method is generally used for linear gradient solutions of physical quantities, but the VOF function is a discontinuous function, the present invention uses the MPS method for gradient interpolation.

[0027] The basic idea of ​​the MPS method is to first interpolate the volume fraction of each vertex of this unit using the unit that has points in common with the current unit, and then use the values ​​of these three vertices to calculate the volume fraction gradient. Figure 3 As shown, the formula is as follows: ; ; ; Where r represents the position of the point, n0 represents a scalar function, d represents the dimension of the problem to be solved, and ω represents a weighted interpolation function. is the volume fraction value of the vertex of this unit obtained by interpolation of the common point units, as shown below: ; Therefore, the normal direction in the two-dimensional case can be written as follows: ; The normal direction in three dimensions can be written as follows: ; Determination of interface position: In two-dimensional space, the reconstructed piecewise linear interface is composed of line segments. The present invention uses the two endpoints of the line segment to represent each line segment. The position of the interface can be determined by finding the coordinates of the two endpoints of each line segment. For the three-dimensional case, the linear interface reconstructed in each tetrahedral grid unit is a polygonal plane. The position of the interface can be determined by finding the vertices of the polygon. After obtaining the volume fraction of the unit and the normal vector of the interface, the vertices of the interface can be obtained through a series of geometric operations. The whole process of geometric operations will be described in detail below.

[0028] In order to facilitate the classification of interfaces and simplify subsequent programming calculations in tetrahedral meshes, the four vertices of the unit must first be sorted and located. First, assume that the four vertices of the unit before sorting are A0, A1, A2, and A3, and the vertices of the unit after sorting are represented as A, B, C, and D. Assume a plane P passing through the unit vertices, and the normal of the plane is the normal of the interface , traverse all the vertices of the unit so that all the vertices are on one side of the unit, that is, satisfy the following formula, as shown in the attached Figure 4 Shown on the left.

[0029] ; In the formula, is the normal vector of the interface, is the direction vector from any vertex of the unit to another vertex. After ensuring that the assumed plane P satisfies all vertices on one side, the distance d from all vertices to the assumed plane P is solved. To make each distance satisfy the formula, the vertices of the original unit can be sorted. The unit after the vertex sorting and positioning is shown in the attached figure. Figure 4 As shown on the right, plane P passes through point A of the unit. That is: ; After sorting and locating all the vertices of the unit, the intersection points of the reconstructed interface and the edges of the unit can be calculated through geometric operations, so as to obtain the vertices of the interface and determine the position of the interface to prepare for flux calculation. According to the volume ratio of the main phase fluid inside the unit, the half grid can be divided into three categories, such as Figure 5 As shown, the shaded part in the figure represents the main phase fluid, ensuring that point A is inside the main phase fluid and the interface normal is positive from inside the main phase fluid to outside the main phase fluid.

[0030] With the given normal vector n Draw a plane through vertex B as shown in the attached figure Figure 5 As shown on the left. The intersection points of this plane with AC and AD are F and E respectively. Then the volume ratio of tetrahedron ABEF to ABCD can be written as follows.

[0031] ; Similarly, with the given normal vector Draw a plane through vertex C. The plane intersects DA and DB at L and K, respectively, as shown in the attached figure. Figure 5 As shown on the right, the volume ratio of tetrahedron DCKL to ABCD can be written as follows.

[0032] ; If the fluid volume fraction of this unit , then the linear interface must cut the tetrahedron ABEF, such as Figure 5 As shown on the left. The present invention calls this the first type of interface unit. , then the linear interface must cut the tetrahedron DCKL, such as Figure 5 As shown on the right. The present invention calls it the third type of interface unit. Otherwise, the linear interface must be between two planes passing through B and C, that is, fB<F<1-fC, as Figure 5 As shown in the middle, the present invention refers to it as the second type of interface unit.

[0033] There are also some special cases in the calculation. ,So , points E and F collapse to vertex A. The following formula is also valid for these special cases. The only thing to note is that in order to avoid the denominator being 0, in actual programming, you need to explicitly The situation is classified as the third category and , and The situation is classified as the second category, and .

[0034] For the first type of interface unit, the interface polygon is a triangle, and the coordinates of its vertices G, H, and I can be calculated by the following formula.

[0035] ; For the third type of interface unit, the interface polygon is also a triangle, and the coordinates of its vertices G, H, and I can be calculated by the following formula.

[0036] ; For the second type of interface unit, the situation is more complicated. The interface polygon is an irregular quadrilateral, and the coordinates of its vertices G, H, I and J can be calculated by the following formula.

[0037] ; in: ; ; in: ; The geometric PLIC method is simple to implement. The only parameters required are the volume fraction value F of the unit and the interface normal vector. Although the geometric construction process is relatively complicated, it saves computing resources compared with the iterative method. At present, the method in the present invention only uses the grid center format algorithm, and the grid format needs to be improved. The volume fraction on the grid point can be interpolated to the grid center through the interpolation factor and then calculated.

[0038] Solution for fluid volume flux: To calculate the interface position at each moment, it is necessary to know the fluid volume fraction F of the unit at each moment and the normal direction of the unit interface After the volume fraction at time k is known, the fluid volume flux of each control surface of the unit is first calculated, and then the volume fraction at time k+1 is obtained by subtracting the net outflow volume flux from the fluid volume of the unit at time k and dividing it by the total volume of the unit.

[0039] In the present invention, the flux of each control surface adopts the geometric cutting method, assuming a plane parallel to the control surface to be solved ,flat The intersection points with the unit edge are P1, P2 and P3. The distance from the solution plane is ,in is the normal velocity of the plane to be solved, is the time step, only In the case of the out-of-plane normal direction, the main phase fluid cut by the plane is the volume flux of the unit flowing out through this plane. There are as many as 40 relative positions between the cutting plane and the main phase fluid, but they are generally derived from six cases. Therefore, the present invention selects six typical cases for three types of half grids for illustration. Figure 6 -Attached Figure 8 Schematic diagrams of these six situations are given.

[0040] For the first type of half-grid, when solving the flux of surface ABD, let is the distance from point H to surface ABD. , then the outflow flux is the total main phase fluid ,like , then the flux of outflow surface ABD is Figure 6 The volume of the shaded area in the left figure. The flux calculation method for surface ABC and surface ACD is similar to that for surface ABD. When solving the flux for surface BCD, let is the distance from point I to surface BCD, when When , the flux outflowing from surface ABD is Figure 6 As shown in the shaded area on the right, When , the flux of outflow surface ABD is 0.

[0041] For the second type of half grid, when solving the flux of surface ABD, let and are the distances from point I and point H to surface ABD respectively. and Are greater than When , the flux outflowing from plane ABD is Figure 7 The volume of the shaded area in the left figure. and All less than When the flux outflowing from surface ABD is the entire main phase fluid. The solution method for the flux of surface ABC is similar to that of surface ABD. When solving the flux of surface BCD, let and are the distances from point G and point H to surface ABD respectively. and Are greater than When , the flux outflowing from plane BCD is Figure 7 The volume of the shaded area in the right figure. The solution method for the flux of surface ACD is similar to that of surface ABD.

[0042] For the third type of half grid, Figure 8 The schematic diagram of flux solution for two typical cases of surface ABD and surface ABC is similar to the first and second half-grid methods. There are several other cases of relative relationships, but they are all solved by solving the volume of the cut main phase fluid (irregular polyhedron). The method is the same as above, so I will not go into details here.

[0043] Based on the above method, the fluid transport of all grid cells in the computational domain can be obtained. Then, the volume fraction distribution of each cell can be updated through them, and the interface can be reconstructed in combination with the obtained interface normal. Therefore, the entire process of reconstructing the interface on the unstructured tetrahedral grid cell using the PLIC method can be described as: (1) Before the calculation begins, the volume fraction of each interface unit is initialized according to the initial conditions, with the unit volume fraction of the main phase fluid being 1 and the unit volume fraction of the other phase being 0.

[0044] (2) The local volume fraction gradient is used to solve the interface normal, and then the vertices of the interface are determined by the PLIC method to obtain the cross-sectional position.

[0045] (3) The fluid transport through each control surface of the unit is calculated using the flux calculation method described above, and then the volume fraction in the unit is updated to obtain the fluid volume fraction of each unit at the next moment.

[0046] (4) Repeat steps 2 and 3 until the specified end time.

[0047] Embodiment 2: Four verification examples are used in the present invention to verify the accuracy of the PLIC-VOF method of the present invention in capturing the two-phase flow interface. They are the translation of a hollow sphere, the rotation of a notched cylinder, the deformation of a sphere in a vortex flow field, and the deformation of a sphere in a shear flow field, all of which have achieved good results.

[0048] 1. Translation of a hollow sphere: The computational domain of this example uses a 1×1×1 hexahedral region, and the number of unstructured tetrahedral grids is 1.778 million. Assume that the initial center position of the hollow sphere is (0.2, 0.2, 0.2), the inner diameter is 0.1, and the outer diameter is 0.15. The part between the inner diameter and the outer diameter is defined as the main phase fluid, and the volume fraction is set to 1, and the volume fraction of the rest is 0. The velocity field is given at the initial moment: ; Fig. 9 The phase diagram of the diagonal section of the computational domain at each moment is given. Fig. 9 (a) is the initial state.

[0049] Under the action of the velocity field, the hollow sphere moves along the diagonal line for half a cycle. Fig. 9 As shown in (b), after moving one cycle, the cloud graph is as follows Fig. 9 As shown in (c), it can be seen from the cloud map that no obvious deformation occurs during the movement of the hollow sphere, and the original state can be well maintained. The PLIC-VOF method adopted in the present invention can capture the interface position more accurately in the translational flow field.

[0050] 2. Rotation of notched disc in flow field: The calculation domain used in this example is 1×1×0.3, the calculation grid is 1.082 million, and the center of the circle is assumed to be at (0,0,0.15). The radius of the disk is defined as 0.4, the width of the notch is defined as 0.15, the thickness of the cylinder is defined as 0.1, and the interior of the notched disk is the main phase fluid. Fig.10 The vof isosurface diagram of the interface at each moment is given. Fig.10 (a) is the initial state. At the initial moment, the velocity field of the entire flow field is given as shown in the following equation. Assume that the notched disk rotates at a constant angular velocity and returns to its original position after one rotation. The state of the rotation process is shown in Fig.10 (b) and 10(c).

[0051] ; Where, t x and t yare the x and y coordinates of the center point of the tetrahedral unit, respectively. The notched disk rotates under the action of a constant velocity field, and its shape remains basically unchanged during the rotation process. Since the PLIC-VOF method approximates the two-phase flow interface through a plane, a small cut angle will appear at the open angle interface. After a rotation cycle, the interface returns to its original state. The interface capture algorithm used in the present invention can better maintain the shape of the interface in a rotating flow field.

[0052] 3. The ball moves in the vortex flow field: This is a classic example to verify the interface capture algorithm. , the number of grids is 2.167 million. First, define a spherical region and make its interior the main phase fluid, such as Fig.11 As shown in (a), the coordinates of the center of the sphere are (0.5, 0.3, 0.5), and the radius of the sphere is 0.2. The assumed velocity field is shown below.

[0053] ; Under the violent action of the vortex flow field, the sphere will be stretched and deformed greatly. When the specified time is calculated, the flow field will be reversed with the same cycle. After half a cycle of the sphere's movement, the interface is obviously stretched into a spiral shape, such as Fig.11 As shown in (b), after half a cycle of the reverse velocity field, the interface basically returns to its initial state. Fig.11 As shown in (c), due to the irregular shape of the non-structural tetrahedral unit, the total volume of the main phase fluid is slightly reduced, but this is still within a reasonable range.

[0054] 4. The ball moves in the shear flow field: In order to verify the ability of this method to capture the deformation interface in three-dimensional space, the deformation of the sphere in the three-dimensional velocity field is simulated. The calculation domain of this example adopts a 1×1×1 unit calculation domain, a calculation grid of 2.167 million, and a sphere with a radius of r=0.15 that is deformed. The given velocity field is shown in the following formula.

[0055] ; In this example, the period T = 3. Similar to the motion of a sphere in a vortex flow field, when the sphere is deformed to T / 2, the velocity field reverses, causing the sphere to return to its initial state at t = T. Fig.12 (a)-(i) describe the evolution of the sphere. It can be seen from the cloud map that under the action of the shear flow field, the sphere undergoes a drastic deformation. It is first compressed into a crescent shape and then, over time, the thickness gradually decreases and curls up.

[0056] After the velocity field of the half-cycle flow field reverses direction, the shape of the interface gradually recovers and basically returns to a spherical state after one cycle, which is consistent with the calculation results. The PLIC-VOF algorithm used in the present invention can capture the three-dimensional large deformation interface in the complex flow field with good effect.

[0057] 5. VOF conservation verification: The verification examples for the interface volume flux all have theoretical solutions. In order to verify the accuracy of the VOF algorithm, the present invention introduces the mass conservation deviation E m , volume fraction distribution deviation E r and the geometric deviation E g Three quantities are used for evaluation, and the calculation formula is shown below.

[0058] ; In the above formula is the theoretical solution of the volume fraction in the i-th unit at time t, F i (t) represents the calculated value of the volume fraction in the i-th unit at time t, is the volume fraction value of the i-th unit at the initial moment. V i is the volume of the ith unit. For the translation and rotation cases, the theoretical values ​​of the volume fractions of each unit are known, because the exact position of the rotation and translation interface is known at any point in time, and the interface does not deform during the entire process. In the other two examples, when the time reaches T / 2, the velocity field reverses direction, and the interface returns to its original initial state, so there is also a theoretical solution. Tables 1 to 4 give the error values ​​of the first four examples.

[0059] Table 1 Translation interface error of hollow sphere

[0060] Table 2 Rotation interface error of notched disk

[0061] Table 3 Interface error of a ball moving in a vortex flow field

[0062] Table 4 Interface error of a ball moving in a shear flow field

[0063] It can be seen from the table data that as the grid becomes denser, the calculation error becomes smaller. When the grid is encrypted to a certain extent, the error reduction will become unobvious. The three error values ​​of the present invention are all within a reasonable range, proving that the PLIC-VOF method of the present invention has high interface capture accuracy and good conservation properties.

[0064] The above is a detailed description of an embodiment of the present invention, but the content is only a preferred embodiment of the present invention and cannot be considered to limit the scope of implementation of the present invention. All equivalent changes and improvements made within the scope of the present invention should still fall within the scope of the patent coverage of the present invention.

Claims

1. A method for capturing the two-phase flow interface of an underwater gas jet flow field, characterized in that: The following steps are involved: S1: Embed the fluid interface into the static grid through the scalar value, and define the volume fraction scalar field function f(x, t), where x represents the position variable, t represents the time variable, f(x, t) = 1 represents the primary phase fluid, and f(x, t) = 0 represents the secondary phase fluid; S2: Calculate the volume fraction F of the main phase fluid in each grid, that is, the ratio of the integral of the labeling function f(x, t) in the entire grid to the grid volume; S3: Reconstruct the interface through the PLIC method, determine the interface normal vector, and solve the intersection points of the interface and each edge of the mesh through the geometric method to obtain the interface position; S4: Calculate the volume flux flowing through each control surface of the grid, and update the volume fraction F of each grid at the next moment to advance the flow field in time; S5: Repeat S3 and S4 until the set end time.

2. The method for capturing the two-phase flow interface of an underwater gas jet flow field according to claim 1, characterized in that: In S2, the process of calculating the volume fraction F of the main phase fluid in each grid, that is, the ratio of the integral of the labeling function f(x, t) in the entire grid to the grid volume is: The calculation formula of volume fraction F is: ； where V c represents the volume of a single grid. If F = 0, it represents an empty grid with no primary-phase fluid in the current grid. If F = 1, it represents a full grid filled with the primary-phase fluid. If 0 < F < 1, it represents a half grid with the fluid interface located in the current grid.

3. The method for capturing the two-phase flow interface of an underwater gas jet flow field according to claim 1, characterized in that: In S3, the process of reconstructing the interface by the PLIC method and determining the interface normal vector is as follows: The interface is reconstructed in each grid, and the fluid interface in the grid is reconstructed into a plane. The intercept form of the interface is expressed as: ； In the formula, is the normal vector of the two-phase flow interface, Represents the position coordinate vector of the intersection of the interface and the grid, a is a constant, and the interface normal is determined by calculating the gradient of the volume fraction. The calculation formula is: ； In the formula, represents the gradient of the volume fraction F, and the interface normal is set to be positive when it points from the inside of the main phase fluid to the outside of the main phase fluid.

4. The method for capturing the two-phase flow interface of an underwater gas jet flow field according to claim 1, characterized in that: In S3, the gradient of the volume fraction F is calculated by gradient interpolation based on the MPS method.

5. The method for capturing the two-phase flow interface of an underwater gas jet flow field according to claim 1, characterized in that: In S3, the process of solving the intersection of the interface and each edge of the grid by geometric method to obtain the position of the interface is as follows: The mesh is a tetrahedron, and the four vertices of the tetrahedron are marked as A0, A1, A2, and A3. A plane P is set, and the normal of the plane P is the same as the normal of the interface, and the plane P passes through any vertex of the tetrahedron, which is expressed as: ； In the formula, Represents the direction vector from any vertex of the mesh to another vertex, obtains the distance from all vertices to plane P, sorts all vertices from near to far according to the distance, marks them as A, B, C, D, and obtains the distance d between vertices A, B, C, D and plane P A ,d B ,d C ,d D ; Set planes B and C with the same normal direction as the interface normal through point B and point C respectively, mark the volume ratio of the tetrahedron separated from the mesh by plane B to the whole mesh as fb, and mark the volume ratio of the tetrahedron separated from the mesh by plane C to the whole mesh as fc; For any grid, when the fluid volume fraction F≤fb, the interface is a triangle, and the interface in the grid is marked as the first type of interface unit; when the fluid volume fraction F≥1-fc, the interface is a triangle, and the interface in the grid is marked as the third type of interface unit; when fb<F<1-fc, the interface is a quadrilateral, and the interface in the grid is marked as the second type of interface unit; Get the volume fraction F and distance d for each interface type A ,d B ,d C ,d D , select the corresponding geometric method according to the interface type, and calculate the position coordinate vector of the intersection of the interface and the grid.

6. The method for capturing the two-phase flow interface of an underwater gas jet flow field according to claim 5, characterized in that: In S4, the process of calculating the volume flux flowing through each mesh interface is: The volume flux through each control surface of the mesh is calculated based on the geometric cutting method. The control surface is the boundary plane of the tetrahedral mesh. For the control surface whose volume flux is to be calculated, a plane μ parallel to the control surface is set. The intersection points of the plane μ and the mesh edge are marked as P1, P2 and P3. The distance between the plane μ and the corresponding control surface is V n Δt, where V n is the normal velocity of the fluid in the direction of the control surface, Δt is the time step at each moment; When the interface is a triangular interface, the control surface is divided into a main control surface that intersects with the fluid interface and a secondary control surface that does not intersect with the fluid interface. For the main control surface, obtain the interface grid intersection point that is not on the main control surface, mark it as point X, and mark the distance from point X to the main control surface as h. X , if h X ≤V n ·Δt, then the volume flux outflow in the time step is all the main phase fluid in the grid, if h X >V n ·Δt, then the volume flux flowing out in the time step is the main phase fluid between the main control surface and plane μ; For the secondary control surface, obtain the interface grid intersection point closest to the secondary control surface, marked as point I; mark the distance from point I to the main control surface as h Y , if h Y ≤V n ·Δt, then the volume flow out of the time step is the main phase fluid between the auxiliary control surface and the plane μ. If h Y >V n ·Δt, then the outflow volume flux in the time step is 0; When the interface is a quadrilateral interface, for any control surface, obtain the two interface mesh intersection points that are not on the control surface, mark them as point Z and point W, and mark the distances from point Z and point W to the main control surface as h. Z and h W ; The control surface with a quadrilateral intersection between the interface and the control surface is marked as the main control surface, and the control surface with a triangle intersection between the interface and the control surface is marked as the secondary control surface. For the main control surface, when h Z and h W Both are greater than V n ·Δt, the flux outflowing from the main control surface is the main phase fluid between the main control surface and plane μ; when h Z and h W All are less than V n When Δt, the volume flux flowing out in the time step is all the main phase fluid in the grid; For the secondary control surface, when h Z and h W Both are greater than V n ·Δt, the flux outflowing from the secondary control surface is the main phase fluid between the secondary control surface and plane μ. Z and h W All are less than V n When Δt is equal to 0, the outflow volume flux in the time step is 0.

7. The method for capturing the two-phase flow interface of an underwater gas jet flow field according to claim 6, characterized in that: In S5, the volume fraction F of each unit at the next moment is updated, and the process of time advancement of the flow field is: The volume fraction F of each grid at the kth moment is known k and normal , subtract the net outflow volume flow from the fluid volume at the kth moment to get the volume difference, and divide the volume difference by the total volume of the grid to get the volume fraction F at the k+1th moment k+1 , according to the volume fraction F k+1 Change the flow field properties, advance the flow field in time, and obtain the true conservation variable values ​​and original variable values ​​of the flow field at time k+1.

Images