Two-phase channel turbulence simulation method

Abstract

The invention relates to a two-phase channel turbulence simulation method, which comprises the following steps that: an initial channel turbulence field is established, and the initial channel turbulence field comprises a first direction initial speed, a second direction initial speed and a third direction initial speed; under the condition that at least one liquid drop is added at a specific position, a liquid drop channel turbulent flow field of the at least one liquid drop is initialized; solving the phase field model by using a discrete unified gas kinematics format (DUGKS), and processing a convection term by using a weighted essential oscillation-free format (WENO) to describe two-phase evolution of the turbulent flow field of the liquid drop channel; the first direction target speed, the second direction target speed and the third direction target speed of any volume unit and the sequence parameters of the two-phase flow are obtained through calculation. According to the method, the phase field model is solved by using the DUGKS method, and the WENO format is combined, so that the interface physical and microscopic physical effects can be better captured, and the accuracy and reliability of two-phase flow numerical simulation are improved.

Description

Technical Field

[0001] The invention relates to the technical field of computer simulation, and in particular to a method for simulating two-phase channel turbulence. Background Art

[0002] At present, most numerical simulations of fluid mechanics are based on the Navier-Stokes (NS) equations, which are based on the continuum assumption and mainly focus on macroscopic fluid behavior, but cannot effectively capture the fluid interface dynamics at the mesoscopic scale. In two-phase flow, physical phenomena such as interfacial tension and interfacial dynamics play an important role at the mesoscopic scale. The traditional NS equations cannot specifically deal with physical phenomena at this scale, so there are difficulties in simulating complex interfaces in two-phase flow. The existing technology lacks effective tools to deal with physical phenomena at this scale, and usually needs to rely on additional interface tracking or reconstruction techniques. This not only increases the complexity of the simulation, but also easily introduces numerical errors, resulting in serious impacts on the accuracy and stability of the simulation results.

[0003] To deal with turbulence problems under high Reynolds number conditions, researchers usually use three main methods: Reynolds average method (RANS), large eddy simulation (LES) and direct numerical simulation (DNS). Each method has different capabilities in resolving turbulence scales. The Reynolds average method simplifies calculations by modeling turbulent stresses, but reduces accuracy; large eddy simulation models small-scale pulsations and can achieve a certain balance between computational effort and accuracy; while direct numerical simulation does not require modeling and can resolve all scales of turbulence, thus providing the most accurate flow field information. However, the computational cost of DNS is extremely high, which limits its scope of application.

[0004] From a microscopic perspective, fluids are composed of a large number of molecules, so the behavior of fluids can be simulated by molecular kinetic theory. However, due to the large number of molecules, direct simulation of molecular motion is almost impossible with current computing power. In order to find a balance between the macroscopic and microscopic, mesoscopic methods (such as the lattice Boltzmann method LBM, the gas kinetic method GKS, and the unified gas kinetic method UGKS) came into being. These methods capture the details of the flow field from a mesoscopic perspective by solving the Boltzmann equation, which not only avoids the modeling of physical quantities, but also can handle complex flow problems such as multiphase flow and multiscale flow.

[0005] Discrete unified gas kinetic method (DUGKS) is a newer mesoscopic method. Since it was proposed in 2013, it has gradually attracted widespread attention because it combines the advantages of LBM and UGKS. DUGKS uses the finite volume method to deal with particle transport and collision problems, overcoming the limitations of LBM on irregular grids and multi-scale problems. For example, DUGKS can capture multi-frequency vortex shedding more accurately than LBM, and shows superior performance in dealing with single-phase flow and turbulence problems. However, despite the excellent performance of DUGKS in single-phase flow and turbulence simulation, its application in droplet channel turbulence is still in the exploratory stage. In the simulation of multiphase flow and long-term evolution, many existing two-phase flow simulation algorithms still have shortcomings in mass conservation and numerical stability, and urgently need to be further optimized to cope with complex droplet dynamic behavior and its interaction with turbulence.

