A novel two-phase flow interface capturing method

Abstract

The application discloses a novel two-phase flow interface capturing method, relates to the field of two-phase flow interface, and comprises the following steps: selecting a smooth characteristic function of a fluid region as an indicator function of two fluids and initializing; based on the initialization result, updating the smooth characteristic function after the interface position changes by using reverse tracking along a characteristic line and meshless function interpolation; and re-initializing the updated smooth characteristic function to complete the capturing of the two-phase flow interface. The application solves the numerical difficulties of the existing method when solving the irregular calculation region and the problem of interface capturing precision deterioration in long-time simulation.

Description

Technical Field

[0001] This invention relates to the field of two-phase flow interfaces, and in particular to a novel method for capturing two-phase flow interfaces. Background Technology

[0002] Two-phase flow interface problems are widely found in chemical production, petroleum industry, metallurgy and mining, thermal energy, and energy engineering. In chemical production, studying two-phase flow interface problems can lead to a better understanding of flow phenomena within reactors, helping to optimize operating conditions and improve production efficiency and product quality. In the petroleum industry, studying 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 technologies.

[0003] When two physically different and incompatible fluids come into contact, an interface forms between them. Because fluids cannot maintain their shape under any or a given shear force, the interface between them will move and deform during flow. The dynamic behavior of the interface is usually very complex, depending not only on the properties of the fluids but also on the flow conditions.

[0004] To successfully simulate flow interface problems, the primary task is to accurately track or capture the motion and deformation of the interface during fluid flow. Over the past few decades, dozens of numerical techniques have been proposed to track / capture interfaces. From a coordinate system perspective, these methods can be broadly categorized into two groups: Lagrangian methods (variable computational domain) and Eulerian methods (fixed computational domain).

[0005] However, existing methods often face numerical difficulties when solving for irregular computational domains, and the accuracy of interface capture deteriorates during long-term simulations. Summary of the Invention

[0006] To address the aforementioned shortcomings in existing technologies, this invention provides a novel two-phase flow interface capture method that solves the problems of numerical difficulties in solving irregular computational domains and the deterioration of interface capture accuracy during long-term simulations.

[0007] To achieve the above-mentioned objectives, the technical solution adopted by this invention is: a novel two-phase flow interface capture method, comprising the following steps:

[0008] S1: Select a smooth characteristic function of the region occupied by one fluid as the indicator function for both fluids, and initialize it;

[0009] S2: Based on the initialization results, the smooth feature function after the interface position changes is updated by using reverse tracing along the feature 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] Where χ(·) is the smoothness characteristic function, x is the position of the fluid element, φ(·) is the sign distance function, and η is the fluid smoothness intensity.

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

[0015] A1: The new position reached by the computational fluid element along the characteristic line of the flow field;

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

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

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

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

[0020] Where, x a x represents the new position reached by the fluid element along the characteristic line of the flow field. d Let v(·) be the starting point of the fluid element, v(·) be the flow velocity, and x be the starting point of the fluid element. m The midpoint of the feature path is t, where t is 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 uniformly 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 an estimated value of acceleration.

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

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

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

[0030] The fixed point iteration is obtained as follows:

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

[0032] Where, d [k+1] Let d be the displacement after the (k+1)th iteration. [k] This represents the displacement after the k-th iteration;

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

[0034] x d =x a -d.

[0035] Furthermore, the formula for updating the smooth feature function in A3 is as follows:

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

[0037] Where χ(·) is the smooth characteristic function.

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

[0039] B1: Assume that in x d The support domain has n x A randomly distributed gridless node {x i The characteristic function value at this node is} To obtain an approximate value for the Taylor polynomial;

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

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

[0042]

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

[0044] in, The value is a Taylor polynomial approximation, and α is the multiindex vector in the n-dimensional case. These are the function values ​​in the Taylor polynomial approximation. Let α1, α2, ..., α be the derivative values ​​in the Taylor polynomial approximation. n For the multiple index components 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 defined vectors, and the superscript T indicates the transpose of the matrix.

[0050] Furthermore, the approximate value of the Taylor polynomial obtained by using the weighted least squares method in B2 is:

[0051]

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

[0053] set up It is a unit coordinate vector, and its k-th component is equal to 1. Using e k The vector c(x) d The analytical representation of the k-th component of ) is: Therefore:

[0054]

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

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

[0057] Where, α k Let M be the k-th multiindex vector arranged in ascending order, where M is the total number of components and m is the order of the Taylor polynomial.

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

[0059]

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

[0061] Furthermore, in step S3, the updated smooth feature function is reinitialized using the following formula:

[0062]

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

