Two-phase flow gas kinetic algorithm optimization method, device and equipment and storage medium

Abstract

The invention discloses a two-phase flow gas kinematics format algorithm optimization method, device and equipment and a storage medium. The method comprises the steps of constructing a gas kinematics format algorithm for phase field solution; the phase field interface flux of the phase field gas kinetic format algorithm is constructed on the basis of Chapman-Enskog analysis; constructing a gas kinetic format algorithm for solving the flow field; the flow field interface flux of the flow field gas kinetic format algorithm is constructed on the basis of Chapman-Enskog analysis; and on the basis of the Cahn-Hilliard equation and the Navier-Stokes equation, a two-phase flow numerical algorithm for coupling a phase field and a flow field is realized. According to the method, starting from a Boltzmann-BGK equation, a phase field and flow field coupled gas kinetic format algorithm is constructed under a finite volume framework, and the numerical flux of a macroscopic equation can be accurately recovered through Chapman-Enskog expansion, so that the calculation process is simplified, the calculation efficiency and accuracy are remarkably improved, and the method is suitable for large-scale popularization and application. The method can effectively simulate interface dynamics and flow field evolution in a complex two-phase flow system.

Description

Technical Field

[0001] The present invention belongs to the technical field of fluid mechanics, and specifically relates to a two-phase flow gas kinetic algorithm optimization method, device, equipment and storage medium. Background Art

[0002] Two-phase flow problems are widely present in nature and engineering applications, especially in the fields of energy development and utilization, carbon capture, utilization and storage. Its complex physical phenomena and basic laws have an important impact on engineering design and optimization. At present, common two-phase flow numerical simulation methods are mostly based on the traditional macroscopic Navier-Stokes equations or discrete velocity type mesoscopic methods such as LBM. However, these methods have limitations in dealing with the capture of two-phase interfaces and the simulation of complex flow fields, especially in the accurate description of phase interfaces and the simulation of complex multiphase flows. The gas kinetic scheme (GKS), as a finite volume format numerical method based on the mesoscopic Boltzmann equation, has the advantages of high accuracy, good stability and strong grid flexibility, and has been widely used in the numerical simulation of compressible flows. However, the existing GKS method still has shortcomings in dealing with two-phase flow problems with significant changes in density fields, especially in the application of interface dynamics, and cannot effectively capture and describe complex phase interfaces. This limitation restricts the wide application of the GKS method in the field of multiphase flows. Summary of the invention

[0003] In view of the deficiencies in the prior art, an object of the present invention is to provide a two-phase flow gas kinetic algorithm optimization method, device, equipment and storage medium to solve the problems existing in the prior art.

[0004] According to one aspect of the present application, a two-phase flow gas kinetic algorithm optimization method is disclosed, the method comprising:

[0005] Constructing a gas kinetic format algorithm for phase field solution, wherein the gas kinetic format algorithm for phase field solution includes a continuous Boltzmann-BGK equation with a source term;

[0006] Constructing the phase field interface flux of the gas kinetic format algorithm for the phase field solution based on Chapman-Enskog analysis, wherein the phase field interface flux can be restored to the Cahn-Hilliard phase field equation;

[0007] Constructing a gas kinetic format algorithm for flow field solution, wherein the gas kinetic format algorithm for flow field solution includes a continuous Boltzmann-BGK equation containing a source term;

[0008] Constructing the flow field interface flux of the gas kinetic format algorithm for solving the flow field based on Chapman-Enskog analysis, wherein the flow field interface flux can be restored to the Navier-Stokes equation;

[0009] Based on the Cahn-Hilliard equation and the Navier-Stokes equation, a two-phase flow numerical algorithm is implemented by coupling the phase field and flow field.

[0010] In some embodiments, the gas kinetic format algorithm for phase field solution is constructed, and the gas kinetic format algorithm for phase field solution includes a continuous Boltzmann-BGK equation containing a source term based on the following formula:

[0011]

[0012] Where c is the velocity of the particle distribution function g, τ g is the relaxation time, and G is the phase field source term.

[0013] In some embodiments, the phase field interface flux of the gas kinetic format algorithm for phase field solution is constructed based on Chapman-Enskog analysis, wherein the phase field interface flux is restored to the Cahn-Hilliard phase field equation including:

[0014] Obtaining the phase field equilibrium distribution function and phase field source term of each order moment satisfying the first condition;

[0015] The moments of equilibrium distribution function and source term satisfy the following conditions:

[0016]

[0017] Where η is the relaxation time adjustment coefficient, and the mobility can be expressed as parameter Set to 0, B α and C 0 is the coefficient to be determined. First, use operator, to obtain the expression of the phase field equilibrium distribution function;

[0018] The phase field equilibrium distribution function is expressed based on Chapman-Enskog:

[0019]

[0020] Taking the zero-order moment of the phase field equilibrium distribution function to obtain a zero-order flux equation;

[0021] The zero-order flux equation is:

[0022]

[0023] In the formula, F = <c α ,g>, F is the phase field interface flux of the gas kinetic format algorithm for phase field solution;

[0024] Substituting the phase field equilibrium distribution function and each order moment of the phase field source term into the zero-order flux equation to obtain the first-order phase field flux;

[0025] The first-order phase field flux is:

[0026]