[0006] On the other hand, in the numerical simulation of long-term evolution, how to maintain the mass conservation of droplets or interfaces has always been a major challenge. Although the Cahn-Hilliard (CH) equation can guarantee mass conservation in theory, in actual numerical simulations, factors such as the truncation operation of the order parameter, the dissipative properties of the diffusion interface, and numerical format errors often lead to the phenomenon of mass non-conservation. In order to compensate for this defect, some existing methods have introduced means such as source term adjustment, order parameter profile correction, or improved numerical algorithms. However, these methods usually rely on additional parameter adjustments and will significantly increase the computational overhead. In addition, the mass conservation effect is often highly sensitive to the choice of parameters and cannot be stably applied in various simulation scenarios.

[0007] The numerical simulation of two-phase fluids faces the problem of how to accurately deal with complex fluid interfaces, especially under turbulent conditions, where the deformation, fragmentation and merging of interfaces are particularly complex. Traditional macroscopic methods mainly include the forward tracking (FT) method, the volume-of-fluid (VOF) method and the level set (LS) method. In the FT method, the dynamic behavior of the droplet or bubble interface is described by tracking its position change, but when the interface breaks or merges, manual intervention is often required, which increases the computational complexity and error. The VOF method can better maintain mass conservation when dealing with interface tracking, but its capture accuracy of the interface shape is heavily dependent on the grid resolution, so it is difficult to deal with complex interface dynamics in turbulence, especially when the grid is coarse, the interface distortion problem is significant. The LS method describes the interface by introducing a distance function and can accurately track the evolution of complex interfaces. However, due to the frequent reinitialization process, it is easy to cause mass conservation problems. Although this method has a good performance on the smoothness of the interface, in large-scale turbulent two-phase flow simulations, it will still increase the computational overhead due to the complex interface reconstruction steps.

[0008] In view of the shortcomings of traditional macroscopic methods, mesoscopic methods have gradually been applied to the numerical simulation of two-phase flow in recent years. By describing the evolution of fluid interface from the mesoscopic scale, these methods can better deal with the complex physical phenomena of the two-phase interface. At present, common mesoscopic methods include color gradient model, pseudopotential model and free energy model. These methods deal with the two-phase interface problem from a multi-scale perspective by rationally designing the Boltzmann equation.

[0009] Among them, the free energy-based phase field model shows significant advantages in dealing with interface problems. This model introduces an order parameter to describe the interface between two fluids, and combines the Cahn-Hilliard (CH) equation to characterize the evolution of the order parameter. It has thermodynamic consistency and a robust physical basis, and can accurately describe the evolution of the two-phase flow interface in complex flows. However, even with kinetic methods, maintaining the mass conservation of droplets in long-term simulations is still a major challenge, especially when dealing with the deformation, breakup and merging of droplets. Existing methods often have the problem of numerical error accumulation. Summary of the invention

[0010] In order to overcome the shortcomings of the prior art, according to the present invention, a method for simulating two-phase channel turbulence is provided, characterized in that it includes the following steps:

[0011] Establishing an initial channel turbulence field, wherein the initial channel turbulence field includes an initial velocity u0 in a first direction, an initial velocity v0 in a second direction, and an initial velocity w0 in a third direction;

[0012] In case of adding at least one droplet at a specific position, initializing a droplet channel turbulent flow field of the at least one droplet;

[0013] The phase field model is solved by using a discrete unified gas kinetic scheme, and the convection term is processed by using a weighted essentially non-oscillatory scheme to describe the evolution of the turbulent flow field in the droplet channel;

[0014] The first direction target velocity u, the second direction target velocity v, the third direction target velocity w and the sequence parameter φ of the two-phase flow of any volume unit are calculated.

[0015] According to one embodiment, the method further comprises:

[0016] The first order parameter of the initial channel turbulence field is φ A The second order parameter of the droplet channel turbulence field is φ B In the case of two fluids, the two-phase interface is defined, where the two-phase interface is, and its third-order parameter is equivalent to (φ A +φ B ) / 2; and

[0017] Using the first-order parameter φ A and the second-order parameter φ B , the density ρ and viscosity μ of the target volume unit are calculated in the phase field model:

[0018]

[0019] According to one embodiment, the method further comprises:

[0020] Using the first-order parameter φ A and the second-order parameter φ B , the volume free energy density Γ(φ) and the total free energy F are calculated in the phase field model:

