Novel two-phase flow interface capturing method

Abstract

The invention discloses a novel two-phase flow interface capturing method, which relates to the field of two-phase flow interfaces, and comprises the following steps: selecting a smooth characteristic function of an area occupied by a fluid as an indicator function of two fluids, and initializing; on the basis of an initialization result, utilizing reverse tracking along a feature line and meshless function interpolation to update a smooth feature function after the interface position is changed; and re-initializing the updated smooth characteristic function to complete the capture of the two-phase flow interface. According to the method, the problems that numerical values are difficult to solve when an irregular calculation area is solved and the interface capture precision becomes poor during long-time simulation in an existing method are solved.

Description

Technical Field

[0001] The invention relates to the field of two-phase flow interface, in particular to a novel two-phase flow interface capturing method. Background Art

[0002] Two-phase flow interface problems are widely found in chemical production, petroleum industry, metallurgy and mining, thermal energy and energy engineering, etc. In chemical production, the study of two-phase flow interface problems can better understand the flow phenomena in the reactor, help optimize operating conditions, and improve production efficiency and product quality; in the petroleum industry, the study of two-phase flow interface problems can optimize oil and gas transportation and reduce energy consumption; in heat exchangers and cooling systems, the application of two-phase flow interface problems can improve heat exchange efficiency. Therefore, efficient numerical simulation methods for two-phase flow interface problems have broad application prospects and can promote the development of industrial production technology.

[0003] When two fluids with different physical properties and are incompatible come into contact, an interface appears between them. Because fluids cannot maintain their shape under any or a certain shear force, the interface between the fluids moves and deforms during the flow. The dynamic behavior of the interface is often complex and depends not only on the properties of the fluids but also on the flow conditions.

[0004] In order to successfully simulate flow interface problems, the first task is to accurately track or capture the moving and deforming interface during the fluid flow process. In the past few decades, dozens of numerical techniques have been proposed to track / capture the interface. From the perspective of the coordinate system, these methods can be roughly divided into two groups: Lagrangian methods (variable computational domain) and Eulerian methods (fixed computational domain).

[0005] However, existing methods often have numerical difficulties in solving irregular calculation areas, and the interface capture accuracy deteriorates during long simulation times. Summary of the invention

[0006] In view of the above-mentioned deficiencies in the prior art, the present invention provides a novel two-phase flow interface capture method that solves the numerical difficulties of the prior art when solving irregular calculation areas and the problem that the interface capture accuracy deteriorates during long-term simulation.

[0007] In order to achieve the above-mentioned purpose, the technical solution adopted by the present invention is: a new two-phase flow interface capture method, comprising the following steps:

[0008] S1: Select the smooth characteristic function of the area occupied by one fluid as the indicator function of the two fluids and initialize it;

[0009] S2: Based on the initialization results, the smooth characteristic function is updated after the interface position changes by using reverse tracing along the characteristic line and meshless function interpolation;

[0010] S3: Reinitialize the updated smooth characteristic function to complete the capture of the two-phase flow interface.

[0011] Furthermore, the initialization formula in S1 is:

[0012]

[0013] Among them, χ(·) is the smooth characteristic function, x is the position of the fluid micro-group, φ(·) is the signed distance function, and η is the fluid smoothness intensity.

[0014] Furthermore, the reverse tracking along the characteristic line in S2 includes the following sub-steps:

[0015] A1: Calculate the new position of the fluid micro-group along the characteristic line of the flow field;

[0016] A2: Based on the new position reached by the fluid micro-group and the displacement along the characteristic path, the starting point position of the fluid micro-group is obtained;

[0017] A3: Update the smooth characteristic function based on the relationship between the starting position of the fluid cluster and the new position reached along the characteristic line of the flow field.

[0018] Furthermore, the new position reached by the fluid micro-group in A1 along the characteristic line of the flow field is:

[0019] x a =x d +v(x m ,t+Δt / 2)Δt