[0064] The beneficial effects of this invention are as follows: 1) The new method does not require explicit solution of partial differential equations, making it simpler and more practical. 2) The new method is a completely meshless method, inheriting the advantage of meshless methods that have no mesh limitations. 3) The new method has inherent conservation properties and features low dissipation and high resolution. 4) The new method can be fully parallelized, exhibiting high computational efficiency and is suitable for large-scale simulations. 5) The new method is applicable to complex geometric regions and irregular node distributions, demonstrating excellent performance in solving complex interface problems. 6) For some highly challenging publicly available test examples, the new method can provide better results. 7) The overall performance of the new method is superior to most existing methods. Attached Figure Description

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

[0066] Figure 2 This is a schematic diagram of a two-phase flow problem.

[0067] Figure 3 Test image for reinitialization effect.

[0068] Figure 4 The image shows the simulation results before reinitialization.

[0069] Figure 5 The image shows the simulation results after re-initialization.

[0070] Figure 6 This is a comparison chart of simulations before and after reinitialization.

[0071] Figure 7 This is a schematic diagram of the notched disk problem in Example 1.

[0072] Figure 8 This is a comparison chart of simulation results for the notched disk problem during one rotation cycle.

[0073] Figure 9 This is a comparison of simulation results for the notched disk problem under a uniform 257×257 point configuration during 2, 3, and 4 rotation cycles.

[0074] Figure 10 This is a schematic diagram of the single-vortex shear problem in Example 2.

[0075] Figure 11 The figure shows a comparison of simulation results for the single-vortex shear problem with periods of 2π, 3π, and 4π.

[0076] Figure 12 This is a comparison of simulation results for a three-dimensional notched spherical droplet problem during one rotation cycle.

[0077] Figure 13 This is a comparison of simulation results for a three-dimensional single-vortex shear problem with a period of 3. Detailed Implementation

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

[0079] like Figure 1 As shown, a novel method for capturing two-phase flow interfaces includes the following steps:

[0080] S1: Select a smooth characteristic function of the region occupied by one fluid as the indicator function for both fluids, and initialize it;

[0081] S2: Based on the initialization results, the smooth feature function after the interface position changes is updated by using reverse tracing along the feature 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 ​​this invention is to select a smooth characteristic function of the region occupied by one fluid as an indicator function for two fluids. By tracing back along the characteristic line and using meshless function interpolation, the characteristic function is updated after the interface position changes. A re-initialization technique is used to maintain the properties of the characteristic function during the flow process. Since the fluid type remains unchanged for a given fluid mass during the flow process, the value of the characteristic function is constant along the characteristic path / line. This is the core foundation of the method, giving it a clear physical meaning or interpretation. For a new time step, the starting point is found by tracing back along the corresponding characteristic path, and the function value is interpolated to the starting point to obtain the characteristic function value at each calculation point. Furthermore, to reduce the interpolation difficulties encountered by mesh-based methods, we use meshless function interpolation to obtain the characteristic function value at the starting point, thus making this method a completely meshless method.

[0084] like Figure 2 As shown, consider a computational domain Ω filled with fluid 1 (x∈Ω1) and fluid 2 (x∈Ω2). Assume the two fluids are physically immiscible, with a clear interface Γ between them. Under a given flow field, the two fluids will undergo motion and deformation. To track the motion and deformation of the two fluids, the standard characteristic function is taken as:

[0085]

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

[0087] When directly applying equations to identify two-phase interfaces on a Cartesian grid, the results are often very coarse, exhibiting a jagged shape. To overcome this difficulty, this invention proposes a smoothing feature function to reduce the dependence of interface identification results on the grid. During the initialization of the feature function, the signed distance function φ(x), widely used in level set methods, is introduced as an auxiliary variable. The smoothing intensity η (interface thickness) is set, and the initial value of the smoothing feature function is defined.

[0088] The initialization formula in S1 is:

[0089]

[0090] Where χ(·) is the smoothness characteristic function, x is the position of the fluid element, φ(·) is the sign distance function, and η is the fluid smoothness intensity.

[0091] Typically, η is a small positive constant that depends on the spatial step size h of the grid, and can be taken as η = h or η = 2h. Compared to the standard eigenfunction, the smooth eigenfunction only modifies the function values ​​within a narrow band of width 2η near the interface. As η approaches zero, the smooth eigenfunction converges to the standard eigenfunction.

[0092] The reverse tracing along the feature line in S2 includes the following steps:

[0093] A1: The new position reached by the computational fluid element along the characteristic line of the flow field;

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