[0021] Γ(φ)=β(φ-φ A ) 2 (φ-φ B ) 2

[0022]

[0023] Here, κ is a constant representing the relative size of bulk free energy, and β is a constant representing the relative size of interfacial free energy.

[0024] According to one embodiment, the method further comprises:

[0025] The calculated total free energy F is subjected to first-order variation processing:

[0026]

[0027] The total chemical potential μ is calculated φ .

[0028] According to one embodiment, the method further comprises:

[0029] Using the total free energy F and the calculated total chemical potential μ φ , and process it as follows:

[0030]

[0031] Where W represents the thickness of the interface between the two phases, σ represents the surface tension between the two phases, and φ eq represents the distribution of the equilibrium order parameter φ.

[0032] According to one embodiment, the step of using a discrete unified gas kinetic format and a coupled phase field model to describe the two-phase evolution of the initial channel turbulence field and the droplet channel turbulence field comprises:

[0033] The two-phase evolution is driven by the chemical potential gradient:

[0034]

[0035] Among them, M CH is the migration rate.

[0036] According to one embodiment, the method further comprises:

[0037] The two-phase evolution driven by the chemical potential gradient is coupled with the hydrodynamic equations:

[0038]

[0039] To describe the evolution of the hydrodynamic force and the order parameter φ, where p is the pressure, F is the volume force, and γ = (ρ A -ρ B ) / (φ A ρ B -φ B ρ A ) is a constant associated with the order parameter and density of the two phases.

[0040] According to one embodiment, the distribution function f of the hydrodynamic force and order parameter is α , g α Described as:

[0041]

[0042] in, Indicates equilibrium state, represents the source terms, which are obtained by inverse design.

[0043] According to one embodiment, the method further comprises:

[0044] The optimization algorithm is used to improve the mass conservation and the capture accuracy of the two-phase interface in the two-phase channel turbulence simulation.

[0045] According to an embodiment, in the two-phase channel turbulent flow, the volume units are divided non-uniformly or evenly.

[0046] According to one embodiment, the convection terms in the discrete unified gas kinetic scheme are processed using a weighted essentially non-oscillatory (WENO) scheme.

[0047] In view of the problem that traditional methods cannot accurately simulate the physics of the two-phase interface, the present invention introduces the discrete unified gas kinetics (DUGKS) method to deal with the physical phenomena of the two-phase interface at the mesoscopic scale, thus overcoming the limitations of traditional methods in simulating complex interfaces. The DUGKS method combines the kinetic properties of the Boltzmann equation with the advantages of the finite volume method, and can achieve high-precision interface tracking in irregular grids and multi-scale flows. Compared with traditional macro-scale-based simulation methods, the present invention avoids additional tracking or reconstruction operations on the interface, effectively reduces numerical errors, and enhances the accuracy of interface physical phenomena (such as interfacial tension and interface dynamics).

[0048] The present invention uses a mesoscopic method combined with a phase field model to design a two-phase flow mesoscopic model, which can better capture the interface physics and microscopic physical effects from a mesoscopic perspective and improve the accuracy and reliability of the simulation; the present invention also uses the WENO format in the convection term to improve the stability of the algorithm and the conservation of mass in long-term multi-scale simulations; through the present invention, the D3Q15A5 model is used to reduce the amount of calculation, thereby improving the calculation efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] The specific embodiments of the present invention are further described in detail below with reference to the accompanying drawings, wherein:

[0050] Figure 1 A schematic diagram of the process of simulating a two-phase channel turbulence according to the present invention is shown;

[0051] Figure 2 A schematic flow chart of the DUGKS algorithm for two-phase channel simulation based on the phase field model according to the present invention is shown;

[0052] Figure 3 Result diagram showing the morphology and coherent structure of a single droplet in channel turbulence with and without WENO format processing according to the present invention;

[0053] Figure 4 The flow field cross-section diagram of a single droplet at the same time and position in the channel turbulence with and without WENO format processing according to the present invention is shown;

[0054] Figure 5 A graph showing the change in droplet mass over time of a single droplet in channel turbulence using and not using the WENO format according to the present invention;