[0020] Among them, x a is the new position reached by the fluid cluster along the characteristic line of the flow field, x d is the starting point of the fluid microgroup, v(·) is the flow velocity, x m is the midpoint of the characteristic path, t is the time, and Δt is the time interval;

[0021] Under the assumption of uniform motion, the midpoint x of the characteristic path m for:

[0022]

[0023] Under the assumption of uniform accelerated motion, the midpoint x of the characteristic path m for:

[0024]

[0025] a=(v(xa ,t+Δt)-v(x d ,t)) / Δt

[0026] Where a is the estimated value of acceleration.

[0027] Furthermore, the starting point position of the fluid micro-group in A2 is:

[0028] The displacement d of the fluid particle along the characteristic path is:

[0029] d=v(x a -d / 2,t+Δt / 2)Δt

[0030] The fixed point iteration is:

[0031] d [k+1] =v(x a -d [k] / 2,t+Δt / 2)Δt

[0032] Among them, d [k+1] is the displacement after the k+1th iteration, d [k] is the displacement after the kth iteration;

[0033] The starting point of the fluid cluster is:

[0034] x d =x a -d.

[0035] Furthermore, the formula for updating the smooth characteristic function in A3 is:

[0036] χ(x a ,t+Δt)=χ(x d ,t)

[0037] Among them, χ(·) is the smooth characteristic function.

[0038] Furthermore, the gridless function interpolation in S2 includes the following sub-steps:

[0039] B1: Assume that at x d There are n in the support domain of x randomly distributed gridless nodes {x i}, the characteristic function value at this node is Get Taylor polynomial approximation;

[0040] B2: Use the weighted least squares method to solve the Taylor polynomial approximation and obtain the smooth characteristic function at the starting point x of the fluid micro-group. d The value of χ(x d ,t), and through χ(x d ,t)Update the smooth characteristic function.

[0041] Furthermore, the Taylor polynomial approximation in B1 is:

[0042]

[0043] α=(α1,…,α n )

[0044] in, is the Taylor polynomial approximation, α is the multi-index vector in the n-dimensional case, is the function value in the Taylor polynomial approximation, are the derivative values ​​in the Taylor polynomial approximation, α1, α2, …, α n is the multi-index component in the n-dimensional case, i represents the i-th gridless node;

[0045] definition:

[0046]

[0047] Then the Taylor polynomial approximation is:

[0048]

[0049] Where c(·) and p(·) are set vectors, and the superscript T represents the transpose of the matrix.

[0050] Furthermore, the Taylor polynomial approximation solved by weighted least squares method in B2 is:

[0051]

[0052] Among them, J(·) is the objective function of the least squares problem, min represents the minimum value, and w i (·) is the weight function associated with the node, P is the basis function matrix, W(·) is the weight function matrix, is the characteristic function node value vector, Φ is the shape function matrix, p(·) is the basis function vector, i=1,…,n x ;

[0053] set up is a unit coordinate vector, whose kth component is equal to 1, using e k The vector c(x d ) is expressed analytically as Thus:

[0054]

[0055] Φ=(P T W(x d )P) -1P T W(x d )

[0056] The first row of Φ corresponds to the function value The kth row corresponds to the derivative value

[0057] Among them, α k is the kth multi-index vector in ascending order, M is the total number of components, and m is the order of the Taylor polynomial;

[0058] Then the smooth characteristic function is at the starting point x of the fluid micro-group d The value of χ(x d ,t) is:

[0059]

[0060] Among them, e1 is the first component of the unit coordinate vector.

[0061] Furthermore, in S3, the updated smooth characteristic function is reinitialized, and the formula is:

[0062]

[0063] in, is the reconstructed signed distance function, η is the fluid smoothness intensity, is the characteristic function of dissipation, l is the distance between any point x and the interface, is the gradient.