[0027] The first-order phase field flux is subjected to a Chapman-Enskog expansion, and the phase field interface flux is restored to the Cahn-Hilliard phase field equation.

[0028] In some embodiments, the gas kinetic format algorithm for constructing the flow field solution includes a continuous Boltzmann-BGK equation containing a source term based on the following formula:

[0029]

[0030] In the formula, τ f is the relaxation time, S is the flow source term, c is the velocity of the particle distribution function f,

[0031] f eq =ρΓ(u,U) is the equilibrium distribution function;

[0032] ρ is the density of the gas, ρ = <f>;

[0033] The expression of Γ(u,U) is:

[0034]

[0035] In some embodiments, the flow field source term is expressed as:

[0036] The flow field interface flux of the gas kinetic format algorithm for solving the flow field is constructed based on Chapman-Enskog analysis, wherein the flow field interface flux is restored to the Navier-Stokes equations including:

[0037] Sure The moment of the equilibrium distribution function of the flow field and the moment of the flow field source term at :

[0038] Among them, the moment of the flow field equilibrium distribution function is:

[0039]

[0040] The moment of the flow field source term is:

[0041]

[0042] Based on Chapman-Enskog analysis, the zero-order moment and the first-order moment of the gas kinetic format algorithm for solving the flow field are obtained respectively;

[0043] Among them, the zero-order moment of the gas kinetic format algorithm for solving the flow field is:

[0044]

[0045] The first-order moment of the gas kinetic format algorithm for solving the flow field is:

[0046]

[0047] Define a new distribution function:

[0048]

[0049] in,

[0050] Based on the new distribution function, the Boltzmann-BGK expression, the moments of the equilibrium distribution function of the flow field and the moments of the flow field source term are determined.

[0051] The Boltzmann-BGK expression is:

[0052]

[0053] The moments of the equilibrium distribution function of the flow field are:

[0054]

[0055] The moments of the flow field source terms are:

[0056]

[0057] Substituting each order moment of the flow field equilibrium distribution function and the flow field source term into the zero-order moment to obtain the first-order flow field flux;

[0058] Substituting the moments of the flow field equilibrium distribution function and the flow field source term into the first-order moment to obtain the second-order flow field flux;

[0059] The first-order flow field flux and the second-order flow field flux are subjected to Chapman-Enskog expansion, and the flow field interface flux is reduced to the Navier-Stokes equation.

[0060] According to one aspect of the present application, a two-phase flow gas kinetic algorithm optimization device is also disclosed, the device comprising:

[0061] A first building module is used to build a phase-field gas kinetic format algorithm for phase-field solution, wherein the phase-field gas kinetic format algorithm includes a continuous Boltzmann-BGK equation containing a source term;

[0062] A first analysis module is used to construct a phase field interface flux of the gas kinetic format algorithm of the phase field solution based on Chapman-Enskog analysis, wherein the phase field interface flux can be restored to the Cahn-Hilliard equation;

[0063] A second construction module is used to construct a gas kinetic format algorithm for flow field solution, wherein the gas kinetic format algorithm for flow field solution includes a continuous Boltzmann-BGK equation containing a source term;

[0064] A second analysis module is used to construct a flow field interface flux of the gas kinetic format algorithm for solving the flow field based on Chapman-Enskog analysis, wherein the flow field interface flux can be restored to the Navier-Stokes equation;

[0065] The coupling module is used to realize a two-phase flow numerical algorithm by coupling a phase field and a flow field based on the Cahn-Hilliard equation and the Navier-Stokes equation.

[0066] According to another aspect of the present application, an electronic device is also disclosed, which includes a memory and at least one processor, wherein the memory stores instructions; the at least one processor calls the instructions in the memory so that the electronic device executes each step of the two-phase flow gas kinetic algorithm optimization method as described above.

[0067] According to another aspect of the present application, a computer-readable storage medium is also disclosed, on which instructions are stored, characterized in that when the instructions are executed by a processor, each step of the two-phase flow gas kinetic algorithm optimization method as described in any of the above items is implemented.

[0068] The present invention includes but is not limited to the following beneficial effects: Starting from the Boltzmann-BGK equation, the present invention constructs a gas kinetic format algorithm for phase field and flow field coupling under a finite volume framework. Through the Chapman-Enskog expansion, the numerical flux of the macroscopic equation can be accurately restored, the calculation process is simplified, and the calculation efficiency and accuracy are significantly improved. This method can effectively simulate the interface dynamics and flow field evolution in complex two-phase flow systems. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art are briefly introduced below.

[0070] Figure 1 is a flow chart of a two-phase flow gas kinetic algorithm optimization method according to an embodiment of the present application;

[0071] Figure 2 is the density distribution on the centerline of the stationary droplet;

[0072] Figure 3 Verification of Laplace's law for different surface tensions of stationary droplets

[0073] Figure 4 is the interface distribution of three-dimensional Rayleigh-Taylor instability at different times (At=0.5, Re=128);

[0074] Figure 5 The time variation of the bubble front, saddle point and spike position of the three-dimensional Rayleigh-Taylor instability (At=0.5, Re=1024);

[0075] Figure 6 It is a structural block diagram of a two-phase flow gas kinetic algorithm optimization device according to an embodiment of the present application;