[0055] Figure 6 The result diagrams showing the initialization morphology, position and background flow velocity cloud diagram of multiple droplets in channel turbulence according to the present invention;

[0056] Figure 7Result diagrams showing the droplet morphology and background flow velocity cloud diagram at different times using different discrete velocity models according to the present invention; and

[0057] Figure 8 The graph shows the change of droplet mass over time using different discrete velocity models according to the present invention. DETAILED DESCRIPTION

[0058] The present invention will be described in detail below in conjunction with the accompanying drawings and specific embodiments. The exemplary embodiments and descriptions of the present invention are used to explain the present invention, but are not intended to limit the present invention.

[0059] The method of the present invention uses the DUGKS method combined with the phase field method for simulating two-phase channel turbulence. The embodiments of the present invention will be further described in detail below in conjunction with the accompanying drawings. Figure 1-2 As shown, they respectively show a flow chart of the two-phase channel turbulence simulation method according to the present invention and a flow chart of the DUGKS algorithm for two-phase flow simulation based on the phase field model.

[0060] The present invention adopts a phase field model to simulate the dynamic behavior of droplets in two-phase channel turbulence. In this model, different order parameters φ are used to mark different fluids, for example, the order parameter φ A and φ B represents two fluids A and B respectively. That is to say, in an ideal case, the internal order parameters of fluids A and B should be φ A and φ B , while the order parameter changes monotonically and smoothly at the interface region between the two phases. Here, the hyperbolic tangent function is used, and the order parameter is equal to (φ A +φ B ) / 2 is taken as the accurate interface between the two phases. In the phase field model, the density ρ and viscosity μ can be expressed as a linear function of the order parameter:

[0061]

[0062] The key to using phase field models in two-phase flow problems is to calculate the order parameter value corresponding to the grid to characterize the interaction between the two phases.

[0063] The phase field model has a good thermodynamic foundation. Here we introduce several thermodynamic concepts related to the phase field model. The free energy F includes two parts: body free energy and interface free energy. The free energy F of the entire two-phase flow system can be expressed as:

[0064]

[0065] The volume free energy density in this patent adopts a double potential well function:

[0066] Γ(φ)=β(φ-φ A ) 2 (φ-φ B ) 2

[0067] Among them, κ and β are both positive constants, reflecting the relative size of the volume free energy and the interface free energy, respectively. Their values ​​depend on the thickness W of the two-phase interface and the size of the surface tension σ in the simulation example. The result of the first-order variation of the free energy is the chemical potential μ of the two-phase flow system. φ , whose expression is:

[0068]

[0069] In a two-phase flow system, the chemical potential distribution of different components (i.e., fluid A and fluid B) is uneven, which will produce a chemical potential gradient. Driven by the chemical potential gradient, the order parameters of the fluids at each location will change, resulting in interaction between the two phases.

[0070] For a flat interface in equilibrium, its chemical potential μ φ is 0, the distribution of the order parameter φ follows the form of the hyperbolic tangent function, and the specific expression is:

[0071]

[0072] in, Indicates the vertical distance of the interface. represents the thickness of the interface between the two phases. In addition, the surface tension between the two phases is

[0073] The present invention uses the above formula as the initialization of the phase field order parameter.

[0074] The Cahn-Hilliard (CH) equation in the phase field model can be understood as the control equation for the change of the order parameter driven by the chemical potential gradient. The present invention couples the CH equation satisfied by the order parameter with the hydrodynamic equation to describe the evolution of the two-phase flow in the following form:

[0075]

[0076] Where γ=(ρ A -ρ B ) / (φ A ρ B -φ B ρ A ) is a constant, which is related to the order parameters and density of the two phases. CH is the mobility in the CH equation.

[0077] The above control equations are solved using DUGKS. This method uses two Boltzmann equations to describe the evolution of hydrodynamic forces and order parameters, and improves mass conservation and interface capture accuracy during the simulation process through an optimization algorithm.

[0078] In this method, the hydrodynamic force is given by the distribution function f α Description, the order parameter is given by the distribution function g α The specific form is as follows:

[0079]