[0064] The beneficial effects of the present invention are: 1) The new method does not require explicit solution of partial differential equations, and is simpler and more practical. 2) The new method is a completely meshless method, which inherits the advantage of meshless methods without mesh restrictions. 3) The new method has natural conservation properties and has the characteristics of low dissipation and high resolution. 4) The new method can be fully parallel, has high computational efficiency, and can be suitable for large-scale simulation. 5) The new method is suitable for complex geometric areas and irregular node distributions, and exhibits excellent performance in solving complex interface problems. 6) For some extremely challenging public test examples, the new method can give better results. 7) The overall performance of the new method is better than most existing methods. BRIEF DESCRIPTION OF THE DRAWINGS

[0065] Figure 1 This is a flow chart of a new two-phase flow interface capture method.

[0066] Figure 2 Schematic diagram of a two-phase flow problem.

[0067] Figure 3 This is a test chart for the reinitialization effect.

[0068] Figure 4 This is the simulation result diagram before reinitialization.

[0069] Figure 5 This is the simulation result diagram after reinitialization.

[0070] Figure 6 The following is a comparison chart of the simulation before and after reinitialization.

[0071] Figure 7 Schematic diagram of the notched disk problem in Example 1.

[0072] Figure 8 This is a comparison chart of the simulation results of the notched disk problem when it rotates one cycle.

[0073] Fig. 9 The figure is a comparison chart of the simulation results of the notched disk problem under uniform 257×257 points when rotating 2, 3, and 4 cycles.

[0074] Fig.10 Schematic diagram of the single vortex shear problem in Example 2.

[0075] Fig.11 This is a comparison chart of simulation results for single vortex shear problem when the period is 2π, 3π, and 4π.

[0076] Fig.12 This is a comparison chart of the simulation results of the three-dimensional notched spherical droplet problem when rotating for one cycle.

[0077] Fig.13 This is a comparison chart of simulation results for the three-dimensional single vortex shear problem when the period is 3. DETAILED DESCRIPTION

[0078] The present invention will be further described below in conjunction with the accompanying drawings and specific embodiments.

[0079] like Figure 1 As shown, a novel two-phase flow interface capture method comprises the following steps:

[0080] S1: Select the smooth characteristic function of the area occupied by one fluid as the indicator function of the two fluids and initialize it;

[0081] S2: Based on the initialization results, the smooth characteristic function is updated after the interface position changes by using reverse tracing along the characteristic line and meshless function interpolation;

[0082] S3: Reinitialize the updated smooth characteristic function to complete the capture of the two-phase flow interface.

[0083] The basic idea of ​​the present invention is to select a smooth characteristic function of the area occupied by a fluid as the indicator function of the two fluids, update the characteristic function after the interface position changes by reverse tracing along the characteristic line and gridless function interpolation, and use a reinitialization technique to maintain the properties of the characteristic function during the flow. Since for a given fluid group, the fluid type remains unchanged during the fluid flow process, the value of the characteristic function is unchanged along the characteristic path / line. This is the core foundation of the method, giving the method a clear physical meaning or explanation. For a new time step, reverse trace along the corresponding characteristic path to find its starting point, and interpolate the function value to the starting point to obtain the characteristic function value at each calculation point. In addition, in order to reduce the interpolation difficulties encountered by the grid-based method, we use gridless function interpolation to obtain the characteristic function value of the starting point, thereby making this method a completely gridless method.

[0084] like Figure 2 As shown, consider a computational region Ω filled with fluid 1 (x∈Ω1) and fluid 2 (x∈Ω2). Assume that the two fluids are physically immiscible and there is a clear interface Γ between them. Under a given flow field, the two fluids will move and deform. In order to track the movement and deformation of the two fluids, the standard characteristic function is taken as:

[0085]

[0086] The above formula is used as the indicator function of fluid 1, and 1-χ(x, t) is used as the indicator function of fluid 2.

[0087] When the two-phase interface is directly identified on a Cartesian grid, the result is often very rough and will appear jagged. In order to overcome this difficulty, the present invention proposes a smooth characteristic function to reduce the dependence of the interface identification result on the grid. When initializing the characteristic function, the signed distance function φ(x) widely used in the level set method is introduced as an auxiliary variable, the smooth intensity η (interface thickness) is set, and the initial value of the smooth characteristic function is defined.