[0076] Figure 7 It is a diagram of an electronic device. DETAILED DESCRIPTION

[0077] The terms "first", "second", "third", "fourth", etc. (if any) in the specification and claims of the present invention and the above-mentioned drawings are used to distinguish similar objects, and are not necessarily used to describe a specific order or sequence. It should be understood that the data used in this way can be interchanged where appropriate, so that the embodiments described herein can be implemented in an order other than that illustrated or described herein. In addition, the terms "including" or "having" and any variations thereof are intended to cover non-exclusive inclusions, for example, a process, method, system, product or device that includes a series of steps or units is not necessarily limited to those steps or units that are clearly listed, but may include other steps or units that are not clearly listed or inherent to these processes, methods, products or devices.

[0078] For ease of understanding, the specific process of the embodiment of the present invention is described below. Specifically, Figure 1 This is a flow chart of the two-phase flow gas kinetic algorithm optimization method of the present application, and the steps include:

[0079] S100, constructing a gas kinetic format algorithm for phase field solution, the gas kinetic format algorithm for phase field solution includes a continuous Boltzmann-BGK equation with source term;

[0080] S102, constructing the phase field interface flux of the gas kinetic format algorithm for phase field solution based on Chapman-Enskog analysis, wherein the phase field interface flux can be restored to the Cahn-Hilliard equation;

[0081] S104, constructing a gas kinetic format algorithm for solving the flow field, wherein the gas kinetic format algorithm for solving the flow field includes a continuous Boltzmann-BGK equation containing a source term;

[0082] S106, constructing a flow field interface flux of a gas kinetic format algorithm for solving the flow field based on Chapman-Enskog analysis, wherein the flow field interface flux can be restored to the Navier-Stokes equation;

[0083] S108. Based on the Cahn-Hilliard equation and the Navier-Stokes equation, a two-phase flow numerical algorithm is implemented by coupling the phase field and flow field.

[0084] It can be understood that the present invention starts from the Boltzmann-BGK equation and constructs a gas kinetic format algorithm for phase field and flow field coupling under a finite volume framework. Through the Chapman-Enskog expansion, it can accurately restore the numerical flux of the macroscopic equation, simplify the calculation process, and significantly improve the calculation efficiency and accuracy. This method can effectively simulate the interface dynamics and flow field evolution in complex two-phase flow systems.

[0085] The following is a detailed description of the application scheme in combination with a specific formula:

[0086] Since the method proposed in the present invention is a two-phase flow GKS method based on the phase field model, the phase field model is first introduced here, and then how to construct the numerical flux of the phase field and flow field in the GKS format is introduced, and the simulation of complex two-phase flow problems is realized by coupling the phase field and flow field. The technical solution of the present invention is described in detail below:

[0087] In the phase field theory of two-phase systems, the different fluid phases are characterized by the order parameter To identify, the order parameter is closely related to the free energy of the system. For an isothermal binary fluid system, the total free energy of the system consists of two parts: bulk free energy and interface free energy

[0088]

[0089] Where V is the control volume, and the bulk free energy can be expressed by the double potential well function ψ(φ) = β(φ-φ A ) 2 (φ-φ B ) 2 Indicates that κ and β are parameters related to surface tension σ and interface thickness W, respectively, and their expressions are as follows:

[0090]

[0091] Where φ A ,φ B are the gas phase and liquid phase order coefficients respectively, and the phase interface can be expressed as (φ A -φ B ) / 2. The chemical potential μ is defined as the variational derivative of the system free energy with respect to the order parameter,

[0092]

[0093] When the system is in equilibrium, the chemical potential is constant, especially for planar interfaces, the distribution of order coefficients can be described as:

[0094]

[0095] Where z is the coordinate perpendicular to the interface. The key to capturing the interface of two-phase fluids lies in the calculation of the order coefficient. The macroscopic control equation of the captured phase interface can be described by the Cahn-Hilliard equation.

[0096]

[0097] Where u is the fluid velocity and M is the migration coefficient. The Cahn-Hilliard equation satisfies the properties of mass conservation and energy not increasing with time. In incompressible two-phase flow, the flow field changes can be described by the Navier-Stokes equation:

[0098]

[0099] In two-phase flow, the density of the fluid can be calculated from the order parameter:

[0100]

[0101] where ρ 1 and ρ 2 is the density of the two-phase fluid, p is the hydrodynamic pressure, v is the kinematic viscosity, is the Navier-Stokes strain rate tensor of the velocity gradient, and ξ is the volume viscosity coefficient. total is the total force, including the surface tension F s There are different forms of surface tension, and here we use the form of potential function to minimize the generation of false velocities.

[0102]

[0103] Furthermore, we can redefine the pressure p h =p / ρ, and we get a non-conservative momentum equation based on velocity.

[0104]

[0105] In the mesoscopic method, the particle velocity distribution function is usually used to characterize the probability of particle velocity distribution at the corresponding time position, and its distribution satisfies the Boltzmann-BGK equation.

[0106]