[0080] Among them, the equilibrium state of the distribution function and source term It is necessary to perform inverse design to obtain the equilibrium distribution function and source term used in the inverse design of the present invention in the specific form:

[0081]

[0082] in,

[0083] The relaxation times corresponding to the above two Boltzmann equations are For the above Boltzmann equation, the present invention adopts the DUGKS method to solve it discretely. Here, for the convenience of description, the Boltzmann equation is written in a unified form:

[0084]

[0085] Among them, the collision operator In the actual calculation process, several simplified distribution functions are introduced to facilitate the calculation. and These distribution functions are linear combinations of the original distribution function, the equilibrium distribution function, and the source term:

[0086]

[0087] The DUGKS method has the advantages of good stability and the ability to use non-uniform grids. However, in order to facilitate the description of droplet morphology, a uniform grid is still used. In order to improve the numerical stability and solution accuracy of DUGKS in multiphase flow simulation, the present invention adopts the classic third-order WENO format with t n At time x j+1 / 2 and perpendicular to the x-direction. If the x-component of the discrete particle velocity is j+1 / 2 The position is positive, then the approximate mesh interface results based on different templates are:

[0088]

[0089] Among them, β represents The smoothness, nonlinear weight in ε is a small number used to prevent the denominator from being 0. Here, ε=10 -6 Finally, the grid interface x j+1 / 2 The distribution function at can be reconstructed as the following combination:

[0090]

[0091] At the beginning of the simulation, it is necessary to first import a fully developed channel turbulence field (i.e., fluid A) as the initial flow field. Then, a number of droplets of a specified size (i.e., fluid B) are added at the specified positions in the flow field, which is mainly implemented using the equilibrium state order parameter distribution mentioned above. The DUGKS method is then used to evolve the flow field. The macroscopic physical quantities of interest in the simulation, such as pressure p, velocity u, and order parameter φ, can be calculated by taking the moment of the distribution function:

[0092]

[0093] Here, the third-order WENO is used to process the distribution function on the grid interface, and the code after this processing is called the PF(CH)-DUGKS-WENO code. Here, the difference between the calculation results of the PF(CH)-DUGKS-WENO code and the PF(CH)-DUGKS code is mainly compared.

[0094] Table 1: Parameter settings in the WENO format test of the channel turbulence example of the initial single droplet

[0095]

[0096] As shown in Table 1, the two-phase flow field parameters, N x 、N y 、N z are the number of grids in the x, y, and z directions, Re τ =u τ H / v is the shear Reynolds number, the half channel width H is the selected characteristic length scale, u τ is the friction velocity, v and ρ are the kinematic viscosity coefficient and density of the fluid respectively. is the minimum grid in the wall direction after dimensionless transformation, δt is the simulation time interval, and D0 is the droplet diameter at initialization. H ) and Cahn number (Ch H ) are defined as: As mentioned above, σ is the surface tension and W is the interface thickness.

[0097] First, using single-phase DUGKS, under the same single-phase channel flow physical parameters (i.e., single-phase flow field parameters in Table 1), a steady-state channel turbulence field is obtained; then this velocity field is imported into the PF(CH)-DUGKS code and PF(CH)-DUGKS-WENO code as the initial velocity field of the two-phase channel flow; at the same time, a droplet with a diameter of D0 = 0.8H is placed in the center of the channel. This moment is recorded as the dimensionless time t * =u τ t / H = 0. Coherent structure using velocity gradient tensor The second invariant Q is identified by:

[0098]

[0099] in, They are The antisymmetric part (called the rotation tensor) and the symmetric part (called the strain rate tensor) of is the curl (vorticity) of u. The dimensionless forms of ω and Q are:

[0100]

[0101] Figure 3 The result diagram shows the morphology and coherent structure of a single droplet in channel turbulence with and without WENO format processing according to the present invention, which shows that the present invention has a good performance in t * = 10.8 seconds, that is, the simulation program ran for 6×10 5 The droplet morphology and flow field coherent structure under the step simulation, where the coherent structure is Q * = 0.008 isosurface. Gao Chang et al. found that the coherent structure near the droplet surface is rich when simulating the same example. Figure 3 Rich coherent structures also appear near the middle droplet, which verifies the correctness of the algorithm of the present invention.