[0088] The initialization formula in S1 is:

[0089]

[0090] Among them, χ(·) is the smooth characteristic function, x is the position of the fluid micro-group, φ(·) is the signed distance function, and η is the fluid smoothness intensity.

[0091] Usually, η is a small positive constant, which depends on the spatial step size h of the grid, and can be taken as η = h or η = 2h. Compared with the standard characteristic function, the smooth characteristic function only modifies the function value in a narrow band area with a width of 2η near the interface. When η approaches zero, the smooth characteristic function converges to the standard characteristic function.

[0092] The reverse tracking along the characteristic line in S2 includes the following sub-steps:

[0093] A1: Calculate the new position of the fluid micro-group along the characteristic line of the flow field;

[0094] A2: Based on the new position reached by the fluid micro-group and the displacement along the characteristic path, the starting point position of the fluid micro-group is obtained;

[0095] A3: Update the smooth characteristic function based on the relationship between the starting position of the fluid cluster and the new position reached along the characteristic line of the flow field.

[0096] In the two-phase interface flow problem, from time t to t+Δt, a d The fluid particle at the point ∈Ω(t) will leave its position and arrive at a new position x along the characteristic line of the flow field. a ∈Ω(t+Δt). Therefore, the time evolution of the characteristic function can be expressed as:

[0097] χ(x a ,t+Δt)=χ(x d ,t)

[0098] To update the characteristic function according to the above formula, we must deduce x d and x a Because these two points are on the same feature line, we have:

[0099]

[0100] Generally speaking, the flow velocity v(x, t) along the characteristic line is related to time, so the integral on the right side of the equation is difficult to calculate analytically. In order to obtain second-order accuracy, the present invention adopts the midpoint quadrature formula when calculating the integral on the right side of the equation. The quadrature formula assumes that the velocity within the entire time interval [t, t+Δt] is constant and takes the value as the velocity at the intermediate moment.

[0101] The new position reached by the fluid micro-group in A1 along the characteristic line of the flow field is:

[0102] x a =x d +v(x m ,t+Δt / 2)Δt

[0103] Among them, x ais the new position reached by the fluid cluster along the characteristic line of the flow field, x d is the starting point of the fluid microgroup, v(·) is the flow velocity, x m is the midpoint of the characteristic path, t is the time, and Δt is the time interval;

[0104] Under the assumption of uniform motion, the midpoint x of the characteristic path m for:

[0105]

[0106] Under the assumption of uniform accelerated motion, the midpoint x of the characteristic path m for:

[0107]

[0108] a=(v(x a ,t+Δt)-v(x d ,t)) / Δt

[0109] Where a is the estimated value of acceleration.

[0110] Assume that the position x of the fluid particle at time t+Δt is a With grid point x i (arrival point) is consistent. At time t before a time step, each fluid particle is located at its respective starting point x d That is, the fluid particle changes from x to d After a time step Δt, it arrives at the grid point x along its characteristic path a Under the assumption of uniform motion, if the displacement along the characteristic path is defined as d = x a -x d ,but:

[0111] The starting point position of the fluid microgroup in A2 is:

[0112] The displacement d of the fluid particle along the characteristic path is:

[0113] d=v(x a -d / 2,t+Δt / 2)Δt

[0114] The fixed point iteration is:

[0115] d [k+1] =v(x a -d [k] / 2,t+Δt / 2)Δt

[0116] Among them, d [k+1] is the displacement after the k+1th iteration, d [k] is the displacement after the kth iteration;

[0117] When setting d = 0 as the initial condition, usually only 2 to 3 iterations are sufficient to obtain a sufficiently accurate result. Once d is obtained, it can be obtained by d =x a -d to get the starting point, so that by the formula χ(x a ,t+Δt)=χ(x d ,t) Update the characteristic function. For the case of uniformly accelerated motion, the calculation is very similar and will not be repeated here.