[0107] Where f is the velocity distribution function, f eq is the equilibrium velocity distribution function, and τ is the relaxation time. According to the Chapman-Enskog analysis, the above-mentioned NS equation can be accurately restored. The technical solution proposed in the present invention starts from the BGK equation containing external force terms, and constructs numerical fluxes through Chapman-Enskog analysis, and the corresponding macro equations can be obtained, including the Cahn-Hilliard phase field equation and the NS equation, that is, the Low Mach Number Approximation (LMNA) method. The following will introduce how to construct the numerical flux.

[0108] Firstly, an aerodynamic scheme for phase field solution is constructed, starting from the continuous Boltzmann-BGK equation including source terms.

[0109] That is, a detailed introduction to steps S100-S102, specifically:

[0110] The gas kinetic format algorithm for phase field solution includes the continuous Boltzmann-BGK equation with source term based on the following formula:

[0111]

[0112] Where c is the velocity of the particle distribution function g, τ g is the relaxation time, and G is the phase field source term.

[0113] Furthermore, in order to accurately restore the Cahn-Hilliard equation, the equilibrium distribution function and the moments of the phase field source term satisfy the following conditions:

[0114]

[0115] Where η is the relaxation time adjustment coefficient, and the mobility can be expressed as parameter Set to 0, B α and C 0 is the coefficient to be determined. First, use Operator, we get the expression of the distribution function, which can be expressed as an equilibrium state and a source term using the Chapman-Enskog approximation.

[0116]

[0117] Then take the zero-order moment of the above formula and we get: Where F = <c α ,g> is the phase field interface flux of GKS. Substituting the phase field distribution function and the moments of the phase field source term into it, we can get the first-order phase field flux.

[0118]

[0119] The first-order phase field flux in the above equation is expanded by Chapman-Enskog, and the constructed first-order phase field flux expression can be restored to the Cahn-Hilliard equation. The GKS flux equation can be obtained by equation (13):

[0120]

[0121] Where n jk is the unit normal vector at the center jk of the interface between unit j and unit k. The distribution function is Taylor expanded at t = 0 with second-order time accuracy to obtain the function distribution function of position and time:

[0122]

[0123] And perform Chapman-Enskog expansion based on the expression of the distribution function

[0124]

[0125] Then, combined with the distribution function expression in the above formula, the following representation of the distribution function is obtained:

[0126]

[0127] By integrating the flux over time, we can finally get the flux expression of the phase field equation:

[0128]

[0129] where φ h and They are obtained by evolving the equations in half a time step, which is approximately

[0130] The time step Δt in GKS is given by the CFL condition:

[0131]

[0132] Where Δx is the grid step size.

[0133] Since the diffusion term in the CH (Cahn-Hilliard) equation contains the gradient of the chemical potential, it is necessary to calculate the fourth-order derivative of the order parameter. In order to characterize the influence of nearby units on the target unit in the diffusion term, this paper calculates this term using the following formula:

[0134]

[0135] The above process introduces in detail how to construct the interface flux of GKS. According to the flux formula, the GKS algorithm can be used to calculate the phase field part. The gas kinetic format algorithm based on the finite volume format framework is as follows:

[0136] (1) Determine the calculation area, divide or read the grid information, fluid parameters and control parameters.

[0137] (2) Based on the information read, an array is created and the corresponding parameters of the calculation area are initialized. After completion, a time step loop is entered. In the gas kinetic format, the macroscopic quantities will be directly evolved.

[0138] (3) Set boundary conditions and assign values ​​to the corresponding macro quantities at all boundaries in the flow field. This paper uses a virtual grid to process the boundaries of the flow field.

[0139] (4) Solve the flow field NS (Navier-Stokes) equations and first calculate the pressure and velocity (or momentum) of the flow field unit interface.

[0140] Flux, and then update the macroscopic quantities in the unit through the flux information of the interface. If there is external force, the external force is applied by splitting and the macroscopic quantities are updated at the same time.

[0141] (5) Determine whether the time step meets the loop condition, output data, and end the loop.

[0142] The above process constructs a flux based on the GKS format based on the Cahn-Hilliard phase field theory, and based on this flux, an interface capture GKS method based on phase field solution can be realized, in which the interface flux is given by the first-order Chapman-Enskog approximation, and the method has second-order time accuracy.

[0143] Next, we will start from the continuous Boltzmann-BGK equation and construct the numerical flux of the flow field with low Mach approximation in the finite volume framework. Through the Chapman-Enskog approximation expansion, we can restore to the Navier-Stokes equation. According to the assumption of constant density nearly incompressible flow, the momentum-based flux is simplified to a velocity-based flux. Different from the GKS method for dealing with high-speed compressible flow, by constructing a simpler interface flux for low-speed two-phase flow problems, the program calculation process can be simplified and the program's calculation efficiency can be improved. That is, the specific introduction of steps S104-S106: The expression of the continuous Boltzmann-BGK equation containing source terms for the flow field is:

[0144]

[0145] In the formula, τ f and S are the relaxation time of the flow field and the flow field source term respectively.

[0146] f eq =ρΓ(u,U) is the equilibrium distribution function, ρ is the density of the gas, ρ= <f>, momentum is ρu= <uf>.

[0147] The expression of Γ(u,U) is:

[0148]

[0149] In some examples, the flow source term can be expressed as:

[0150]

[0151] when When , the moment of the equilibrium distribution function of the flow field can be obtained as:

[0152]

[0153] The moment of the flow field source term is:

[0154]