[0095] A3: Update the smooth characteristic function based on the relationship between the starting position of the fluid micro-particle 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 phase located at x... d The fluid particle at point ∈Ω(t) will leave its position and move along the characteristic lines of the flow field to a new position x. 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 based on the above formula, x must be derived. d and x a The relationship between them. Because these two points lie on the same characteristic line, therefore:

[0099]

[0100] Generally, the flow velocity v(x,t) along the characteristic line is time-dependent, making the integral on the right-hand side of the equation difficult to calculate analytically. To achieve second-order accuracy, this invention employs a midpoint quadrature formula when calculating the integral on the right-hand side of the equation. This quadrature formula assumes that the velocity is constant throughout the entire time interval [t, t+Δt], and takes the velocity at the midpoint of the time interval.

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

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

[0103] Where, x ax represents the new position reached by the fluid element along the characteristic line of the flow field. d Let v(·) be the starting point of the fluid element, v(·) be the flow velocity, and x be the starting point of the fluid element. m The midpoint of the feature path is t, where t is 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 uniformly 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 an estimated value of acceleration.

[0110] Assume the position x of the fluid element at time t+Δt. a With grid point x i (Arrival Point) Consistent. At time t, one time step prior, each fluid particle is located at its respective origin point x. d In other words, the fluid element at time t moves from x... d Starting from the point, after one time step Δt, it reaches 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 of the fluid element in A2 is:

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

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

[0114] The fixed point iteration is obtained as follows:

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

[0116] Where, d [k+1] Let d be the displacement after the (k+1)th iteration. [k] This represents the displacement after the k-th iteration;

[0117] When setting d=0 as the initial condition, usually only 2 to 3 iterations are sufficient to obtain results with adequate accuracy. Once d is obtained, it can be obtained through x. d =x a -d gives the starting point's location, thus allowing us to use the equation χ(x) a ,t+Δt)=χ(x d The characteristic function is updated by (t). The calculation is very similar for the case of uniformly accelerated motion, and will not be repeated here.

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

[0119] x d =x a -d.

[0120] The formula for updating the smooth feature function in A3 is as follows:

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

[0122] Where χ(·) is the smooth characteristic function.

[0123] Typically starting point x d With grid point x i They do not overlap. At time t, the characteristic function only exists at grid point x. i There is a value at this location, therefore, when using χ(x) a ,t+Δt)=χ(x d When updating the characteristic function, interpolation must first be used to obtain the starting point x. d The function value at that point. There are no special restrictions on the choice of interpolation method here. However, for grid-based methods, for each starting point x... d Previously, it was necessary to traverse all grid cells to determine which grid cell contained the point before interpolation, a process that was very difficult to implement. To reduce the interpolation difficulty caused by the grid, this 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 calculation point x d The function value and derivative value at point m are approximated by the following m-th order Taylor polynomial:

[0125]

[0126] To determine the function values ​​in the Taylor polynomial approximation and derivative value Select the point x d appropriate subdomain This subdomain is usually called x d The support domain. Assume there are n in the support domain. x A randomly distributed "meshless" node {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 in x d The support domain has n x A randomly distributed gridless node {x i The characteristic function value at this node is} To obtain an approximate value for the Taylor polynomial;

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

[0130] The Taylor polynomial approximation value in B1 is:

[0131]

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

[0133] in, The value is a Taylor polynomial approximation, and α is the multiindex vector in the n-dimensional case. These are the function values ​​in the Taylor polynomial approximation. Let α1, α2, ..., α be the derivative values ​​in the Taylor polynomial approximation. n Let i represent the i-th gridless node in the n-dimensional case, where n is the multi-index component. x Not less than To ensure the problem is in R n The well-stability of the medium.

[0134] At this point, the function value and 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 choose this weight function; for example, it can be a cubic spline weight function as follows:

[0137]

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

[0139] For convenience, the multiple indicators are arranged as {α1, α2, ..., α...} M The order of}, 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.

[0142] definition:

[0143]

[0144] Then the Taylor polynomial approximation is:

[0145]

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

[0147] At this point, the unknown function value and its derivative value are... It is given by the following formula:

[0148]

[0149] Where c k It is a vector c(x) d The k-th component of ). Using the above notation, the weighted least squares problem becomes:

[0150]

[0151] The approximate value of the Taylor polynomial obtained by using the weighted least squares method in B2 is:

[0152]

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

[0154] set up It is a unit coordinate vector, and its k-th component is equal to 1. Using e k The vector c(x) d The analytical representation of the k-th component of ) is: Therefore:

[0155]

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

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

[0158] Where, α k Let M be the k-th multiindex vector arranged in ascending order, where M is the total number of components and m is the order of the Taylor polynomial.

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

[0160]

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

