CN116680762A - Reinitialization geometric cutting correction method for free surface capture - Google Patents

Abstract

The invention provides a method for correcting reinitialized geometric cutting of free surface capture, which comprises the following steps: 1) Dividing a calculation domain into two-dimensional/three-dimensional finite element surface/volume grids; 2) Splitting a finite element grid in a calculation domain into ridge units, and storing ridge unit information; 3) Solving a convection control equation in a Level Set method to obtain LS distance function values of the current time steps on the nodes of each ridge line unit; 4) Retrieving and identifying interface ridge lines; 5) Calculating and storing initial interpolation coefficients, and judging the original position of an interface; 6) Performing an outer loop process: solving a reinitialization equation to obtain a node current distance function value, and 7) executing an inner loop process: correcting the position of the disturbed interface by using the initial geometric range and the form of the interface; 8) Judging whether the internal circulation meets the calculation requirement and executing corresponding operation; 9) Checking whether the virtual iteration step meets the calculation requirement and executing corresponding operation; 10 Repeating steps 3 to 9 until all time steps are calculated.

Description

Reinitialization geometric cutting correction method for free surface capture

Technical Field

The invention relates to a numerical calculation method for interface capture, in particular to a reinitialization geometric cutting correction method based on a geometric correction strategy and a free surface capture technology.

Background

Numerical modeling of two-phase flow mobile interfaces is a critical issue in both engineering and computational science. Among the various numerical methods, the Level Set (LS) method proposed by Osher and Sethian in literature l (Osher S., sethian J.A. front propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations [ J ]. Journal of Computational Physics,1988, 79 (1): 12-49.) makes LS method an important tool for capturing mobile interfaces because of its ability to handle topology merging and fragmentation problems of mobile interfaces.

The LS method is to represent the interface using a set of zero values for the distance function. In order to determine the distribution of physical properties, the LS function needs to be a sign distance function. Sussman proposes a standard reinitialization method in document 2 (Sussman M., smereka P., osher S.Alevel set approach for computing solutions to incompressible two-phase flow [ J ]. Journal of Computational Physics,1994, 114 (1): 146-159.) whereby any scalar field can be converted to a sign distance function by iteratively solving the Hamilton-Jacobi equation to steady state, and the LS method of implicit reinitialization proposed by Sussman is the standard LS method. However, the standard LS method has a major disadvantage in that lack of mass conservation causes excessive mass loss/increase in the calculation process, possibly leading to disturbance of the zero level set interface, and possibly even divergence and collapse of the numerical calculation. This mass non-conservation results mainly from the discretization and re-initialization process of the numerical equations.

Many researchers have made various attempts to overcome the mass non-conservation problem in the standard LS method. In addition to using higher order discrete schemes to compensate for the discretization error of the equation, one straightforward way is to encrypt the grid, but this is an inefficient and costly way. Another method is to couple the LS method with other numerical methods such as CLSVOF method (see document 3 (Liu f., xu y., li Y.A coupled level-set and volume-of-fluid method for simulating axi-symmetric incompressible two-phase flows [ J ]. Applied Mathematics Computation,2017, 293:112-130 ])) and ICLS method (see document 4 (Mao J., zhao l., liu x., et al a wire-phases model for the simulation of landslide-generated waves using the improved conservative level set method [ J ]. Computers & Fluids,2017, 159:243-253 ]). These two methods take advantage of the mass conservation properties of the fluid Volume (VOF) and Conservative LS (CLS) methods, respectively, to avoid the problem of standard re-initialization, but they lose the simplicity of the standard LS method, making its computational cost generally higher than that of the standard LS method, especially in large-scale numerical simulations.

On the other hand, researchers have improved on iterative re-initialization procedures in view of the simplicity and accuracy of interface descriptions. Based on the strategy, there are two main methods, one is to modify the iterative reinitialization equation, and the other is to directly modify the disturbance interface. In contrast to the former, the latter can directly limit interface perturbations for the purpose of capturing accurate interface contours, as the Particle Level Set (PLS) method proposed by Enright in document 5 (Enright D., fedkiw R., ferziger J., et al, A hybrid particle level set method for improved interface capturing [ J ]. Journal of Computational Physics,2002, 183 (1): 83-116.). Although the method has good conservation of mass, the calculation efficiency is obviously reduced due to the introduction of excessive Lagrange particles, and the problems of complex interface correction format, complex particle regeneration process and the like exist.

Disclosure of Invention