[0102] Figure 4 The flow field cross-section diagram of a single droplet at the same time and position in the channel turbulence with and without the WENO format according to the present invention is shown, in particular, also at t * =10.8, which is the section of the droplet at the same time and the same position when using and not using the WENO format. In this embodiment, the Weber number used is We H =0.71, the surface tension is relatively large at this time, and the droplet will not deform too much, so using the WENO format can make the result more reasonable.

[0103] The total mass of the droplet based on the order parameter is defined as:

[0104]

[0105] Figure 5 The graph shows the change of the droplet mass over time in the channel turbulence with and without the WENO format according to the present invention, where the solid line and the dashed line are the evolution of the droplet mass over time with and without the WENO method, respectively. It can be seen that after changing PF(CH)-DUGKS to PF(CH)-DUGKS-WENO, the mass conservation is greatly improved. For example, for the t selected above, * = 10.8, due to the diffusion of the interface, the result of PF(CH)-DUGKS increases by 70.5% compared with the initial value, while the result of PF(CH)-DUGKS-WENO increases by 6.15% compared with the initial value, a difference of more than 10 times. Therefore, it is necessary to use PF(CH)-DUGKS in combination with WENO format to solve two-phase channel turbulence problems.

[0106] Table 2: Parameter settings for the lattice velocity model test of the initial multi-droplet channel turbulence case

[0107]

[0108] Re τ = 100 shows that it is necessary to increase the Reynolds number to obtain a more fully developed turbulence. However, increasing the Reynolds number will significantly increase the amount of calculation. For example, if the same dimensionless time is required, Re τ =150 The calculation amount is about Re τ =5.06 times of 100. At the same time, Re τ =150In the statistical quasi-steady state, the average velocity fluctuation ratio Re τ = 100 is small, the statistic is more stable over time, so the dimensionless time required for statistics may be smaller than Re τ =100 is short. In the mesoscopic method, one of the main reasons for the large amount of calculation is that there are many velocities in the discrete velocity concentration, and the amount of calculation is almost proportional to the number of velocities. In order to speed up the program calculation, it is necessary to reduce the number of grid velocities while keeping the accuracy unchanged. Here, the fifth-order D3Q15A5 model is considered to replace the fifth-order D3Q19A5 model to reduce the number of discrete velocities and thus reduce the amount of calculation.

[0109] The initialization velocity field here adopts the same shear Reynolds number, namely Re τ =150 for the fully developed single-phase channel turbulent flow field. The computational domain, grid and single-phase flow parameters are shown in Table 2.

[0110] Initialize multiple droplets in two-phase channel turbulence, the number is N d= 373 to observe more complex droplet-turbulence interactions. The parameter settings are shown in Table 2 for the two-phase parameters. The two-phase turbulence initialization results are shown in Figure 6 As shown, the initial droplets are evenly distributed in the channel turbulence, which is recorded as t * =0, and Re τ = 150. The PF(CH)-DUGKS-WENO codes of the two discrete velocity sets D3Q19A5 and D3Q15A5 were used for calculations. The difference between these two codes lies only in the discrete particle velocity models.

[0111] The results of the two sets of codes can be found in Figure 7 In the first and second columns, Re τ =150,We H = 1. At t * =2.4,4.8,7.2, the corresponding running steps are 2×10 5 ,4×10 5 ,6×10 5 At the same time, from the results of the total mass of the droplets obtained by using the two speed sets in the present invention (see Figure 8 ), the mass conservation performance of the D3Q15A5 model is slightly better than that of the D3Q19A5. Due to the chaotic characteristics of turbulence and the complex interaction with droplets, the instantaneous results obtained by the two velocity models cannot overlap at these moments. However, from the calculation results obtained using these two different discrete velocity models, namely the order parameter φ and the hydrodynamic parameters (p, u, v, w), the two have obvious similarities: (1) At the same moment, the number and size distribution of droplets are similar; (2) Small droplets are not easily deformed and are close to spherical; large droplets are relatively easier to deform. In general, the D3Q19A5 model and the D3Q15A5 model have no essential effect on the results, and the results are in line with physical laws. The calculations of both are correct, which verifies the correctness of the code.