[0155] According to the Chapman-Enskog analysis, the zero-order moment and first-order moment of the gas kinetic format algorithm equation for solving the flow field in the flow field are:

[0156]

[0157] The above moment calculation results show that the kinetic model can accurately recover the Navier-Stokes equations. In order to further derive the expression of the interface flux, the numerical flux of the flow field will be derived based on the low Mach approximation assumption.

[0158] Define a new distribution function as:

[0159]

[0160] where Γ 0 (u) is:

[0161]

[0162] The BGK equation expression defined according to the new distribution function is:

[0163]

[0164] The corresponding flow field equilibrium state and flow field source term expressions are:

[0165]

[0166] The moments of the equilibrium distribution function of the flow field are:

[0167]

[0168] Moments of the flow field source terms:

[0169]

[0170] Substitute the moments of the flow field equilibrium distribution function and the flow field source term into the zero-order moment to obtain the first-order flow field flux;

[0171] Substitute the moments of the flow field equilibrium distribution function and the flow field source term into the first-order moment to obtain the second-order flow field flux;

[0172] By performing Chapman-Enskog expansion on the first-order flow field flux and the second-order flow field flux, the flux formulas of pressure p and momentum ρU can be obtained:

[0173]

[0174] For constant density incompressible flow problems, a simplified velocity-based flux formula can be obtained:

[0175]

[0176] In real two-phase flow problems, the NS equations are usually used to describe the evolution of the flow field. How to couple the phase field and flow field in the framework of the gas kinetic format is a key issue. The following will introduce how to achieve phase field and flow field coupling in the GKS framework, that is, the specific introduction of step S108: the numerical algorithm steps of the gas kinetic format for simulating complex two-phase flow are given.

[0177] The interface flux formula in the phase field part is denoted by J φ According to formula (19),

[0178]

[0179] in, is the velocity of half a time step of the flow field evolution. The pressure and velocity interface flux formulas in the flow field are respectively denoted as J p and J u , according to formula (34), we can get:

[0180]

[0181] The gradient term in the source term can be calculated using the second-order isothermal center format:

[0182]

[0183] Where Φ is an arbitrary scalar and the coefficient a = θ r θ l ,b=θ r +θ l ,c=θ r -θ l ,θ r =(x j+1 -x j ) / δ x ,θ l =(x j -x j-1 ) / δ x Represents the scaling factors for the backward and forward step sizes, respectively.

[0184] It can be understood that the calculation process of the two-phase flow GKS algorithm in one time step is:

[0185] (1) After initializing the macroscopic quantity in each cell, the boundary conditions are applied to the calculation area, and the flux J at the cell interface is calculated according to formula (36) and formula (37) respectively: p and J u .

[0186] (2) In the phase field part, firstly according to the flux J in formula (37) u Calculate the velocity u of the flow field evolution for half a time step h , and then calculate the unit interface flux J according to formula (35) φ Then, according to formula (35), the macroscopic quantity in the unit at the next moment is calculated, that is, the update order parameter φ n+1 .

[0187] (3) In the flow field, the flux J calculated in (1) p and J u Evolving the flow field for one time step and applying the source term by the splitting method, the macroscopic quantity p of the flow field can be obtained. n+1 and u n+1 .

[0188] (4) The CH equation and the NS equation are solved by the GKS method, and the two-phase flow process of the coupled velocity field can be calculated.

[0189] Therefore, it can be summarized into the following steps:

[0190] Step 1: Determine the calculation area, divide or read the grid information, fluid parameters and control parameters, open an array based on the read information, initialize the corresponding parameters of the calculation area, enter the time step loop after completion, and enter Step 2.

[0191] Step 2: Set boundary conditions, assign values ​​to all boundary macro quantities in the flow field, and proceed to Step 3.

[0192] Step 3: Solve the Navier-Stokes equations of the flow field, calculate the pressure and velocity (or momentum) flux of the flow field unit interface according to the flow field flux formula, and then update the flow field macroscopic quantities in the unit through the flux information of the interface. If there is an external force, apply the external force through splitting, and update the macroscopic quantities, and enter Step 4.

[0193] Step 4: Calculate the macroscopic information of the half time step of the flow field flux evolution calculated in Step 3, calculate the phase field macroscopic quantity passing through the interface within the time step according to the phase field flux formula, that is, the order parameter, and update the macroscopic quantity in the unit, and enter Step 5.

[0194] Step5: Determine whether the time step meets the loop condition, if yes, go to Step2, if not, go to Step6.

[0195] Step 6: Output flow field data and end the simulation.

[0196] The grid information in the above Step 1 includes grid length, width, number of grids, unit grid length and grid coordinates; the fluid parameters include the two-phase density, velocity, pressure and viscosity parameters of each grid point; the control parameters include dimensionless parameters such as time step and CFL number.

[0197] The boundary conditions in the above Step 2 are similar to those in the macro method, and the boundary conditions are set directly by the macro values ​​at the given boundaries.

[0198] The Navier-Stokes equation in Step 3 above is:

[0199]

[0200] In the above formula, ρ is the fluid density, t is the time, u is the fluid velocity, p is the fluid dynamics pressure, S is the strain rate tensor, ξ is the volume viscosity coefficient, and F total It is the external force term acting on the fluid. The interface flux in this step includes the flux of pressure and velocity (or momentum). The unit interface flux can be used to update the macroscopic information in the unit.