Aiming at the defects in the prior art, the problem of high difficulty in accurately describing a free surface of a two-phase flow and the problem of mass non-conservation in the process of re-initializing an LS method cannot be simultaneously met to provide accurate interface description, ensure high calculation efficiency and other conditions, the invention provides a method for correcting the re-initialized geometric cutting by capturing the free surface, overcomes the problem of mass non-conservation in the process of re-initializing the level set method, solves the problem of describing a moving interface in a two-dimensional or three-dimensional space, avoids the problem of mass loss caused by interface disturbance in the process of re-initializing the standard LS method through a simple and efficient correction strategy, and can be applied to accurately describing the moving interface in the two-dimensional or three-dimensional space.

In order to achieve the above purpose, the present invention adopts the following technical scheme:

a method of reinitializing geometric cut correction for free-face capture, comprising the steps of:

step 1, dividing a calculation domain into a two-dimensional finite element surface grid or a three-dimensional finite element body grid;

step 2, dividing the finite element grid in the calculation domain into ridge units, and storing ridge unit information;

step 3, solving a convection control equation in a Level Set method to obtain LS distance function values of the current time steps on the two end nodes of each ridge unit;

step 4, checking whether LS distance function values on two end nodes of the ridge line unit are different in sign, and if yes, identifying the ridge line unit as an interface ridge line unit;

step 5, calculating initial interpolation coefficients of two end points of the interface ridge line unit to obtain coordinates of an interface interpolation point, and further obtaining an initial position of an interface;

step 6, executing an outer circulation process: solving a reinitialization equation in a Level Set method to obtain current LS distance function values of two end nodes of an interface ridge unit;

step 7, executing an inner loop process: updating LS distance function values of two end nodes of the interface ridge unit by using the initial geometric range and the form of the interface, and correcting the disturbance interface position;

step 8, repeating the process of iterating the inner loop, judging whether the inner loop reaches a convergence condition, if so, updating the LS distance function value of the node and entering step 9, otherwise, adding one in the inner iteration step, and returning to step 7;

step 9, checking whether the virtual iteration step meets the maximum iteration times, if so, entering a step 10, otherwise, repeating the steps 6 to 8;

and step 10, repeating the steps 3 to 9 until all time steps are calculated.

In order to optimize the technical scheme, the specific measures adopted further comprise:

further, in step 3, the convection control equation is specifically:

wherein ,represents LS distance function, t represents calculation time, u _f Representing the fluid velocity vector.

Further, the step 5 specifically comprises:

calculating initial interpolation coefficients of two end nodes of interface ridge line unitThe specific formula is as follows:

in the formula ,an initial interpolation coefficient representing a first node of the interface ridge element,/>An initial interpolation coefficient representing a second node of the interface ridge element,/>LS distance function value representing the initial virtual time step of the first node at re-initialization, +.>LS distance function values representing the initial virtual time step of the second node at the reinitialization;

obtaining coordinates of an interface interpolation point P based on linear interpolation

wherein ,X₁ and X₂ The coordinate vectors of the first node and the second node are respectively obtained, and then the initial position of the interface is obtained.

Further, the step 6 specifically includes:

solving a reinitialization equation in the Level Set method:

where τ represents the virtual time, the superscript n represents the virtual time step,representing a smooth sign function, defined as:

wherein Deltax represents the interface width parameter, and the LS distance function value of the ith endpoint in the kth virtual time step after solving the reinitialization equation in the current time step is obtained by the above method

Further, the step 7 specifically includes:

step 7.1, determining a new LS distance function value of the original interface interpolation point PObtained by linear interpolation:

wherein j is the number of inner iteration steps, and j is more than or equal to 0;for representing interface interpolation point displacement errors;

if it isThe P is still located at the position of the original interface interpolation point, the original interface is not disturbed, the step 7 is exited, and the step 9 is directly executed; if->The interface interpolation point generates disturbance, and the step 7.2 is executed;

step 7.2, according to the initial interpolation coefficient, the displacement error is calculatedThe new LS distance function value is redistributed to the nodes at the two ends of the interface ridge to eliminate the interface displacement error when the interface interpolation points are disturbed, the new LS distance function value of the nodes at the two ends of the interface ridge is obtained, the interface position is corrected, and the formula is expressed as follows:

further, the step 8 specifically comprises:

repeating steps 7.1 to 7.2 until the following convergence conditions are reached:or j is larger than or equal to eta, lambda depends on the grid size and the required interface precision, and the value range is 10- ³ ～10- ² The method comprises the steps of carrying out a first treatment on the surface of the η is the number of inner iterations specified, not exceeding 20;