[0118] The starting point of the fluid cluster is:

[0119] x d =x a -d.

[0120] The formula for updating the smooth characteristic function in A3 is:

[0121] χ(x a ,t+Δt)=χ(x d ,t)

[0122] Among them, χ(·) is the smooth characteristic function.

[0123] Usually starting point x d With grid point x i At time t, the characteristic function is only at the grid point x i There is a value at , so when using χ(x a ,t+Δt)=χ(x d ,t) When updating the characteristic function, you must first use interpolation to get the starting point x d The function value at x. There is no special restriction on the choice of interpolation method. However, for the grid-based method, for each starting point x d , it is necessary to traverse all grid cells before interpolation to determine which grid cell contains this point, and this process is very difficult to implement. In order to reduce the difficulty of interpolation caused by the grid, the present invention uses a gridless interpolation method. The implementation process of this method is given below.

[0124] For a given time t, in order to approximate the characteristic function At the calculation point x d The function value and derivative value at , here consider the following m-order Taylor polynomial approximation:

[0125]

[0126] To determine the function value in the Taylor polynomial approximation and the derivative value Select the point containing x d The appropriate subdomain of This subdomain is often referred to as x d Suppose there are n in the support domain. x randomly distributed "gridless" nodes {x i}, the characteristic function values ​​at these nodes are

[0127] The gridless function interpolation in S2 includes the following steps:

[0128] B1: Assume that at x d There are n in the support domain of x randomly distributed gridless nodes {x i}, the characteristic function value at this node is Get Taylor polynomial approximation;

[0129] B2: Use the weighted least squares method to solve the Taylor polynomial approximation and obtain the smooth characteristic function at the starting point x of the fluid micro-group. d The value of χ(x d ,t), and through χ(x d ,t)Update the smooth characteristic function.

[0130] The Taylor polynomial approximation in B1 is:

[0131]

[0132] α=(α1,…,α n )

[0133] in, is the Taylor polynomial approximation, α is the multi-index vector in the n-dimensional case, is the function value in the Taylor polynomial approximation, are the derivative values ​​in the Taylor polynomial approximation, α1, α2, …, α n is the multi-index component in the n-dimensional case, i represents the i-th gridless node; here n x Not less than To ensure that the problem is in R n The well-posedness of .

[0134] At this time, the function value and the derivative value It can be obtained by minimizing the functional in the sense of weighted least squares:

[0135]

[0136] Among them, {w i (x d)} is the weight function associated with the node. There are many ways to take the weight function, for example, it can be taken as the following cubic spline weight function:

[0137]

[0138] Where r = |x i -x d | / r i , r i For node x i The corresponding influence radius.

[0139] For convenience, the multiple indices are arranged as {α1,α2,…,α M}, where all neighbors α k and α k+1 (k=1,…,M-1) satisfies one of the following two conditions:

[0140] i)|α k |<|α k+1 |;

[0141] ii)|α k |=|α k+1 |, and there exists i such that and in Represents |α k The i-th component of .

[0142] definition:

[0143]

[0144] Then the Taylor polynomial approximation is:

[0145]

[0146] Where c(·) and p(·) are set vectors, and the superscript T represents the transpose of the matrix.

[0147] At this point, the unknown function value and its derivative value Given by:

[0148]

[0149] where c k is the vector c(x d ). Using the above notation, the weighted least squares problem becomes:

[0150]

[0151] The Taylor polynomial approximation solved by weighted least squares method in B2 is:

[0152]

[0153] Among them, J(·) is the objective function of the least squares problem, min represents the minimum value, and w i (·) is the weight function associated with the node, P is the basis function matrix, W(·) is the weight function matrix, is the characteristic function node value vector, Φ is the shape function matrix, p(·) is the basis function vector, i=1,…,n x ;

[0154] set up is a unit coordinate vector, whose kth component is equal to 1, using e k The vector c(x d ) is expressed analytically as Thus:

[0155]

[0156] Φ=(P T W(x d )P) -1 P T W(x d )