[0162] Numerical errors caused by backtracking along feature lines and meshless function interpolation can lead to numerical dissipation of feature functions. Therefore, during motion and deformation, feature functions gradually lose their characteristic properties, making it easy to identify erroneous interfaces. To maintain long-term accuracy, a feature function reinitialization technique is needed to restore the dissipated feature functions to their original form without altering the current interface shape and position.

[0163] Since the smooth feature function is defined using the signed distance function as an auxiliary variable, the re-initialization approach of this invention is: first, reconstruct the signed distance function near the interface, and then redefine the feature function.

[0164] In step S3, the updated smooth feature function is reinitialized using the following formula:

[0165]

[0166] The signed distance function is an auxiliary variable; in smooth eigenfunctions, only its value near the interface is truly essential. Therefore, the signed distance function can be reconstructed as:

[0167]

[0168] Consider a dissipative characteristic function Assuming the interface can still be accessed 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 This can be calculated directly using the meshless function interpolation above. Note that this formula only holds true in the interface region, because in the fluid region... or At this point, the denominator

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

[0172] It can be proven that the characteristic function reinitialization method proposed in this invention completely preserves the interface position and fluid conservation. For the interface... The reconstructed symbolic distance function is easily obtained. This results in a reinitialized characteristic function χ(x,t) = 0.5, implying that the interface position remains unchanged. For fluids... have Therefore, χ(x,t)>0.5, which means that the type of fluid has not changed. The situation is similar for the other fluid.

[0173] In one embodiment of the present invention, in order to verify the effect of the re-initialization method, the following is considered: Figure 3The model problem is illustrated below. Assume that initially, there is a circular droplet surrounded by another fluid. Under a given periodic flow field, the two fluids begin to flow, and after one period, they return to their initial positions. Figure 4-6 The simulation results before and after reinitialization are presented. It can be seen that without reinitialization, the simulation results are severely distorted; after reinitialization, the simulation results are almost completely consistent with the theoretical results (initial state).

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

[0175] Test case 1.

[0176] The question is as follows Figure 7 As shown. Initially, a notched, disk-shaped droplet is enveloped by another fluid.

[0177] Two fluids in a given velocity field

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

[0179] It rotates counterclockwise around the center of the region and returns to its initial position after a time of T=2. Figure 8 The simulation results are presented for one rotation under different node distributions. From left to right, they are: simulation results for a uniform 65×65 point distribution, simulation results for a uniform 129×129 point distribution, and simulation results for a uniform 257×257 point distribution. Figure 9 The simulation results presented are for 2, 3, and 4 rotations. For this example, many methods struggle to produce good results even after 1 rotation. However, this invention still yields excellent results after 4 rotations.

[0180] Test case 2.

[0181] The question is as follows Figure 10 As shown. Initially, a circular droplet is enveloped by another fluid. The two fluids are in 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] The device is subjected to shearing and stretching. It rotates clockwise for the first T / 2 time period and counterclockwise for the next T / 2 time period. After T time period, it returns to the initial position. Figure 11Simulation results are presented for T=2, 3, and 4 (rotation periods of 2π, 3π, and 4π). For this example, most methods struggle to produce good results at T=2. However, this invention still yields very good results at T=4.

[0184] Test case 3.

[0185] This problem can be considered as Example 1 in a three-dimensional scenario. Initially, a spherical droplet with a notch is enveloped by another fluid. The two fluids are 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 a time interval T = 2. Figure 12 Simulation results for one rotation cycle are presented. From left to right: simulation results for a uniform 65×65×65 point layout, a uniform 129×129×129 point layout, and a uniform 257×257×257 point layout. For this example, most methods struggle to produce good results. This invention, however, not only yields excellent results but is also highly efficient.

[0188] Test case 4.

[0189] This problem can be considered as Example 2 in a three-dimensional scenario. Initially, a spherical droplet is enveloped by another fluid. The two fluids are in a given velocity field...

[0190]

[0191] It undergoes shearing and stretching, and returns to its initial position after time T. Figure 13 The simulation results for T=3 are presented. From top to bottom, they are: simulation results under uniform 129×129×129 points, simulation results under uniform 257×257×257 points, and simulation results under uniform 513×513×513 points. For this example, most methods struggle to produce good results. However, this invention not only yields excellent results but is also highly efficient.