when one of the above convergence conditions is satisfied, ending the iterative loop process in step 7, and updating the LS distance function value of the ith node in the kth virtual time stepAnd step 9 is performed; otherwise, step 7 is repeated, and the j+1th inner iteration of the kth virtual time step is performed.

Further, in step 9, the maximum number of iterations is not less than 2 and not more than 5.

The beneficial effects of the invention are as follows:

since the splitting of the cell into ridges only needs to be performed once during the whole calculation. Moreover, the number of interface ridges is small compared to the entire computational grid. In addition, the interface correction process does not need a plurality of inner iteration steps, and is only implemented on the nodes of the interface ridge, so that the interface correction process is not limited by the shape of the free surface set, and the calculation on the nodes has the advantage of high calculation speed. Therefore, the method has simplicity and high efficiency, remarkably improves the mass conservation of the standard level set method, and can effectively meet the requirement of accurately describing the unstable free surface in the two-phase flow numerical value calculation process.

Drawings

FIG. 1 is a schematic diagram of a finite element mesh split into ridge cells;

FIG. 2 is a schematic diagram of identifying interface ridge cells and interface interpolation points;

FIG. 3 is a schematic diagram of interface perturbation during a reinitialization process;

FIG. 4 is a schematic diagram of disturbance interface correction;

FIG. 5 is a graph of initial moment interface position for the problem of deformed circles in a two-dimensional shear flow field;

FIG. 6 is a graph of the problem of deformed circles in a two-dimensional shear flow field, and the evolution process of the interface profile;

fig. 7 is a graph comparing interfacial profiles at typical moments for the problem of deformed circles in a two-dimensional shear flow field.

Detailed Description

The invention will now be described in further detail with reference to the accompanying drawings.

The invention provides a free surface capturing reinitialization geometric cutting correction method, which comprises the following steps:

step 1, dividing a calculation domain into a two-dimensional finite element surface grid or a three-dimensional finite element body grid;

step 2, dividing the finite element grid in the calculation domain into ridge units, and storing ridge unit information; taking a quadrilateral mesh as an example, splitting the quadrilateral mesh unit in fig. 1 can obtain 4 ridge line units E ₁₂ 、E ₂₃ 、E ₃₄ and E₄₁ 。

Step 3, solving a convection control equation in a Level Set method to obtain LS distance function values of the current time step (nth time step) on the two end nodes of each ridge line unit; the convection control equation is specifically:

wherein ,represents LS distance function, t represents calculation time, u _f Representing the fluid velocity vector.

Step 4, checking whether LS distance function values on two end nodes of the ridge line unit are different in sign, and if yes, identifying the ridge line unit as an interface ridge line unit; taking a quadrilateral grid as an example, after the step 3 is executed, the nodes at the two ends of each ridge line can be obtainedValues, as shown in FIG. 2-> and />Is positive sign, and-> and />Negative sign, then edge line L ₁₂ and L₃₄ Then it is an interfacial ridge unit, and P ₁ and P₂ Respectively are interfacial ridge lines L ₁₂ and L₃₄ Is defined.

Step 5, calculating initial interpolation coefficients of two end points of the interface ridge line unit to obtain coordinates of an interface interpolation point, and further obtaining an initial position of an interface; the method comprises the following steps:

calculating initial interpolation coefficients of two end nodes of interface ridge line unitThe specific formula is as follows:

in the formula ,an initial interpolation coefficient representing a first node of the interface ridge element,/>An initial interpolation coefficient representing a second node of the interface ridge element,/>LS distance function value representing the first node at time step 0,/I>An LS distance function value representing the second node at time step 0;

obtaining coordinates of an interface interpolation point P based on linear interpolation

wherein ,X₁ and X₂ The coordinate vectors of the first node and the second node are respectively obtained, and then the initial position of the interface is obtained.

Step 6, executing an outer circulation process: solving a reinitialization equation in a Level Set method to obtain current LS distance function values of two end nodes of an interface ridge unit; the method comprises the following steps:

solving a reinitialization equation in the Levelset method:

where τ represents the virtual time, the superscript n represents the virtual time step,representing a smooth sign function, defined as:

wherein Deltax represents the interface width parameter, and the LS distance function value of the ith endpoint in the kth virtual time step after solving the reinitialization equation in the current time step is obtained by the above method