[0201] The splitting method in Step 3 is implemented as follows:

[0202]

[0203] The speed u can be obtained by the first formula in the above solution. * , and then apply external force according to the splitting method to obtain the velocity u of the next time step n+1 .

[0204] The interface flux in Step 4 is the flux of the phase field order parameter. To calculate this flux, you need to first obtain the value of the half-time step of the flow field velocity evolution at the interface, and calculate it using the flow field flux formula in Step 3; then calculate it according to the phase field flux formula and update the order parameter information in the unit.

[0205] It is understandable that GKS, as a mesoscopic method under the finite volume framework, has the following main advantages: the flow field information is evolved and updated through macroscopic quantities, and the effective processing directly based on macroscopic quantities in traditional numerical methods can be directly applied, such as the reconstruction of physical quantities within the unit or at the interface, boundary conditions, and acceleration techniques for steady flow, while the kinematic format of the gas is constructed based on the continuous Boltzmann equation. There is no need to discretize and store the distribution function in velocity space, which greatly reduces the amount of calculation and storage of the method, which is far less than the cost of directly solving the BGK equation. The numerical flux under the finite volume framework is constructed by the continuous Boltzmann-BGK equation, and the flux constructed by the Chapman-Enskog expansion analysis method can be reduced to the macroscopic control equations (Cahn-Hilliard equations and Navier-Stokes equations), and the momentum-based flow field flux can be simplified to the velocity-based flow field flux based on the assumption of equal density and nearly incompressible flow. Once the form of the flux is obtained, a numerical algorithm based on the gas momentum format can be constructed. The algorithm is a direct evolution and update of the macroscopic quantities of the flow field. The calculation process is simple and the calculation accuracy is high. It can perform efficient numerical simulation of physical processes in low-speed flows.

[0206] Furthermore, the present invention constructs a numerical flux based on the Cahn-Hilliard phase field model in the GKS method, and successfully achieves accurate capture of the phase interface in the two-phase flow. The two-phase flow algorithm in the present invention significantly improves the interface capture capability when dealing with complex two-phase flow problems. The algorithm can accurately simulate complex phenomena such as phase separation and interface movement, overcoming the shortcomings of the traditional GKS method in interface capture.

[0207] Furthermore, the present invention achieves low Mach number adaptability. Specifically, by improving the traditional GKS method, the present invention constructs a flow field numerical flux that can accurately restore the macroscopic Navier-Stokes equation based on the low Mach number approximation assumption. By coupling the phase field flow field, an efficient and accurate two-phase flow algorithm is proposed. Compared with the traditional method, the present invention not only simplifies the construction process of the numerical flux, but also significantly improves the calculation efficiency and numerical stability, and can accurately simulate the complex phenomena of interface motion and flow field evolution in two-phase flow.

[0208] The technical solution of the present invention successfully overcomes the limitations of the traditional GKS method in dealing with variable density low-speed flow and two-phase flow problems through the above-mentioned technical innovations, and provides an efficient and accurate mesoscopic numerical algorithm. By adopting the low Mach number two-phase flow GKS method based on the phase field model proposed by the present invention, the interface dynamics and complex flow field evolution in the two-phase flow system can be efficiently and accurately simulated.

[0209] In order to verify the correctness of the algorithm proposed in the present invention, the following example is performed.

[0210] In the simulation of the stationary droplet example, a circular droplet with a radius of R is placed in the center of an L×L calculation area, and the boundary conditions around it are set to periodic boundaries. The distribution of the order parameter is:

[0211]

[0212] Among them, (x c ,y c ) is the position coordinate of the center of the calculation area (L / 2, L / 2), the calculation grid is 200×200, and the droplet radius is R=50. The other parameters of the model are set to ρ A =500,ρ B =1, v=0.2, W=4.0, σ=1.0×10 -5 ,M=0.5,CFL=0.5. Order parameter distribution of stationary droplet in equilibrium state. Figure 2 It can be seen that the density distribution along the horizontal center line is in good agreement with the standard solution. When the system reaches equilibrium, the relationship between the droplet radius R and the pressure difference ΔP inside and outside the droplet satisfies the description of Laplace's law, that is, ΔP = σ / R, and the pressure P is calculated as in Figure 3 The fitting curves of the pressure difference between the inside and outside of the droplet and the inverse of the radius are given under three different surface tension conditions. The results show that the model conforms to Laplace's law.

[0213] In the simulation of Rayleigh-Taylor instability, the computational domain is a cube with a square cross-section and an aspect ratio of 4:1. The grid of the computational area is 128×128×512. No-slip rebound boundary conditions are applied to the upper and lower boundaries, and periodic boundary conditions are applied to the four side boundaries. Other parameters are set as: W=4,σ=1×10 -4 , CFL = 0.1, At = 0.5. The Peclet number is fixed at 744. The initial interface function of the two-phase fluid and the order parameter distribution at the initial moment are as follows:

[0214]

[0215] The two dimensionless parameters closely related to this problem are the Atwood number and the Reynolds number, which are defined as follows:

[0216]