[0157] The first row of Φ corresponds to the function value The kth row corresponds to the derivative value

[0158] Among them, α k is the kth multi-index vector in ascending order, M is the total number of components, and m is the order of the Taylor polynomial;

[0159] Then the smooth characteristic function is at the starting point x of the fluid micro-group d The value of χ(x d ,t) is:

[0160]

[0161] Among them, e1 is the first component of the unit coordinate vector.

[0162] Numerical errors caused by back-tracing along the characteristic lines and mesh-free function interpolation may lead to numerical dissipation of the characteristic function. Therefore, during the motion and deformation process, the characteristic function will gradually lose its characteristic properties and easily identify the wrong interface. In order to maintain long-term accuracy, a characteristic function reinitialization technique is required to restore the dissipated characteristic function to its original form without changing the current interface shape and position.

[0163] Since the smooth characteristic function is defined by using the signed distance function as an auxiliary variable, the reinitialization idea of ​​the present invention is: first reconstruct the signed distance function near the interface, and then redefine the characteristic function.

[0164] In S3, the updated smooth characteristic function is reinitialized, and the formula is:

[0165]

[0166] The signed distance function is an auxiliary variable, and only its value near the interface is necessary when smoothing the characteristic function. Therefore, the signed distance function can be reconstructed as:

[0167]

[0168] Consider a dissipative characteristic function Assuming the interface is still accessible via To identify, in the interface area Using the definition of the first-order directional derivative along the gradient direction, the distance between any point x and the interface can be calculated by l;

[0169]

[0170] The gradient term It can be directly calculated by the above meshless function interpolation. Note that this formula is only valid in the interface region because in the fluid region or At this time, the denominator

[0171] in, is the reconstructed signed distance function, η is the fluid smoothness intensity, is the characteristic function of dissipation, l is the distance between any point x and the interface, is the gradient.

[0172] It can be proved that the characteristic function reinitialization method proposed in the present invention completely preserves the position of the interface and the conservation of the fluid. It is easy to obtain the reconstructed signed distance function The reinitialized characteristic function χ(x, t) = 0.5 is obtained, which means that the interface position remains unchanged. have Therefore, χ(x,t)>0.5, which means that the type of fluid has not changed. For the other fluid, the situation is similar.

[0173] In one embodiment of the present invention, in order to verify the effect of the reinitialization method, the following is considered: Figure 3The model problem shown. Assume that at the initial moment, there is a circular droplet, which is wrapped by another fluid. Under a given periodic flow field, the two fluids begin to flow, and after one cycle, they return to their respective initial positions. Figure 4-6 The simulation results before and after reinitialization are given. It can be seen that when reinitialization is not used, the simulation results are seriously distorted; after reinitialization, the simulation results are almost completely consistent with the theoretical results (initial state).

[0174] In this embodiment, in order to verify the effectiveness of this solution, the following four test examples are performed:

[0175] Test case 1.

[0176] The problem is Figure 7 As shown. At the initial moment, a notched disc-shaped droplet is surrounded by another fluid.

[0177] Two fluids in a given velocity field

[0178] v(x,y,t)=π((0.5-y),(x-0.5))

[0179] The surface rotates counterclockwise around the center of the area and returns to the initial position after T=2 seconds. Figure 8 The simulation results of rotating one cycle under different node distributions are shown in Figure 1. From left to right, they are: the simulation results under uniform 65×65 points, the simulation results under uniform 129×129 points, and the simulation results under uniform 257×257 points. Fig. 9 The simulation results are given when the rotation is 2, 3, and 4 cycles. For this example, many methods have difficulty in giving good results when the rotation is 1 cycle. However, the present invention can still give very good results when the rotation is 4 cycles.

[0180] Test case 2.

[0181] The problem is Fig.10 As shown in Figure 1. At the initial moment, a circular droplet is surrounded by another fluid. The two fluids have a given velocity field.

[0182] v(x,y,t)=πcos(πt / T)(-sin 2 (πx)sin(2πy),sin 2 (πy)sin(2πx))