As shown in fig. 3, the coordinates of the original interface interpolation point P: wherein X₁ and X₂ The coordinate vectors of nodes 1 and 2, respectively. Interpolation Point in reinitialization->Is represented by linear interpolation as: />However, due to the presence of a numerical dissipation, it is generally +.>And indicating that the interpolation point of the interface generates displacement and the interface is disturbed. Therefore, step 7 needs to be performed to tend to return the perturbed interpolation point to its original position.

This step makes the Level Set function satisfy the properties of the symbolic distance function, but the interface is disturbed during the re-initialization process, so the following steps are needed to be executed to drive the disturbed interface back to its original position.

Step 7, executing an inner loop process: updating LS distance function values of two end nodes of the interface ridge unit by using the initial geometric range and the form of the interface, and correcting the disturbance interface position;

step 7.1, determining a new LS distance function value of the original interface interpolation point PObtained by linear interpolation:

wherein j is the number of inner iteration steps, and j is more than or equal to 0;for representing interface interpolation point displacement errors;

if it isThe P is still located at the position of the original interface interpolation point, the original interface is not disturbed, the step 7 is exited, and the step 9 is directly executed; if->The interface obtained by calculation through the Level Set method is inaccurate, disturbance occurs to the interpolation point of the interface, and the step 7.2 is executed;

step 7.2, according to the initial interpolation coefficient, the displacement error is calculatedThe new LS distance function values of the nodes at the two ends of the interface ridge are obtained by being redistributed to the nodes at the two ends of the interface ridge, the interface displacement error during disturbance of the interpolation point of the interface is eliminated, the interface position is corrected, and the new LS distance function values are expressed as follows by a formula:

taking the quadrilateral unit in fig. 4 as an example, the interface disturbance error when the interface interpolation point moves is redistributed to the nodes at two ends of the interface ridge, so that the interface interpolation point is driven to return to the original position, and the disturbed interface can be restored to the original position.

Step 8, repeating the process of iterating the inner loop, judging whether the inner loop reaches a convergence condition, if so, updating the LS distance function value of the node and entering step 9, otherwise, adding one in the inner iteration step, and returning to step 7;

the convergence condition is specifically:or j.gtoreq.eta. Where λ depends on the mesh size and the required interface accuracy. The value range is 10 ^-3 ～10 ^-2 . And eta is the appointed inner iteration times, and good results can be obtained by generally not exceeding 20 according to numerical example experience.

When one of the above convergence conditions is satisfied, ending the iterative loop process in step 7, and updating the LS function value of the ith node in the kth virtual time stepAnd step 9 is performed; otherwise, step 7 is repeated, and the j+1th inner iteration of the kth virtual time step is performed.

Step 9, checking whether the virtual iteration step meets the maximum iteration number, generally taking 2 times, and generally not exceeding 5 times. If yes, entering a step 10, otherwise, repeating the steps 6 to 8;

and step 10, repeating the steps 3 to 9 until all time steps are calculated.

Examples:

the deformation round calculation example in the two-dimensional shear flow field is a classical calculation example for verifying the capability of capturing an interface by a numerical simulation method, and the calculation example can intuitively reflect the capability of capturing severe topological changes and tiny interfaces in the variable flow field processed by the numerical method.

In this example, an initial circular interface centered on (0.50, 0.75) with a radius of 0.15 is located at [0,1 ]]×[0,1]Square computational domain, as shown in fig. 5. All variables in this example are dimensionless basis. In the method of the invention, for the initial moment, the distance function of each point on the circular interface satisfies the following conditionThe circles deform in the following shear flow fields:

where T is the shear period, taking t=8. The circular interface is stretched into a vortex at this velocity field, where the vortex reaches maximum deformation at time t=t/2=4 and returns to the original interface at t=t=8.

By adopting the method of the invention, the calculation domain in the calculation example is divided into a quadrilateral grid of 400 multiplied by 400, and the number of internal iterations is set to 20. Fig. 6 shows the interface contour evolution process of the deformed circular interface obtained by the method of the present invention over time on a 400×400 grid in a variation flow. Fig. 7 is a graph in which the interface profiles at the time of maximum deformation at t=4 and the time of end of calculation at t=8 are counted, respectively, in which the black dotted line represents an accurate solution and the red solid line represents the numerical calculation result obtained by the method of the present invention.

The calculation result shows that the method can accurately capture the fine interfaces of the vortex tail parts at different moments in the change flow field, and the fine interfaces are basically restored to the initial shape and the initial position at the simulation ending moment, the interfaces at the calculation ending moment have no obvious interface disturbance and deformation, and the problem of capturing the two-phase flow moving interfaces can be effectively solved.