[0192] This invention addresses the widespread two-phase flow problems in chemical production, petroleum industry, metallurgy and mining, thermal energy and energy engineering, proposing a novel method for capturing two-phase flow interfaces. This method does not require explicit solution of partial differential equations, involving only some algebraic and interpolation calculations, making it simpler and more practical than existing methods. The new method comprises four core steps: initialization of characteristic functions, backtracking along characteristic lines, meshless function interpolation, and reinitialization of characteristic functions, combining the advantages of both the characteristic line method and the meshless method. Backtracking along characteristic lines endows the new method with inherent conservation properties, while meshless function interpolation inherits the advantages of high accuracy, flexibility, adaptability, high parallel efficiency, and applicability to complex regions found in meshless methods. These advantages enable the new method to outperform existing methods in solving flow interface problems, thus providing a simple and practical new approach for solving various flow interface problems in scientific research and engineering applications.

[0193] Those skilled in the art will recognize that the embodiments described herein are intended to help the reader understand the principles of the invention, and should be understood that the scope of protection of the invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific modifications and combinations based on the technical teachings disclosed in this invention without departing from the spirit of the invention, and these modifications and combinations are still within the scope of protection of the invention.

Claims

1. A novel method for capturing the interface of a two-phase flow, characterized in that, Includes the following steps: S1: Select a smooth characteristic function of the region occupied by one fluid as the indicator function for both fluids, and initialize it; The initialization formula in S1 is: in, For smooth characteristic functions, The location of the fluid element. The symbolic distance function, For fluid smoothness strength; S2: Based on the initialization results, the smooth feature function after the interface position changes is updated by using reverse tracing along the feature line and meshless function interpolation; The reverse tracing along the feature line in S2 includes the following steps: A1: The new position reached by the computational fluid element along the characteristic line of the flow field; A2: Based on the new position reached by the fluid micro-element and its displacement along the characteristic path, the starting position of the fluid micro-element is obtained; A3: Update the smooth characteristic function based on the relationship between the starting position of the fluid micro-particle and the new position reached along the characteristic lines of the flow field; The gridless function interpolation in S2 includes the following steps: B1: Assuming in Supported domains include Randomly distributed gridless nodes The characteristic function value at this node is This yields an approximate value using the Taylor polynomial. B2: Use the weighted least squares method to solve for the Taylor polynomial approximation and obtain the position of the smooth characteristic function at the starting point of the fluid micro-element. value at and through Update the smooth feature function; S3: Reinitialize the updated smooth characteristic function to capture the two-phase flow interface; In step S3, the updated smooth feature function is reinitialized using the following formula: in, For the reconstructed symbolic distance function, The characteristic function of dissipation, For any point Distance between the interface and the screen For gradient, For time.

2. The novel two-phase flow interface capture method according to claim 1, characterized in that, The new position reached by the fluid element in A1 along the characteristic line of the flow field is: in, The new position reached by the fluid particle along the characteristic line of the flow field. The origin of the fluid element. For flow velocity, The midpoint of the feature path. For time, For time intervals; Under the assumption of uniform motion, the midpoint of the characteristic path for: Under the assumption of uniformly accelerated motion, the midpoint of the characteristic path for: in, This is an estimate of the acceleration.

3. A novel two-phase flow interface capture method according to claim 2, characterized in that, The starting point of the fluid element in A2 is: Displacement of fluid element along characteristic path for: The fixed point iteration is obtained as follows: in, For the first The displacement after the next iteration. For the first The displacement after the next iteration; The starting point of the fluid element is: 。 4. A novel two-phase flow interface capture method according to claim 3, characterized in that, The formula for updating the smooth feature function in A3 is as follows: in, It is a smooth characteristic function.

5. A novel two-phase flow interface capture method according to claim 4, characterized in that, The Taylor polynomial approximation value in B1 is: in, This is an approximation of the Taylor polynomial. for Multiple index vectors in the 3D case, These are the function values ​​in the Taylor polynomial approximation. This is the derivative value in the Taylor polynomial approximation. for Multiple indicator components in the case of dimension Indicates the first One gridless node; definition: Then the Taylor polynomial approximation is: in, and For a given vector, the superscript This represents the transpose of a matrix.

6. A novel two-phase flow interface capture method according to claim 5, characterized in that, The approximate value of the Taylor polynomial obtained by using the weighted least squares method in B2 is: in, Let the objective function be the least squares problem. This represents the minimum value. The weight function associated with the node, For the basis function matrix, The weight function matrix, The feature function node value vector, For the shape function matrix, For the basis function vector, ; set up For unit coordinate vectors, its first... Each component equals 1, using vector The Each component is analytically represented as Therefore: The first row corresponds to the function value , No. The corresponding derivative value of the row ; in, The first in ascending order Multiple indicator vectors, The total number of components, Let be the order of the Taylor polynomial; The smooth characteristic function is located at the starting point of the fluid element. value at for: in, It is the first component of the unit coordinate vector.

Images