[0183] Shear stretching is performed, the first T / 2 time is clockwise rotation, the second T / 2 time is counterclockwise rotation, and after T time, it returns to the initial position. Fig.11The simulation results are given when T=2, 3, and 4 (rotation periods 2π, 3π, and 4π). For this example, most methods have difficulty in giving good results when T=2. However, the present invention can still give very good results when T=4.

[0184] Test case 3.

[0185] This problem can be considered as Example 1 in three dimensions. At the initial moment, a spherical droplet with a notch is surrounded by another fluid. The two fluids move in a given velocity field.

[0186] v(x,y,z,t)=π((0.5-y),(x-0.5),0)

[0187] It rotates around the z-axis and returns to its initial position after T=2 seconds. Fig.12 The simulation results of one rotation cycle are given. From left to right, they are: the simulation results under uniform 65×65×65 points, the simulation results under uniform 129×129×129 points, and the simulation results under uniform 257×257×257 points. For this example, most methods are difficult to give good results. However, the present invention can not only give very good results, but also is very efficient.

[0188] Test case 4.

[0189] This problem can be considered as Example 2 in three dimensions. At the initial moment, a spherical droplet is surrounded by another fluid. The two fluids have a given velocity field.

[0190]

[0191] Shear stretching is performed under the condition of T time, and then it returns to the initial position. Fig.13 The simulation results when T=3 are given. From top to bottom, they are: the simulation results under uniform 129×129×129 points, the simulation results under uniform 257×257×257 points, and the simulation results under uniform 513×513×513 points. For this example, most methods are difficult to give good results. However, the present invention can not only give very good results, but also is very efficient.

[0192] The present invention aims at the two-phase flow problems that are widely present in the fields of chemical production, petroleum industry, metallurgy and mining, thermal energy and energy engineering, and proposes a new two-phase flow interface capture method. The method does not require explicit solution of partial differential equations, and only involves some algebra and interpolation calculations, which is simpler and more practical than existing methods. The new method includes four core steps: initialization of characteristic functions, reverse tracing along characteristic lines, gridless-based function interpolation, and reinitialization of characteristic functions, and has the advantages of both characteristic line method and gridless method. The reverse tracing along characteristic lines makes the new method naturally conservative, and the gridless-based function interpolation makes it inherit the advantages of the gridless method, such as high precision, good flexibility, easy adaptability, high parallel efficiency, and applicability to complex areas. These advantages make the overall performance of the new method better than existing methods in solving flow interface problems, thereby providing a simple and practical new method for solving various flow interface problems in scientific research and engineering applications.

[0193] Those skilled in the art will appreciate that the embodiments described herein are intended to help readers understand the principles of the present invention, and should be understood that the protection scope of the present invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific variations and combinations that do not deviate from the essence of the present invention based on the technical revelations disclosed by the present invention, and these variations and combinations are still within the protection scope of the invention.

Claims

1. A novel two-phase flow interface capture method, characterized in that: The following steps are involved: S1: Select the smooth characteristic function of the area occupied by one fluid as the indicator function of the two fluids and initialize it; S2: Based on the initialization results, the smooth characteristic function is updated after the interface position changes by using reverse tracing along the characteristic line and meshless function interpolation; S3: Reinitialize the updated smooth characteristic function to complete the capture of the two-phase flow interface.

2. A novel two-phase flow interface capture method according to claim 1, characterized in that: The initialization formula in S1 is: Among them, χ(·) is the smooth characteristic function, x is the position of the fluid micro-group, φ(·) is the signed distance function, and η is the fluid smoothness intensity.

3. A novel two-phase flow interface capture method according to claim 1, characterized in that: The reverse tracking along the characteristic line in S2 includes the following sub-steps: A1: Calculate the new position of the fluid micro-group along the characteristic line of the flow field; A2: Based on the new position reached by the fluid micro-group and the displacement along the characteristic path, the starting point position of the fluid micro-group is obtained; A3: Update the smooth characteristic function based on the relationship between the starting position of the fluid cluster and the new position reached along the characteristic line of the flow field.