The above is only a preferred embodiment of the present invention, and the protection scope of the present invention is not limited to the above examples, and all technical solutions belonging to the concept of the present invention belong to the protection scope of the present invention. It should be noted that modifications and adaptations to the invention without departing from the principles thereof are intended to be within the scope of the invention as set forth in the following claims.

Claims (7)

1. A method for reinitializing geometric cut correction for free-face capture, comprising the steps of:

step 1, dividing a calculation domain into a two-dimensional finite element surface grid or a three-dimensional finite element body grid;

step 2, dividing the finite element grid in the calculation domain into ridge units, and storing ridge unit information;

step 3, solving a convection control equation in a level set method to obtain LS distance function values of the current time steps on the two end nodes of each ridge unit;

step 4, checking whether LS distance function values on two end nodes of the ridge line unit are different in sign, and if yes, identifying the ridge line unit as an interface ridge line unit;

step 5, calculating initial interpolation coefficients of two end points of the interface ridge line unit to obtain coordinates of an interface interpolation point, and further obtaining an initial position of an interface;

step 6, executing an outer circulation process: solving a reinitialization equation in a level set method to obtain current LS distance function values of two end nodes of an interface ridge unit;

step 7, executing an inner loop process: updating LS distance function values of two end nodes of the interface ridge unit by using the initial geometric range and the form of the interface, and correcting the disturbance interface position;

step 8, repeating the process of iterating the inner loop, judging whether the inner loop reaches a convergence condition, if so, updating the LS distance function value of the node and entering step 9, otherwise, adding one in the inner iteration step, and returning to step 7;

step 9, checking whether the virtual iteration step meets the maximum iteration times, if so, entering a step 10, otherwise, repeating the steps 6 to 8;

and step 10, repeating the steps 3 to 9 until all time steps are calculated.

2. The method of claim 1, wherein in step 3, the convection control equation is specifically:

wherein ,represents LS distance function, t represents calculation time, u _f Representing the fluid velocity vector.

3. The method of claim 1, wherein step 5 is specifically:

calculating initial interpolation coefficients of two end nodes of interface ridge line unitThe specific formula is as follows:

in the formula ,an initial interpolation coefficient representing a first node of the interface ridge element,/>An initial interpolation coefficient representing a second node of the interface ridge element,/>LS distance function value representing the initial virtual time step of the first node at re-initialization, +.>LS distance function values representing the initial virtual time step of the second node at the reinitialization;

obtaining coordinates of an interface interpolation point P based on linear interpolation

wherein ,X₁ and X₂ The coordinate vectors of the first node and the second node are respectively obtained, and then the initial position of the interface is obtained.

4. The method of claim 1, wherein step 6 is specifically:

solving a reinitialization equation in the Levelset method:

where τ represents the virtual time, the superscript n represents the virtual time step,representing a smooth sign function, defined as:

wherein Deltax represents the interface width parameter, and the LS distance function value of the ith endpoint in the kth virtual time step after solving the reinitialization equation in the current time step is obtained by the above method

5. The method of claim 1, wherein step 7 specifically comprises:

step 7.1, determining a new LS distance function value of the original interface interpolation point PObtained by linear interpolation:

wherein j is the number of inner iteration steps, and j is more than or equal to 0;for representing interface interpolation point displacement errors;

if it isThe P is still located at the position of the original interface interpolation point, the original interface is not disturbed, the step 7 is exited, and the step 9 is directly executed; if->The interface interpolation point generates disturbance, and the step 7.2 is executed;

step 7.2, according to the initial interpolation coefficient, the displacement error is calculatedThe new LS distance function value is redistributed to the nodes at the two ends of the interface ridge to eliminate the interface displacement error when the interface interpolation points are disturbed, the new LS distance function value of the nodes at the two ends of the interface ridge is obtained, the interface position is corrected, and the formula is expressed as follows:

6. the method of claim 5, wherein step 8 is specifically:

repeating steps 7.1 to 7.2 until reachingThe following convergence conditions:or j is larger than or equal to eta, lambda depends on the grid size and the required interface precision, and the value range is 10 ^-3 ～10 ^-2 The method comprises the steps of carrying out a first treatment on the surface of the η is the number of inner iterations specified, not exceeding 20;

when one of the above convergence conditions is satisfied, ending the iterative loop process in step 7, and updating the LS distance function value of the ith node in the kth virtual time stepAnd step 9 is performed; otherwise, step 7 is repeated, and the j+1th inner iteration of the kth virtual time step is performed.

7. The method of claim 1, wherein in step 9, the maximum number of iterations is not less than 2 and not more than 5.