[0217] Among them, g is the acceleration due to gravity, and the other parameters are fixed to L 0 =256, W=5,σ=5×10 -5 , CFL = 0.25, the characteristic time is Defined as The Péclet number of . Figure 4 The interface evolution diagrams with Reynolds number of 1024 at different times are given. Figure 5 The time variations of the positions of the bubble front, saddle point and spike are given for a Reynolds number of 1024, which are also in good agreement with the reference results. In summary, it can be proved that this method has high accuracy and reliability in simulating complex two-phase flow problems.

[0218] According to another aspect of the present application, a two-phase flow gas kinetic algorithm optimization device is also disclosed, specifically, as follows Figure 6 As shown, the device comprises:

[0219] The first building block is used to build a gas kinetic format algorithm for phase field solution, which includes a continuous Boltzmann-BGK equation with a source term.

[0220] The first analysis module is used to construct the phase field interface flux of the gas kinetic format algorithm for phase field solution based on Chapman-Enskog analysis, wherein the phase field interface flux can be restored to the Cahn-Hilliard equation;

[0221] The second building block is used to build a gas kinetic format algorithm for flow field solution, and the gas kinetic format algorithm for flow field solution includes a continuous Boltzmann-BGK equation with a source term;

[0222] The second analysis module is used to construct the flow field interface flux of the kinetic format algorithm for flow gas solution based on Chapman-Enskog analysis, wherein the flow field interface flux can be restored to the Navier-Stokes equation;

[0223] The coupling module is used to realize the two-phase flow numerical algorithm based on the Cahn-Hilliard equation and the Navier-Stokes equation by coupling the phase field and flow field.

[0224] The application introduction of the relevant modules of the device in this example can refer to the relevant introduction of the principle of the above method, which will not be repeated here.

[0225] According to another aspect of the present application, the present application also discloses an electronic device, which includes a memory and at least one processor, wherein instructions are stored in the memory; at least one processor calls the instructions in the memory so that the electronic device executes each step of the above-mentioned two-phase flow gas kinetic algorithm optimization method.

[0226] above Figure 6 The apparatus in the embodiment of the present invention is described in detail from the perspective of modular functional entities, and the electronic device in the embodiment of the present invention is described in detail from the perspective of hardware processing.

[0227] Figure 7 7 is a schematic diagram of the structure of an electronic device provided by an embodiment of the present invention. The electronic device 700 may have relatively large differences due to different configurations or performances, and may include one or more processors (central processing units, CPU) 710 (for example, one or more processors) and a memory 720, and one or more storage media 730 (for example, one or more mass storage devices) storing application programs 733 or data 732. Among them, the memory 720 and the storage medium 730 can be temporary storage or permanent storage. The program stored in the storage medium 730 may include one or more modules (not shown in the figure), and each module may include a series of instruction operations in the electronic device 700. Furthermore, the processor 710 can be configured to communicate with the storage medium 730 to execute a series of instruction operations in the storage medium 730 on the electronic device 700.

[0228] The electronic device 700 may also include one or more power supplies 740, one or more wired or wireless network interfaces 750, one or more input and output interfaces 750, and / or one or more operating systems 731, such as Windows Serve, Mac OS X, Unix, Linux, FreeBSD, etc. It will be appreciated by those skilled in the art that Figure 7 The structure of the electronic device shown does not constitute a limitation on the electronic device, and may include more or less components than shown in the figure, or combine some components, or arrange the components differently.

[0229] The present invention also provides a computer-readable storage medium, which may be a non-volatile computer-readable storage medium or a volatile computer-readable storage medium. Instructions are stored in the computer-readable storage medium. When the instructions are executed on a computer, the computer executes the steps of the logistics order control method.

[0230] Those skilled in the art can clearly understand that, for the convenience and brevity of description, the specific working process of the system, device or unit described above can refer to the corresponding process in the aforementioned method embodiment and will not be repeated here.

[0231] If the integrated unit is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or all or part of the technical solution can be embodied in the form of a software product, which is stored in a storage medium and includes several instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) to perform all or part of the steps of the methods of each embodiment of the present invention. The aforementioned storage medium includes: U disk, mobile hard disk, read-only memory (ROM), random access memory (RAM), disk or optical disk, etc., various media that can store program codes.

[0232] The above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit the same. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that the technical solutions described in the aforementioned embodiments may still be modified, or some of the technical features thereof may be replaced by equivalents. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.< / uf> < / f> < / f>

Claims

1. A two-phase gas kinetic format algorithm optimization method, characterized in that: The method comprises: Constructing a gas kinetic format algorithm for phase field solution, wherein the gas kinetic format algorithm for phase field solution includes a continuous Boltzmann-BGK equation with a source term; Constructing the phase field interface flux of the gas kinetic format algorithm for the phase field solution based on Chapman-Enskog analysis, wherein the phase field interface flux can be restored to the Cahn-Hilliard phase field equation; Constructing a gas kinetic format algorithm for flow field solution, wherein the gas kinetic format algorithm for flow field solution includes a continuous Boltzmann-BGK equation containing a source term; Constructing the flow field interface flux of the gas kinetic format algorithm for solving the flow field based on Chapman-Enskog analysis, wherein the flow field interface flux can be restored to the Navier-Stokes equation; Based on the Cahn-Hilliard equation and the Navier-Stokes equation, a two-phase flow numerical algorithm is implemented by coupling the phase field and flow field.

