Simulation method and system for large-density-ratio two-phase flow
By employing irregular nodes and the moving least squares approximation method in two-phase flow simulation, combined with the lattice Boltzmann method, the problems of complex geometries and coupling of time and space steps in traditional grid methods are solved, achieving high-precision and stable two-phase flow simulation.
Patent Information
- Application Number
- CN202511602307.3
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-11-04
- Publication Date
- 2026-02-10
AI Technical Summary
Traditional mesh methods suffer from approximation errors when dealing with complex geometries, cannot effectively handle the coupling between time and space steps, and struggle to maintain high accuracy and computational stability, especially in simulations of two-phase flows with high density ratios.
By employing irregular nodes and moving least squares approximation methods, combined with the lattice Boltzmann method, high-precision simulation of the flow field and phase field is achieved by performing collision operations and migrations of the flow field and phase field distribution functions on irregular nodes and using MLS approximation to reconstruct the distribution functions.
It breaks through the limitations of traditional mesh methods, can flexibly handle complex geometries, improve simulation adaptability, accuracy and stability, and maintains computational efficiency and accuracy at large time steps. It is suitable for numerical simulation of multiphase flow and complex flow fields.
Smart Images
Figure CN121503320A_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of two-phase flow simulation technology, specifically relating to a simulation method and system for two-phase flow with a high density ratio. Background Technology
[0002] Two-phase flow is ubiquitous in nature and engineering applications, such as the surface tension of liquid droplets, natural gas extraction and transportation, sulfidation, absorption, evaporation, condensation and chemical reactions in chemical processes, interface treatment at the melt front in injection molding, and the movement of bubbles in blood vessels. Numerical simulation studies of two-phase flow are of great significance for improving industrial equipment, increasing production efficiency, enhancing product quality, and developing new products. Due to the coupling of multiple physical mechanisms and complex topological changes at the phase interfaces, the dynamic behavior of this type of flow is extremely complex, making it a very challenging problem in fluid mechanics research. Traditional theoretical analysis methods struggle to obtain analytical solutions; while experimental methods face problems such as high cost and limited measurement accuracy. Therefore, developing efficient and reliable numerical simulation methods is crucial for promoting the research and application of two-phase flow.
[0003] Traditional solvers based on the macroscopic Navier-Stokes equations have been the most commonly used methods for two-phase flow simulations since the 1960s. In this field, researchers have developed various numerical methods, such as finite difference, finite volume, and finite element methods, to discretize and solve the nonlinear Navier-Stokes equations. These methods have shown significant advantages in solving many complex problems in computational fluid dynamics. However, they all struggle to handle complex nonlinear convection terms and numerical discretization problems involving viscous terms with higher-order derivatives. Furthermore, solving incompressible flows typically requires additional calculation of the pressure Poisson equation, which is often slow, thus significantly impacting computational efficiency.
[0004] The Lattice Boltzmann Method (LBM or LB method) is a computational fluid dynamics method based on statistical physics principles. Its core idea is to transform the fluid flow problem in a continuous medium into a particle motion problem in discrete space. In recent years, it has been widely used for two-phase flow simulations. As a numerical method originating from the mesoscopic scale, LBM differs from the direct solution of the traditional Navier-Stokes equations by indirectly obtaining macroscopic fluid variables through the evolution of particle distribution functions. Its method has a simple structure, is easy to implement in parallel, and can easily achieve complex boundary conditions, making it a powerful tool for simulating multiphase fluid dynamics. Through the collaborative efforts of many scholars, LBM has achieved great success in areas such as two-phase flow system modeling, fluid-structure interaction modeling, and nonlinear equation system modeling, and has become an important method for solving fluid mechanics problems.
[0005] While LBM has significant advantages in two-phase flow systems, it still has some limitations. Especially when dealing with complex solution domains or irregular solids or boundaries, traditional LBM struggles to realize its advantages. Obtaining high-precision results requires finer meshes or larger solution ranges, which significantly increases computational resources and complexity, posing a challenge for LBMs using regular lattice structures. Furthermore, given the coupling between time steps and mesh spacing, LBM's computational capabilities at large time steps are relatively weak. Summary of the Invention
[0006] The purpose of this invention is to address the problems of traditional mesh methods in handling complex geometries, such as approximation errors in geometry, inability to handle the coupling between time and space steps, and inability to maintain high accuracy and computational stability for complex flows and geometries. This invention proposes a simulation method and system for two-phase flows with a high density ratio.
[0007] The technical solution of the present invention is as follows: Firstly, a method for simulating two-phase flow with a high density ratio, comprising the following steps: Generate irregular nodes within the computational region and initialize them; Based on irregular nodes, collision operations are performed on the distribution functions of the flow field and phase field; At irregular nodes, the distribution functions of the flow field and phase field are traced in reverse along their characteristic directions using the semi-Lagrange method to obtain the starting point; Based on the starting point, the distribution functions of the flow field and phase field at the migration time step at the starting point are reconstructed by interpolation using the moving least squares approximation, so as to obtain the distribution functions of the flow field and phase field after MLS approximation. The migration step is performed based on the distribution functions of the flow field and the phase field after the MLS approximation, and the computation nodes are updated to obtain the distribution functions of the flow field and the phase field after migration. Using the distribution function of the flow field and the distribution function of the phase field after migration, the physical quantities at the current time step are calculated; The equilibrium distribution functions and source term distribution functions of the flow field and phase field after migration are calculated, and the collision migration operation is repeatedly performed on the distribution functions of the flow field and phase field to calculate the physical quantities at each time step until the preset number of simulation time steps or convergence criteria are met, thus completing the simulation of a two-phase flow with a high density ratio.
[0008] Preferably, the collision migration operation on the distribution functions of the flow field and phase field is implemented by the MRT model; the dynamic behavior of the flow field is described by solving the Navier-Stokes equations using the lattice Boltzmann method, and the interface evolution of the phase field is described by solving the Allen-Cahn equations using the lattice Boltzmann method. The MRT evolution of the MRT-LB model of the Navier-Stokes equations is as follows:
[0009] in, The distribution function representing the flow field. Indicates the time step. Indicates the position of irregular nodes. Indicates time, Indicates the starting point. Represents the Kronecker symbol, when hour ,when hour , Represents the invertible collision matrix The Middle Line number The elements of the column are , Denotes a diagonal relaxation matrix. Represents the transformation matrix. The source term distribution function represents the flow field. The equilibrium distribution function representing the flow field; The MRT evolution of the conservative ACE-LB model of the Allen-Cahn equation is as follows:
[0010] in, Indicates the order parameter The distribution function, Represents the invertible collision matrix The Middle Line number The elements of the column are , Denotes a diagonal relaxation matrix. The source term distribution function represents the phase field. This represents the equilibrium distribution function of the phase field.
[0011] As a preferred option, the specific formula for calculating the starting point is:
[0012] in, Indicates the starting point. Indicates the position of irregular nodes. Represents discrete velocity. Indicates the time step.
[0013] Preferably, based on the starting point, the distribution functions of the flow field and phase field at the migration time step at the starting point are reconstructed by interpolation using the moving least squares approximation to obtain the distribution functions of the flow field and phase field after MLS approximation. Specifically, this includes the following steps: Define the generalized distribution function The distribution function representing the flow field or phase field, the generalized distribution function Defined in the region middle, Indicates that it includes all The space of dimensional vectors, the generalized distribution function The approximate function is :
[0014] in, Denotes basis functions. Denotes the coefficients to be determined. Indicates the number of basis functions. Denotes the offset polynomial basis functions, superscript To represent transpose, we have:
[0015] in, Indicates a reference point. and Represents irregular nodes exist direction and directional components, and Indicates reference point In direction and Component of direction; The support domain of the coefficients to be determined is selected, and a weight function corresponding to each meshless discrete node and its influence radius is constructed within the support domain using cubic splines; wherein, the support domain of the coefficients to be determined includes A gridless discrete node ; Based on the weighting function, the coefficients to be determined are determined by minimizing the weighted squared error of all meshless discrete nodes in the support domain:
[0016] in, For Gram matrices, the matrix is... ,matrix , Indicates the corresponding to the first The weight function of a gridless discrete node and its influence radius, matrix superscript Indicates transpose. Indicates the first The generalized distribution function corresponding to a gridless discrete node; Substitute the coefficients to be determined into the approximate function. ,get:
[0017] in, Indicates the number of coefficients to be determined; make The generalized distribution functions of the flow field and phase field after MLS approximation are obtained as follows:
[0018] in, , This represents the flow field distribution function after MLS approximation. This represents the phase field distribution function after the MLS approximation. Representing shape functions The transpose of .
[0019] Preferably, a migration step is performed based on the distribution functions of the flow field and the phase field after the MLS approximation to obtain the migration distribution functions of the flow field and the phase field, specifically as follows: At the starting point The flow field and phase field migration steps are performed above, and the distribution function of the migrated flow field is:
[0020] in, The distribution function represents the flow field after migration. Indicates the starting point The flow field distribution function after MLS approximation; The distribution function of the phase field after migration is:
[0021] in, The distribution function of the phase field after migration. Indicates the starting point The phase field distribution function after MLS approximation.
[0022] Preferably, the physical quantities include order parameters, fluid density, kinematic viscosity coefficient, chemical potential, external force term, flow field velocity, and hydrodynamic pressure.
[0023] As a preferred option, the order parameter The calculation formula is:
[0024] in, The distribution function of the phase field after migration; fluid density The calculation formula is:
[0025] in, Indicates fluid The order parameter, Indicates fluid The order parameter, Indicates fluid density, Indicates fluid density, kinematic viscosity coefficient The calculation formula is:
[0026] in, Indicates fluid The kinematic viscosity coefficient, Indicates fluid The kinematic viscosity coefficient; Chemical potential The calculation formula is:
[0027] in, Indicates the order parameter The average value, , Indicates the order parameter gradient, and Indicates the thickness of the interface and surface tension The relevant physical constants are specifically defined as follows:
[0028] External force The calculation formula is:
[0029] in, Indicates surface tension. Represents volume force; The formula for calculating the flow field velocity is:
[0030] in, Indicates the flow field velocity. Represents discrete velocity. The distribution function represents the flow field after migration. Indicates the time step; Dynamic water pressure The calculation formula is:
[0031] in, Indicates the speed of sound. This represents the weight at position 0. The gradient representing density. As an auxiliary variable, , This indicates the velocity at position 0.
[0032] Preferably, the equilibrium distribution function of the phase field after migration for:
[0033] in, Indicates the weighting coefficient. Indicates free parameters, Indicates the order parameter, Indicates the speed of sound. Represents discrete velocity. Indicates the flow field velocity. ; Equilibrium distribution function of the flow field after migration for:
[0034] in, Indicates dynamic water pressure. Indicates fluid density, As an auxiliary variable, , Represents discrete velocity; Source term distribution function of the phase field after migration for:
[0035] in, For the order parameter The function has , Indicates fluid The order parameter, Indicates fluid The order parameter, Indicates the interface thickness. Represents the unit normal vector. The partial derivative with respect to time t; Source term distribution function of the flow field after migration for:
[0036] in, Indicates the external force term. This indicates the density of the liquid phase. This represents the density of the gas phase. Indicates the order parameter The gradient.
[0037] The beneficial effects of this invention are: 1. This invention extends LBM to the meshless case, breaking through the limitations of the traditional uniform mesh method, and can flexibly handle flow problems with complex geometries, improving the adaptability and accuracy of the simulation.
[0038] 2. This invention can effectively eliminate the coupling problem between time step and spatial step, thereby achieving a stable calculation process and maintaining good stability even at a large time step, thus improving calculation efficiency and accuracy.
[0039] 3. This invention inherits all the advantages of the traditional LBM and gives full play to the advantages of the LB model in interface capture and parallel computing, making it particularly suitable for numerical simulation of multiphase flow and complex flow fields.
[0040] 4. This invention employs a non-uniform spatial discretization method, which successfully solves the two-phase flow problem in complex geometric domains while maintaining high accuracy and stability, overcoming the shortcomings of traditional methods in handling complex geometric regions.
[0041] Secondly, a simulation system based on a high density ratio two-phase flow includes: The first module is used to generate and initialize irregular nodes within the computational region; The second module is used to perform collision operations on the distribution functions of the flow field and phase field based on irregular nodes; The third module is used to trace the distribution functions of the flow field and phase field in reverse along their characteristic directions at irregular nodes based on the semi-Lagrange method to obtain the starting point; The fourth module is used to interpolate and reconstruct the distribution functions of the flow field and phase field at the starting point using the moving least squares approximation, based on the starting point, to obtain the distribution functions of the flow field and phase field after the MLS approximation. The fifth module performs a migration step based on the distribution functions of the flow field and the phase field after the MLS approximation, and updates the computation nodes to obtain the distribution functions of the flow field and the phase field after migration. The sixth module is used to calculate the physical quantities at the current time step using the distribution function of the flow field and the distribution function of the phase field after migration. The seventh module is used to calculate the equilibrium distribution function and source term distribution function of the flow field and phase field after migration, and repeatedly perform collision migration operation on the distribution functions of the flow field and phase field to calculate the physical quantities at each time step until the preset number of simulation time steps or convergence criteria are met, thus completing the simulation of high density ratio two-phase flow.
[0042] Thirdly, a non-transitory computer-readable storage medium is provided that stores computer instructions for causing a computer to perform the method as described in the first aspect. Attached Figure Description
[0043] Figure 1 The diagram shows a flowchart of a simulation method for a two-phase flow with a high density ratio.
[0044] Figure 2 The starting point is shown in the standard grid. and arrival point The positions overlap.
[0045] Figure 3 The image shows an unstructured mesh composed of discrete points (solid gray dots), where the starting point... Not necessarily related to the destination coincide.
[0046] Figure 4 The image shows the shape of the droplet before and after deformation.
[0047] Figure 5 The image shows the phase separation results at different times within the regular region.
[0048] Figure 6 The image shows the phase separation results at different times within the irregular region.
[0049] Figure 7 The figure shows the Reynolds number. Re =20, Weber number We =8000 Droplet impacts on the thin film at different times.
[0050] Figure 8 The figure shows the Reynolds number. Re =500, Weber number We =8000 Droplet impacts on the thin film at different times.
[0051] Figure 9 The figure shows the Reynolds number. Re =1000, Weiber number We =8000 Droplet impacts on the thin film at different times.
[0052] Figure 10 The figure shows the diffusion radius at different Reynolds numbers.
[0053] Figure 11 The image shows droplets passing through a cylindrical array, with a density ratio of 500 and a capillary number of... Ca =0.1 Interface evolution at different times.
[0054] Figure 12 The figure shows the evolution of the bubble rising interface morphology in a periodically contracting-expanding capillary at different times.
[0055] Figure 13 The figure shows the deformation coefficient of the droplet (the distance the droplet's center of mass moves) for different droplet sizes.
[0056] Figure 14 The figure shows the average axial velocity of the droplets (the distance moved relative to the center of mass of the droplet) for different droplet sizes. Detailed Implementation
[0057] Exemplary embodiments of the present invention will now be described in detail with reference to the accompanying drawings. It should be understood that the embodiments shown and described in the drawings are merely exemplary and are intended to illustrate the principles and spirit of the invention, and are not intended to limit the scope of the invention.
[0058] Example 1: This invention provides a simulation method for high density ratio two-phase flow, involving two LB equations: one for solving the Navier-Stokes equations for the flow field, and the other for solving the Allen-Cahn equations (ACE) for interface evolution. In this method, the collision step still occurs at Eulerian points, consistent with the traditional LBM method. However, due to the use of a meshless point distribution, the Eulerian points do not always coincide with the migration target location of the distribution function, thus preventing the migration step from being directly completed as in standard LBM. To address this issue, this invention performs the migration step within the SL framework. Specifically, at each time step, the starting point of the distribution function is traced back along the characteristic line, and then the distribution function from the previous time step is interpolated and reconstructed using the moving least squares method to obtain the updated value for the current time step. Since the propagation velocity in each direction is constant, the tracing of the characteristic line only needs to be performed once, and the starting point can be accurately calculated. The advantage of this method is that it does not require the starting point to strictly coincide with the Eulerian point, thus providing greater freedom for time and space discretization.
[0059] like Figure 1 As shown, a simulation method for a two-phase flow with a high density ratio includes the following steps: S1. Generate irregular nodes within the computational region and initialize them; S2. Based on irregular nodes, perform collision operations on the distribution functions of the flow field and phase field; S3. At irregular nodes, the distribution functions of the flow field and phase field are traced in reverse along their characteristic directions using the semi-Lagrange method to obtain the starting point; S4. Based on the starting point, the distribution functions of the flow field and phase field at the migration time step at the starting point are reconstructed by interpolation using the moving least squares approximation to obtain the distribution functions of the flow field and phase field after MLS approximation. S5. Based on the distribution functions of the flow field and the phase field after the MLS approximation, a migration step is performed, and the computation nodes are updated to obtain the distribution functions of the flow field and the phase field after migration. S6. Calculate the physical quantities at the current time step using the distribution function of the flow field and the distribution function of the phase field after migration; S7. Calculate the equilibrium distribution function and source term distribution function of the flow field and phase field after migration, and repeat the collision migration operation on the distribution functions of the flow field and phase field to calculate the physical quantities at each time step until the preset number of simulation time steps or convergence criteria are met, thus completing the simulation of high density ratio two-phase flow.
[0060] In this embodiment, the D2Q9 LB model of the MRT method is used to study the two-dimensional hydrodynamic problem and solve the conservative ACE within the phase-field theory framework. Under the LB framework of this ACE, the MRT evolution of the conservative ACE-LB model can be expressed as:
[0061] Similar to the conservative ACE LB model, the MRT collision migration operation of the incompressible Navier-Stokes equations MRT-LB model can be expressed as:
[0062] in, Indicates the order parameter The distribution function, Indicates the time step. Indicates the position of irregular nodes. Indicates time, Indicates the starting point. Represents the Kronecker symbol, when hour ,when hour ; Represents the invertible collision matrix The Middle Line number The elements of the column are , Denotes a diagonal relaxation matrix. The distribution function representing the flow field. Represents the invertible collision matrix The Middle Line number The elements of the column are , Denotes a diagonal relaxation matrix. The transformation matrix, in the D2Q9 lattice model, is represented as:
[0063] and Let the equilibrium distribution function be:
[0064]
[0065] in, Indicates the weighting coefficient. Indicates free parameters, Indicates the flow field velocity. Indicates the speed of sound. Indicates dynamic water pressure. Indicates the flow field density. As an auxiliary variable, , Represents discrete velocity; in the D2Q9 model, the weighting coefficients It is given in the following manner: speed of sound Discrete velocity Defined as:
[0066] in, Indicates lattice velocity, , Indicates grid spacing. Indicates the time step. Represents the cosine function. Represents pi (π). This represents the sine function.
[0067] and The source term distribution functions representing the phase field and flow field are:
[0068]
[0069] in, Represents the unit normal vector. This represents the partial derivative with respect to time t. This indicates the density of the liquid phase. This represents the density of the gas phase. Indicates the order parameter The gradient.
[0070] Furthermore, Chapman-Enskog analysis showed that mobility With relaxation time and free parameters related:
[0071] in, This represents the relaxation time of the phase field. Kinematic viscosity of fluids Defined as:
[0072] in, This represents the relaxation time of the flow field.
[0073] In this embodiment, the formula for calculating the starting point is: .
[0074] In this embodiment, the meshless method is characterized by its complete independence from the mesh, making it particularly suitable for handling large deformations and dynamic node problems. It employs a striped sparse matrix structure, resulting in high computational efficiency and facilitating the solution of large-scale problems. Furthermore, it eliminates the need for re-meshing, allowing for adaptive point addition during computation, thus offering greater flexibility and smoother, more continuous results. Unlike traditional methods, the core of the meshless method lies in the construction of shape functions. Since the MRT-LB models of the flow field and phase field are similar in format and identical in the MLS approximation, a generalized distribution function is introduced for simplicity. When describing the flow field, When describing the phase field, it is The distribution functions of the flow field and phase field at the starting point during the migration time step are reconstructed using the moving least squares approximation, resulting in the distribution functions of the flow field and phase field after the MLS approximation. This process includes the following steps: Define scalar functions The distribution function representing the flow field or phase field is defined in the region. In, its approximate value is At point MLS approximation at [location] for:
[0075] in, Denotes basis functions. Denotes the coefficients to be determined. This represents the number of basis functions, whose coefficients depend on... This demonstrates the "shifting" characteristic of the MLS method. To improve stability, offset polynomial basis functions are employed. , Describing the basis functions of the offset polynomial The transpose of is:
[0076] The coefficients are in a selected subdomain As determined in the text, this subdomain is called The supporting domain. This subdomain contains A gridless discrete node its functions The value is known: Its corresponding matrix form is:
[0077] in, This indicates the number of coefficients to be determined.
[0078] In the MLS framework, coefficients By minimizing the norm What we got was... Indicates the corresponding to the first Scalar function of a gridless discrete node Approximate value, inner product for:
[0079] in, Represents the weight function. , Represents a node The radius of influence. Weight function. Corresponding to node and its radius of influence In this embodiment, cubic splines are mainly used for construction, and the specific function definition is as follows:
[0080] weight function Through formula Calculations are performed to aid understanding. The goal of the optimization process is to minimize the following function:
[0081] This function effectively measures the weighted squared error between nodes. To achieve the minimization of the stationarity condition, the following system of equations is obtained:
[0082] The corresponding matrix form can be written as:
[0083]
[0084]
[0085]
[0086] matrix It is called a Gram matrix, and That is In the The projection of Zhang Cheng onto the function space (i.e., the approximate space). If the Gram matrix is nonsingular, then the coefficients... It can be explicitly obtained using the following formula:
[0087] Will Substitution In the middle, we get:
[0088] Similarly, we have:
[0089] in, It is a shape function, defined as:
[0090] In summary, it can be used The distribution functions and related physical information of the phase field and flow field after least squares approximation are calculated.
[0091] In this embodiment, the physical quantities include order parameters, fluid density, kinematic viscosity coefficient, chemical potential, external force term, flow field velocity, and hydrodynamic pressure, specifically: Calculate the migration steps and update the nodes. Specifically: Obtain the flow field distribution function after MLS approximation and phase field distribution function :
[0092]
[0093] in, and Let represent the flow field distribution function and the phase field distribution function after the collision, respectively; At the starting point The flow field and phase field migration steps are performed above, and the distribution function of the migrated flow field is:
[0094] in, The distribution function represents the flow field after migration. Indicates the starting point The flow field distribution function after MLS approximation; The distribution function of the phase field after migration is:
[0095] in, The distribution function of the phase field after migration. Indicates the starting point The phase field distribution function after MLS approximation.
[0096] Calculate the order parameter The specific formula is as follows:
[0097] in, This represents the distribution function of the phase field after migration. Since the distribution function after migration at the current moment is the same as the distribution function used for collision at the next moment, the time is not specified.
[0098] The specific formulas for calculating fluid density and kinematic viscosity are as follows:
[0099]
[0100] in, Indicates fluid density, Indicates fluid The order parameter is usually set to , Indicates fluid The order parameter is usually set to 0, therefore, when At that time, the interface is located at this position. Indicates fluid density, Indicates fluid density, Indicates the kinematic viscosity coefficient. Indicates fluid The kinematic viscosity coefficient, Indicates fluid The kinematic viscosity coefficient.
[0101] Calculating chemical potential and external force The specific formula is as follows:
[0102]
[0103] in, Indicates the order parameter The average value, .parameter and It is related to the interface thickness and surface tension The relevant physical constants are specifically defined as follows: .
[0104] The specific formulas for calculating the flow field velocity and hydrodynamic pressure within the solution region are as follows:
[0105]
[0106] in, Indicates the flow field velocity. Represents discrete velocity. This represents the distribution function of the flow field after migration. Since the distribution function after migration at the current moment is the same as the distribution function used for collision at the next moment, the time will not be specified further. Indicates the time step. Indicates the speed of sound. This represents the weight at position 0, which is 4 / 9 in the D2Q model. The gradient representing density. As an auxiliary variable, .
[0107] In this embodiment, the equilibrium distribution function of the migrated phase field for:
[0108] in, Indicates the weighting coefficient. Indicates free parameters, Indicates the order parameter, Indicates the speed of sound. Represents discrete velocity. Indicates the flow field velocity. Indicators, ; Equilibrium distribution function of the post-migration phase flow field for:
[0109] in, Indicates dynamic water pressure. Indicates fluid density, As an auxiliary variable, , Represents discrete velocity; Source term distribution function of the migrated phase field for:
[0110] in, For the order parameter The function has , Indicates fluid The order parameter, Indicates fluid The order parameter, Indicates the interface thickness. Represents the unit normal vector. The partial derivative with respect to time t; Source term distribution function of the migrated flow field for:
[0111] in, Indicates the external force term. This indicates the density of the liquid phase. This represents the density of the gas phase. Indicates the order parameter The gradient.
[0112]
[0113]
[0114] and The coupling is the macroscopic equation of the two-phase flow model of this invention, wherein For mobility, It is the unit normal vector. It is about The function is defined as:
[0115] in, It is the fluid density. It is dynamic water pressure. It is dynamic viscosity, which can be obtained by... The calculation yielded, where It is the kinematic viscosity coefficient. External force terms, such as those including volume forces. and surface tension .
[0116] like Figure 2 and Figure 3 This is an explanatory diagram of the meshless algorithm of this invention. Figure 4 , Figure 5 , Figure 6 , Figure 7 , Figure 8 , Figure 9 and Figure 10 This is a schematic diagram illustrating a calculation example of two-phase flow within a regular region according to the present invention. Figure 11 , Figure 12 , Figure 13 and Figure 14 This is a schematic diagram illustrating a calculation example of two-phase flow in an irregular region according to the present invention.
[0117] Compared to the limitation of time step size by spatial step size in traditional LBM, this invention can independently adjust temporal and spatial resolution while maintaining stability, thereby improving overall simulation efficiency and adaptability. Furthermore, this invention further improves interface capture accuracy and flow stability by introducing a multiple relaxation time (MRT) model.
[0118] Example 2: Based on Example 1, this embodiment of the invention provides a simulation system for high density ratio two-phase flow, which can be used to implement the simulation method for high density ratio two-phase flow as described in the foregoing embodiments. The system includes: The first module is used to generate and initialize irregular nodes within the computational region; The second module is used to perform collision operations on the distribution functions of the flow field and phase field based on irregular nodes; The third module is used to trace the distribution functions of the flow field and phase field in reverse along their characteristic directions at irregular nodes based on the semi-Lagrange method to obtain the starting point; The fourth module is used to interpolate and reconstruct the distribution functions of the flow field and phase field at the starting point using the moving least squares approximation, based on the starting point, to obtain the distribution functions of the flow field and phase field after the MLS approximation. The fifth module performs a migration step based on the distribution functions of the flow field and the phase field after the MLS approximation, and updates the computation nodes to obtain the distribution functions of the flow field and the phase field after migration. The sixth module is used to calculate the physical quantities at the current time step using the distribution function of the flow field and the distribution function of the phase field after migration. The seventh module is used to calculate the equilibrium distribution function and source term distribution function of the flow field and phase field after migration, and repeatedly perform collision migration operation on the distribution functions of the flow field and phase field to calculate the physical quantities at each time step until the preset number of simulation time steps or convergence criteria are met, thus completing the simulation of high density ratio two-phase flow.
[0119] According to embodiments of the present invention, the present invention also provides an electronic device, a readable storage medium, and a computer program product.
[0120] In an exemplary embodiment, an electronic device includes: at least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to enable the at least one processor to perform a simulation method for a high density ratio two-phase flow as described in Embodiment 1 above.
[0121] In an exemplary embodiment, the readable storage medium may be a non-transient computer-readable storage medium storing computer instructions for causing a computer to execute the simulation method of a high-density-ratio two-phase flow according to Embodiment 1 above.
[0122] In an exemplary embodiment, the computer program product includes a computer program that, when executed by a processor, implements the simulation method for a high density ratio two-phase flow according to Embodiment 1 above.
[0123] The program code used to implement the methods of the present invention can be written in any combination of one or more programming languages. This program code can be provided to a processor or controller of a general-purpose computer, special-purpose computer, or other programmable data processing device, such that when executed by the processor or controller, the program code causes the functions / operations specified in the flowcharts and / or block diagrams to be implemented. The program code can be executed entirely on the machine, partially on the machine, as a standalone software package partially on the machine and partially on a remote machine, or entirely on a remote machine or server.
[0124] In the context of this invention, a machine-readable medium can be a tangible medium that may contain or store a program for use by or in conjunction with an instruction execution system, apparatus, or device. A machine-readable medium can be a machine-readable signal medium or a machine-readable storage medium. Machine-readable media can include, but are not limited to, electronic, magnetic, optical, electromagnetic, infrared, or semiconductor systems, apparatus, or devices, or any suitable combination of the foregoing. More specific examples of machine-readable storage media include electrical connections based on one or more wires, portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fibers, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination of the foregoing.
[0125] To provide interaction with a user, the systems and techniques described herein can be implemented on a computer having: a display device for displaying information to the user (e.g., a CRT (cathode ray tube) or LCD (liquid crystal display) monitor); and a keyboard and pointing device (e.g., a mouse or trackball) through which the user provides input to the computer. Other types of devices can also be used to provide interaction with the user; for example, feedback provided to the user can be any form of sensory feedback (e.g., visual feedback, auditory feedback, or tactile feedback); and input from the user can be received in any form (including sound input, voice input, or tactile input).
[0126] The systems and technologies described herein can be implemented in computing systems that include backend components (e.g., as a data server), or computing systems that include middleware components (e.g., an application server), or computing systems that include frontend components (e.g., a user computer with a graphical user interface or web browser through which a user can interact with embodiments of the systems and technologies described herein), or any combination of such backend, middleware, or frontend components. The components of the system can be interconnected via digital data communication of any form or medium (e.g., a communication network). Examples of communication networks include local area networks (LANs), wide area networks (WANs), and the Internet.
[0127] Computer systems can include clients and servers. Clients and servers are generally located far apart and typically interact via communication networks. Client-server relationships are created by computer programs running on the respective computers and having a client-server relationship with each other. Servers can be cloud servers, servers in distributed systems, or servers incorporating blockchain technology.
[0128] 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 this invention.
Claims
1. A simulation method for high density ratio two-phase flow, characterized in that, Includes the following steps: Generate irregular nodes within the computational region and initialize them; Based on irregular nodes, collision operations are performed on the distribution functions of the flow field and phase field; At irregular nodes, the distribution functions of the flow field and phase field are traced in reverse along their characteristic directions using the semi-Lagrange method to obtain the starting point; Based on the starting point, the distribution functions of the flow field and phase field at the migration time step at the starting point are reconstructed by interpolation using the moving least squares approximation, so as to obtain the distribution functions of the flow field and phase field after MLS approximation. The migration step is performed based on the distribution functions of the flow field and the phase field after the MLS approximation, and the computation nodes are updated to obtain the distribution functions of the flow field and the phase field after migration. Using the distribution function of the flow field and the distribution function of the phase field after migration, the physical quantities at the current time step are calculated; The equilibrium distribution functions and source term distribution functions of the flow field and phase field after migration are calculated, and the collision migration operation is repeatedly performed on the distribution functions of the flow field and phase field to calculate the physical quantities at each time step until the preset number of simulation time steps or convergence criteria are met, thus completing the simulation of a two-phase flow with a high density ratio.
2. The simulation method for high density ratio two-phase flow according to claim 1, characterized in that, Collision migration operations on the distribution functions of the flow field and phase field are implemented using the MRT model; the dynamic behavior of the flow field is described by solving the Navier-Stokes equations using the lattice Boltzmann method, and the interface evolution of the phase field is described by solving the Allen-Cahn equations using the lattice Boltzmann method. The MRT evolution of the MRT-LB model of the Navier-Stokes equations is as follows: in, The distribution function representing the flow field. Indicates the time step. Indicates the position of irregular nodes. Indicates time, Indicates the starting point. Represents the Kronecker symbol, when hour ,when hour , Represents the invertible collision matrix The Middle Line number The elements of the column are , Denotes a diagonal relaxation matrix. Represents the transformation matrix. The source term distribution function represents the flow field. The equilibrium distribution function representing the flow field; The MRT evolution of the conservative ACE-LB model of the Allen-Cahn equation is as follows: in, Indicates the order parameter The distribution function, Represents the invertible collision matrix The Middle Line number The elements of the column are , Denotes a diagonal relaxation matrix. The source term distribution function represents the phase field. This represents the equilibrium distribution function of the phase field.
3. The simulation method for high density ratio two-phase flow according to claim 1, characterized in that, The specific formula for calculating the starting point is: in, Indicates the starting point. Indicates the position of irregular nodes. Represents discrete velocity. Indicates the time step.
4. The simulation method for high density ratio two-phase flow according to claim 1, characterized in that, Based on the starting point, the distribution functions of the flow field and phase field at the migration time step at the starting point are interpolated and reconstructed using the moving least squares approximation to obtain the distribution functions of the flow field and phase field after MLS approximation. The specific steps include: Define the generalized distribution function The distribution function representing the flow field or phase field, the generalized distribution function Defined in the region middle, Indicates that it includes all The space of dimensional vectors, the generalized distribution function The approximate function is : in, Denotes basis functions. Denotes the coefficients to be determined. Indicates the number of basis functions. Denotes the offset polynomial basis functions, superscript To represent transpose, we have: in, Indicates a reference point. and Represents irregular nodes exist direction and directional components, and Indicates reference point In direction and Component of direction; The support domain of the coefficients to be determined is selected, and a weight function corresponding to each meshless discrete node and its influence radius is constructed within the support domain using cubic splines; wherein, the support domain of the coefficients to be determined includes A gridless discrete node ; Based on the weighting function, the coefficients to be determined are determined by minimizing the weighted squared error of all meshless discrete nodes in the support domain: in, For Gram matrices, the matrix is... ,matrix , Indicates the corresponding to the first The weight function of a gridless discrete node and its influence radius, matrix superscript Indicates transpose. Indicates the first The generalized distribution function corresponding to a gridless discrete node; Substitute the coefficients to be determined into the approximate function. ,get: in, Indicates the number of coefficients to be determined; make The generalized distribution functions of the flow field and phase field after MLS approximation are obtained as follows: in, , This represents the flow field distribution function after MLS approximation. This represents the phase field distribution function after the MLS approximation. Representing shape functions The transpose of .
5. The simulation method for high density ratio two-phase flow according to claim 1, characterized in that, Based on the distribution functions of the flow field and phase field after the MLS approximation, a migration step is performed to obtain the distribution functions of the migrated flow field and phase field, specifically: At the starting point The flow field and phase field migration steps are performed above, and the distribution function of the migrated flow field is: in, The distribution function represents the flow field after migration. Indicates the starting point The flow field distribution function after MLS approximation; The distribution function of the phase field after migration is: in, The distribution function representing the phase field after migration. Indicates the starting point The phase field distribution function after MLS approximation.
6. The simulation method for high density ratio two-phase flow according to claim 1, characterized in that, The physical quantities include sequence parameters, fluid density, kinematic viscosity coefficient, chemical potential, external force terms, flow field velocity, and hydrodynamic pressure.
7. The simulation method for high density ratio two-phase flow according to claim 6, characterized in that, Order parameter The calculation formula is: in, The distribution function of the phase field after migration; fluid density The calculation formula is: in, Indicates fluid The order parameter, Indicates fluid The order parameter, Indicates fluid density, Indicates fluid density, kinematic viscosity coefficient The calculation formula is: in, Indicates fluid The kinematic viscosity coefficient, Indicates fluid The kinematic viscosity coefficient; Chemical potential The calculation formula is: in, Indicates the order parameter The average value, , Indicates the order parameter gradient, and Indicates the thickness of the interface and surface tension The relevant physical constants are specifically defined as follows: External force The calculation formula is: in, Indicates surface tension. Represents volume force; The formula for calculating the flow field velocity is: in, Indicates the flow field velocity. Represents discrete velocity. The distribution function represents the flow field after migration. Indicates the time step; Dynamic water pressure The calculation formula is: in, Indicates the speed of sound. This represents the weight at position 0. The gradient representing density. As an auxiliary variable, , This indicates the velocity at position 0.
8. The simulation method for high density ratio two-phase flow according to claim 1, characterized in that, Equilibrium distribution function of the phase field after migration for: in, Indicates the weighting coefficient. Indicates free parameters, Indicates the order parameter, Indicates the speed of sound. Represents discrete velocity. Indicates the flow field velocity. ; Equilibrium distribution function of the flow field after migration for: in, Indicates dynamic water pressure. Indicates fluid density, As an auxiliary variable, , Represents discrete velocity; Source term distribution function of the phase field after migration for: in, For the order parameter The function has , Indicates fluid The order parameter, Indicates fluid The order parameter, Indicates the interface thickness. Represents the unit normal vector. The partial derivative with respect to time t; Source term distribution function of the flow field after migration for: in, Indicates the external force term. This indicates the density of the liquid phase. This represents the density of the gas phase. Indicates the order parameter The gradient.
9. A simulation system based on a high density ratio two-phase flow, characterized in that, include: The first module is used to generate and initialize irregular nodes within the computational region; The second module is used to perform collision operations on the distribution functions of the flow field and phase field based on irregular nodes; The third module is used to trace the distribution functions of the flow field and phase field in reverse along their characteristic directions at irregular nodes based on the semi-Lagrange method to obtain the starting point; The fourth module is used to interpolate and reconstruct the distribution functions of the flow field and phase field at the starting point using the moving least squares approximation, based on the starting point, to obtain the distribution functions of the flow field and phase field after the MLS approximation. The fifth module performs a migration step based on the distribution functions of the flow field and phase field after MLS approximation, and updates the computation nodes to obtain the distribution functions of the flow field and phase field after migration. The sixth module is used to calculate the physical quantities at the current time step using the distribution function of the flow field and the distribution function of the phase field after migration. The seventh module is used to calculate the equilibrium distribution function and source term distribution function of the flow field and phase field after migration, and repeatedly perform collision migration operation on the distribution functions of the flow field and phase field to calculate the physical quantities at each time step until the preset number of simulation time steps or convergence criteria are met, thus completing the simulation of high density ratio two-phase flow.
10. A non-transitory computer-readable storage medium storing computer instructions, characterized in that, The computer instructions are used to cause the computer to perform the method according to any one of claims 1-8.