[0112] The technical solutions provided by the embodiments of the present invention are introduced in detail above. Specific examples are used herein to illustrate the principles and implementation methods of the embodiments of the present invention. The description of the above embodiments is only applicable to help understand the principles of the embodiments of the present invention. At the same time, for those skilled in the art, according to the embodiments of the present invention, there will be changes in the specific implementation methods and application scopes. In summary, the content of this specification should not be understood as limiting the present invention.

Claims

1. A method for simulating two-phase channel turbulence, characterized in that: The following steps are involved: Establishing an initial channel turbulence field, wherein the initial channel turbulence field includes an initial velocity u0 in a first direction, an initial velocity v0 in a second direction, and an initial velocity w0 in a third direction; In case of adding at least one droplet at a specific position, initializing a droplet channel turbulent flow field of the at least one droplet; The phase field model is solved by using a discrete unified gas kinetic scheme, and the convection term is processed by using a weighted essentially non-oscillatory scheme to describe the evolution of the turbulent flow field in the droplet channel; The first direction target velocity u, the second direction target velocity v, the third direction target velocity w and the sequence parameter φ of the two-phase flow of any volume unit are calculated.

2. The method for simulating two-phase channel turbulence according to claim 1, characterized in that: Also includes: The first order parameter of the initial channel turbulence field is φ A The second order parameter of the droplet channel turbulence field is φ B In the case of two fluids, the two-phase interface is defined, where the two-phase interface is, and its third-order parameter is equivalent to (φ A +φ3) / 2; and Using the first-order parameter φ A and the second-order parameter φ B , the density ρ and viscosity μ of the target volume unit are calculated in the phase field model:

3. The method for simulating two-phase channel turbulence according to claim 1, characterized in that: Also includes: Using the first-order parameter φ A and the second-order parameter φ B , the volume free energy density Γ(φ) and the total free energy F are calculated in the phase field model: C(φ)=β(φ-φ A ) 2 (f-f B ) 2 Here, κ is a constant representing the relative size of bulk free energy, and β is a constant representing the relative size of interfacial free energy.

4. The method for simulating two-phase channel turbulence according to claim 3, characterized in that: Also includes: The calculated total free energy F is subjected to first-order variation processing: The total chemical potential μ is calculated φ .

5. The method for simulating two-phase channel turbulence according to claim 4, characterized in that: Also includes: Using the total free energy F and the calculated total chemical potential μ φ , and process it as follows: Where W represents the thickness of the interface between the two phases, σ represents the surface tension between the two phases, and φ eq represents the distribution of the equilibrium order parameter φ.

6. The method for simulating two-phase channel turbulence according to claim 1, characterized in that: The step of using a discrete unified gas kinetic format and a coupled phase field model to describe the two-phase evolution of the initial channel turbulence field and the droplet channel turbulence field comprises: The two-phase evolution is driven by the chemical potential gradient: Among them, M CH is the migration rate.

7. The method for simulating two-phase channel turbulence according to claim 6, characterized in that: Also includes: The two-phase evolution driven by the chemical potential gradient is coupled with the hydrodynamic equations: To describe the evolution of the hydrodynamic force and the order parameter φ, where p is the pressure, F is the volume force, and γ = (ρ A -ρ B ) / (φ A ρ B -φ B ρ A ) is a constant associated with the order parameter and density of the two phases.

8. The method for simulating two-phase channel turbulence according to claim 7, characterized in that: The distribution function of the hydrodynamic force and order parameter is f α , g α Described as: in, Indicates equilibrium state, represents the source terms, which are obtained by inverse design.

9. The method for simulating two-phase channel turbulence according to claim 7, characterized in that: Also includes: The optimization algorithm is used to improve the mass conservation and the capture accuracy of the two-phase interface in the two-phase channel turbulence simulation.

10. The method for simulating two-phase channel turbulence according to claim 1, characterized in that: In the two-phase channel turbulence, the volume unit is divided non-uniformly or uniformly; the convection term in the discrete unified gas kinetic format is processed using the weighted essential non-oscillatory (WENO) format.

Images