2. The two-phase flow gas kinetic algorithm optimization method according to claim 1 is characterized in that: The gas kinetic format algorithm for phase field solution is constructed, and the gas kinetic format algorithm for phase field solution includes a continuous Boltzmann-BGK equation containing a source term based on the following formula: Where c is the velocity of the particle distribution function g, τ g is the relaxation time, and G is the phase field source term.

3. The two-phase flow gas kinetic algorithm optimization method according to claim 1 is characterized in that: The phase field interface flux of the gas kinetic format algorithm for phase field solution is constructed based on Chapman-Enskog analysis, wherein the phase field interface flux can be restored to the Cahn-Hilliard phase field equation including: Obtaining the phase field equilibrium distribution function and phase field source term of each order moment satisfying the first condition; The equilibrium distribution function and the moments of the source term satisfy the following conditions: Where η is the relaxation time adjustment coefficient, and the mobility can be expressed as parameter Set to 0, B α and C0 are coefficients to be determined. First, use operator, to obtain the expression of the phase field equilibrium distribution function; The phase field equilibrium distribution function is expressed based on Chapman-Enskog: Taking the zero-order moment of the phase field equilibrium distribution function to obtain a zero-order flux equation; The zero-order flux equation is: In the formula, F = <c α , g>, F is the phase field interface flux of the gas kinetic format algorithm for phase field solution; Substituting the phase field equilibrium distribution function and each order moment of the phase field source term into the zero-order flux equation to obtain the first-order phase field flux; The first-order phase field flux is: The first-order phase field flux is subjected to a Chapman-Enskog expansion, and the phase field interface flux is restored to the Cahn-Hilliard phase field equation.

4. The two-phase flow gas kinetic algorithm optimization method according to claim 1 is characterized in that: The gas kinetic format algorithm for constructing the flow field solution includes a continuous Boltzmann-BGK equation containing a source term based on the following formula: In the formula, τ f is the relaxation time, S is the flow source term, c is the velocity of the particle distribution function f, f eq =ρΓ(u, U) is the equilibrium distribution function; ρ is the density of the gas, ρ = <f> ；< / f> The expression of Γ(u, U) is:

5. The two-phase flow gas kinetic algorithm optimization method according to claim 1 is characterized in that: The flow field source term is expressed as: The flow field interface flux of the gas kinetic format algorithm for solving the flow field is constructed based on Chapman-Enskog analysis, wherein the flow field interface flux is restored to the Navier-Stokes equations including: Sure The moment of the equilibrium distribution function of the flow field and the moment of the flow field source term at : Among them, the moment of the flow field equilibrium distribution function is: The moment of the flow field source term is: The zero-order moment and the first-order moment of the gas kinetic format algorithm for solving the flow field obtained based on the Chapman-Enskog analysis include: Among them, the zero-order moment of the gas kinetic format algorithm for solving the flow field is: The first-order moment of the gas kinetic format algorithm for solving the flow field is: Define a new distribution function: in, Based on the new distribution function, the Boltzmann-BGK expression, the moments of the equilibrium distribution function of the flow field and the moments of the flow field source term are determined. The Boltzmann-BGK expression is: The moments of the equilibrium distribution function of the flow field are: The moments of the flow field source terms are: Substituting each order moment of the flow field equilibrium distribution function and the flow field source term into the zero-order moment to obtain the first-order flow field flux; Substituting the moments of the flow field equilibrium distribution function and the flow field source term into the first-order moment to obtain the second-order flow field flux; The first-order flow field flux and the second-order flow field flux are subjected to Chapman-Enskog expansion, and the flow field interface flux is reduced to the Navier-Stokes equation.

6. A two-phase gas kinetic algorithm optimization device, characterized in that: The device comprises: A first building module is used to build a gas kinetic format algorithm for phase field solution, wherein the gas kinetic format algorithm for phase field solution includes a continuous Boltzmann-BGK equation containing a source term; A first analysis module is used to construct a phase field interface flux of the gas kinetic format algorithm for phase field solution based on Chapman-Enskog analysis, wherein the phase field interface flux can be restored to the Cahn-Hilliard phase field equation; A second construction module is used to construct a gas kinetic format algorithm for flow field solution, wherein the gas kinetic format algorithm for flow field solution includes a continuous Boltzmann-BGK equation containing a source term; A second analysis module is used to construct a flow field interface flux of the gas kinetic format algorithm for solving the flow field based on Chapman-Enskog analysis, wherein the flow field interface flux can be restored to the Navier-Stokes equation; The coupling module is used to realize a two-phase flow numerical algorithm by coupling a phase field and a flow field based on the Cahn-Hilliard equation and the Navier-Stokes equation.

7. An electronic device, characterized in that: The electronic device includes a memory and at least one processor, wherein the memory stores instructions; the at least one processor calls the instructions in the memory so that the electronic device executes each step of the two-phase flow gas kinetic algorithm optimization method as claimed in claims 1-5.

8. A computer-readable storage medium having instructions stored thereon, characterized in that: When the instructions are executed by the processor, the various steps of the two-phase flow gas kinetic algorithm optimization method as described in any one of claims 1-5 are implemented.

Images