4. A novel two-phase flow interface capture method according to claim 3, characterized in that: The new position reached by the fluid micro-group in A1 along the characteristic line of the flow field is: x a =x d +v(x m ,t+Δt / 2)Δt Among them, x a is the new position reached by the fluid cluster along the characteristic line of the flow field, x d is the starting point of the fluid microgroup, v(·) is the flow velocity, x m is the midpoint of the characteristic path, t is the time, and Δt is the time interval; Under the assumption of uniform motion, the midpoint x of the characteristic path m for: Under the assumption of uniform accelerated motion, the midpoint x of the characteristic path m for: a=(v(x a ,t+Δt)-v(x d ,t)) / Δt Where a is the estimated value of acceleration.

5. A novel two-phase flow interface capture method according to claim 4, characterized in that: The starting point position of the fluid microgroup in A2 is: The displacement d of the fluid particle along the characteristic path is: d=v(x a -d / 2,t+Δt / 2)Δt The fixed point iteration is: d [k+1] =v(x a -d [k] / 2,t+Δt / 2)Δt Among them, d [k+1] is the displacement after the k+1th iteration, d [k] is the displacement after the kth iteration; The starting point of the fluid cluster is: x d =x a -d。 6. A novel two-phase flow interface capture method according to claim 5, characterized in that: The formula for updating the smooth characteristic function in A3 is: x(x a ,t+Δt)=χ(x d ,t) Among them, χ(·) is the smooth characteristic function.

7. A novel two-phase flow interface capture method according to claim 6, characterized in that: The gridless function interpolation in S2 includes the following steps: B1: Assume that at x d There are n in the support domain of x randomly distributed gridless nodes {x i }, the characteristic function value at this node is Get Taylor polynomial approximation; B2: Use the weighted least squares method to solve the Taylor polynomial approximation and obtain the smooth characteristic function at the starting point x of the fluid micro-group. d The value of χ(x d ,t), and through χ(x d ,t)Update the smooth characteristic function.

8. A novel two-phase flow interface capture method according to claim 7, characterized in that: The Taylor polynomial approximation in B1 is: α＝(α1,…,α n ) in, is the Taylor polynomial approximation, α is the multi-index vector in the n-dimensional case, is the function value in the Taylor polynomial approximation, are the derivative values ​​in the Taylor polynomial approximation, α1, α2, …, α n is the multi-index component in the n-dimensional case, i represents the i-th gridless node; definition: Then the Taylor polynomial approximation is: Where c(·) and p(·) are set vectors, and the superscript T represents the transpose of the matrix.

9. A novel two-phase flow interface capture method according to claim 8, characterized in that: The Taylor polynomial approximation solved by weighted least squares method in B2 is: Among them, J(·) is the objective function of the least squares problem, min represents the minimum value, and w i (·) is the weight function associated with the node, P is the basis function matrix, W(·) is the weight function matrix, is the characteristic function node value vector, Φ is the shape function matrix, p(·) is the basis function vector, i=1,…,n x ; set up is a unit coordinate vector, whose kth component is equal to 1, using e k The vector c(x d ) is expressed analytically as Thus: Φ=(P T W(x d )P) -1 P T W(x d ) The first row of Φ corresponds to the function value The kth row corresponds to the derivative value Among them, α k is the kth multi-index vector in ascending order, M is the total number of components, and m is the order of the Taylor polynomial; Then the smooth characteristic function is at the starting point x of the fluid micro-group d The value of χ(x d ,t) is: Among them, e1 is the first component of the unit coordinate vector.

10. A novel two-phase flow interface capture method according to claim 9, characterized in that: In S3, the updated smooth characteristic function is reinitialized, and the formula is: in, is the reconstructed signed distance function, η is the fluid smoothness intensity, is the characteristic function of dissipation, l is the distance between any point x and the interface, is the gradient.

Images