Parallel acceleration method applicable to two-phase color gradient lattice Boltzmann method

Abstract

The present invention discloses a parallel acceleration method applicable to a two-phase color gradient lattice Boltzmann method, which relates to the technical field of computational fluid dynamics. Taking the total number of fluid points and boundary solid points in a porous medium as the data space size, storage spaces are allocated for the fluid points and boundary solid points, and a storage information array is constructed; the storage information array includes the storage positions of each fluid point and each boundary solid point and the mapping relationship with the points in the adjacent directions; the porous medium is divided into regions to obtain multiple flow field regions; a graphics processor is allocated to each flow field region, and according to the storage information array, the parallel numerical simulation of the color gradient lattice Boltzmann method is respectively carried out on each flow field region through each graphics processor; during the collision simulation in the numerical simulation process, the matrix multiplication operation adopts a sparse matrix multiplication operation. This method can improve the efficiency of numerical simulation of multiphase flow by color gradient LBM.

Description

Technical Field

[0001] The present invention relates to the technical field of computational fluid dynamics, and particularly to a parallel acceleration method applicable to a two-phase color-gradient lattice Boltzmann method. Background Art

[0002] Currently, Computational Fluid Dynamics (CFD) is a science and engineering technology that uses numerical analysis and algorithms to solve fluid dynamics problems. CFD simulates fluid flow and its interaction with solid surfaces through computers and is widely applied in multiple fields such as aerospace, automotive engineering, chemical processes, weather forecasting, biomedicine, etc.

[0003] When dealing with multiphase flow problems, due to the advantages of the color-gradient lattice Boltzmann method (LBM) such as high numerical accuracy, strict mass conservation, small false velocity at the phase interface, good numerical stability under large viscosity ratios, and independent adjustment of surface tension parameters, the color-gradient LBM is widely used in multiphase flow simulations; for example, in fields such as chemical engineering, oil extraction, and environmental engineering, the color-gradient LBM is used to accurately simulate multiphase flows (such as oil-water separation, bubble dynamics, etc.).

[0004] However, the efficiency of numerical simulation of multiphase flow by the color-gradient LBM in the prior art is low. Summary of the Invention

[0005] Based on this, it is necessary to provide a parallel acceleration method applicable to a two-phase color-gradient lattice Boltzmann method for the above technical problems, which can improve the efficiency of numerical simulation of multiphase flow by the color-gradient LBM.

[0006] The present invention adopts the following technical solutions:

[0007] The present invention provides a parallel acceleration method applicable to a two-phase color-gradient lattice Boltzmann method, including:

[0008] Count the total number of fluid points and boundary solid points in the porous medium;

[0009] Using the total number as the data space size, allocate storage space for fluid points and boundary solid points, and construct a storage information array; the storage information array includes the storage positions of each fluid point and each boundary solid point and the mapping relationship with the points in the adjacent directions.

[0010] Divide the porous medium into regions to obtain multiple flow field regions;

[0011] A graphics processor is allocated to each flow field region, and based on the storage information array, parallel numerical simulations of the color-gradient lattice Boltzmann method for each flow field region are respectively performed by each graphics processor; during the numerical simulation, the matrix multiplication operation during the collision simulation adopts a sparse matrix multiplication operation.

[0012] Optionally, a storage information array is constructed, including:

[0013] In the original array [x, y, z] space of the porous medium, in the order of x, y, and z, serial numbers are sequentially assigned to the fluid points and boundary solid points in the storage information array;

[0014] A mapping relationship is established between each point in the storage space and the points in the adjacent directions, and the mapping relationship is saved in the storage information array as a 32-bit integer array.

[0015] Optionally, the porous medium is divided into multiple flow field regions, including:

[0016] Obtain the number of graphics processors for performing numerical simulations;

[0017] Based on the computing resources, the part of the porous medium that needs to allocate storage space is divided into multiple flow field regions; the number of flow field regions is the same as the number of graphics processors, and the difference between the total numbers of fluid points and boundary solid points in each flow field region is less than a preset threshold.

[0018] Optionally, parallel numerical simulations of the color-gradient lattice Boltzmann method for each flow field region of the porous medium include:

[0019] For any one flow field region, calculate the wall normal vector in the flow field region and perform parameter initialization on the flow field region; the initialization parameters include the initial velocity, density, distribution function of each phase, and the color function of the flow field region;

[0020] According to the initialization parameters, the color-gradient lattice Boltzmann method is iteratively executed for each flow field region until a preset termination condition is reached, and the distribution information of each phase in the porous medium is obtained.

[0021] Optionally, the color-gradient lattice Boltzmann method includes:

[0022] According to the color function, the unit interface normal vector at the solid boundary is calculated using the wetting boundary condition, and the color gradient and local interface curvature are calculated based on the unit interface normal vector;

[0023] According to the local interface curvature and velocity, the force term is calculated, and the collision process is performed on each lattice point based on the force term, color gradient, and distribution function to obtain the distribution function of each lattice point after collision;

[0024] Re-color each lattice point according to the distribution function after collision of each lattice point to obtain the two-phase distribution function after re-coloring;

[0025] Migrate each lattice point according to the two-phase distribution function after re-coloring to obtain the distribution function after migration;

[0026] Update the density and velocity of each phase according to the distribution function after migration, and update the color function of the flow field region according to the density of each phase;

[0027] Execute the boundary conditions and correct the distribution function of the boundary lattice points.

[0028] Optionally, calculate the unit interface normal vector at the solid boundary using the wetting boundary condition according to the color function, including:

[0029] Calculate the color gradient of the fluid lattice points in the flow field region according to the color function;

[0030] Determine the estimated unit normal vector perpendicular to the interface according to the color gradient;

[0031] Determine the first unit vector and the second unit vector according to the contact angle, the unit normal vector perpendicular to the wall, and the estimated unit normal vector perpendicular to the interface;

[0032] Obtain the Euclidean distances between the first unit vector and the second unit vector and the estimated unit normal vector perpendicular to the interface respectively;

[0033] Determine the unit vector corresponding to the minimum Euclidean distance as the unit interface normal vector at the solid boundary.

[0034] Optionally, the calculation formula for the local interface curvature is:

[0035] ;

[0036] Where, represents the local interface curvature, , and represent the components of the unit interface normal vector in the x , y , z directions respectively.

[0037] Optionally, the calculation formula for the force term is:

[0038] ;

[0039] ;

[0040] ;

[0041] Among them, represents the force term, represents the velocity, represents the multi-relaxation transformation matrix, represents the inverse of, represents the diagonal relaxation matrix, represents along the weight coefficient in the direction of; represents along the discrete velocity vector in the direction of, represents the surface tension, represents the repulsive force between foams when realizing foam flow, represents the lattice sound speed, represents the lattice time step, represents the surface tension coefficient, represents the local interface curvature, represents the color gradient.

[0042] Optionally, the collision process is:

[0043] ;

[0044] Among them, represents the position vector of the fluid point; represents at the moment after collision in space at the total distribution function in the direction of the discrete velocity ; represents at the moment before collision in space at the total distribution function in the direction of the discrete velocity ; represents the collision operator in the direction of; represents the external force term; the matrix multiplication operation in the collision operator adopts the sparse matrix multiplication operation.

[0045] The present invention provides a parallel acceleration device applicable to the two-phase color gradient lattice Boltzmann method, including:

[0046] A storage module, used to count the total number of fluid points and boundary solid points in the porous medium, use the total number as the data space size, allocate storage space for the fluid points and boundary solid points, and construct a storage information array; the storage information array includes the storage positions of each fluid point and each boundary solid point and the mapping relationship with the points in the adjacent directions;

[0047] A slicing module, used to divide the porous medium into regions to obtain multiple flow field regions;

[0048] A processing module is configured to allocate a graphics processor to each flow field region, and respectively perform parallel numerical simulation of the color gradient lattice Boltzmann method on each flow field region through each graphics processor according to the storage information array; during the numerical simulation, the matrix multiplication operation during the collision simulation adopts the sparse matrix multiplication operation.

[0049] The present invention provides a computer-readable storage medium storing a computer program, and when the computer program is executed by a processor, the above-mentioned parallel acceleration method applicable to the two-phase color gradient lattice Boltzmann method is realized.

[0050] The present invention provides a computer device including a memory, a processor, and a computer program stored on the memory and executable on the processor, and when the processor executes the program, the above-mentioned parallel acceleration method applicable to the two-phase color gradient lattice Boltzmann method is realized.

[0051] The above-mentioned at least one technical solution adopted by the present invention can achieve the following beneficial effects:

[0052] In the present invention, when performing parallel numerical simulation of the color gradient lattice Boltzmann method on each flow field region of the porous medium, the porous medium is decomposed into multiple flow field regions, and multiple graphics processors simultaneously perform numerical simulation on each flow field region, and the sparse matrix multiplication operation is adopted during the collision process simulation, which significantly improves the calculation speed; moreover, the porous medium is stored using a one-dimensional storage information array, that is, the way of only storing fluid information in the storage information array can avoid calculating the solid part, greatly reducing the memory required for calculation and improving the calculation speed, thereby greatly improving the efficiency of multi-phase flow numerical simulation of the porous medium through the color gradient LBM. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] The drawings described herein are used to provide a further understanding of the present invention and constitute a part of the present invention. The schematic embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute an improper limitation to the present invention. In the drawings:

[0054] Figure 1 is a schematic flow chart of a parallel acceleration method applicable to the two-phase color gradient lattice Boltzmann method provided by the present invention;

[0055] Figure 2 is a schematic comparison diagram of the memory allocation methods for the geometry of the porous medium by a traditional array and a sparse storage array provided by the present invention;

[0056] Figure 3 is a schematic diagram of the positions of lattice points in a traditional array and a sparse storage array provided by the present invention;

[0057] Figure 4Schematic diagram of the division of the flow field region provided by the present invention;

[0058] Figure 5 Schematic diagram of the process of the color gradient lattice Boltzmann method provided by the present invention;

[0059] Figure 6 Schematic diagram of the classification of lattice points taking a solid circle in the flow field as an example provided by the present invention;

[0060] Figure 7 Schematic diagram of the implementation of the wetting boundary condition provided by the present invention;

[0061] Figure 8 Schematic diagram of the discrete velocity directions of the D3Q19 discrete velocity model provided by the present invention;

[0062] Figure 9 Schematic diagram of the comparison of the simulation results of multi-foam flow in a circular pipe provided by the present invention;

[0063] Figure 10 Schematic diagram of the efficiency comparison between the parallel acceleration method applicable to two-phase color gradient LBM and the current mainstream OpenMP parallel CPU program provided by the present invention;

[0064] Figure 11 Parallel efficiency diagram of the parallel acceleration method applicable to two-phase color gradient LBM provided by the present invention;

[0065] Figure 12 Schematic diagram of the comparison between the simulation results of droplet generation in a microchannel by the present invention and previous experiments and numerical simulations;

[0066] Figure 13 Schematic diagram of the parallel acceleration device applicable to the two-phase color gradient lattice Boltzmann method provided by the present invention;

[0067] Figure 14 Schematic diagram of the computer device for implementing the parallel acceleration method applicable to the two-phase color gradient lattice Boltzmann method provided by the present invention. Detailed implementation manners

[0068] To make the objectives, technical solutions and advantages of the present invention clearer, the technical solutions of the present invention will be clearly and completely described below in conjunction with the specific embodiments of the present invention and the corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0069] In the prior art, there are many challenges in the algorithm implementation of color-gradient multiphase LBM. In terms of computational complexity, the simulation of two-phase flow usually involves complex physical phenomena such as interface capturing and interaction, resulting in high computational costs. Especially in high-resolution simulations, the computational volume increases exponentially. In terms of the efficiency of parallel computing, although LBM itself can be well parallelized, when dealing with large-scale grids and multiphase flows, data communication and load balancing become bottlenecks. How to optimize the parallel algorithm to improve computational efficiency is an important challenge. In terms of memory management, large-scale parallel computing needs to process a large amount of data. Especially in high-dimensional and high-precision simulations, the use and management of memory become particularly critical. The efficient access and caching mechanism of data need to be further optimized. In terms of algorithm stability and accuracy, in parallel computing, the stability and accuracy of the algorithm may be affected, especially when dealing with dynamic interfaces. Effective numerical methods need to be designed to ensure high efficiency and accuracy in a parallel environment. Under the influence of these factors, the computational scale of color-gradient multiphase LBM is often greatly limited, making it difficult to achieve numerical simulations of large-scale cases.

[0070] To overcome the above problems, the present invention proposes a parallel acceleration method applicable to the two-phase color-gradient lattice Boltzmann method, which has the following characteristics: First, it can be applied to two-dimensional and three-dimensional flows simultaneously and is applicable to complex two-phase flow simulations including flows in complex geometries and foam flows. Second, the program design of the two-phase color-gradient LBM is completed under the JAX framework, realizing large-scale automatic parallelization under multiple graphics processors. Third, for complex porous media geometries, one-dimensional sparse matrix storage is used, combined with optimization methods such as matrix multiplication optimization, greatly reducing the required video memory for calculation and improving the calculation speed.

[0071] The following will, in conjunction with the accompanying drawings, detail the technical solutions provided by each embodiment of the present invention.

[0072] Figure 1 The following is a schematic flowchart of a parallel acceleration method applicable to the two-phase color-gradient lattice Boltzmann method in the present invention, which specifically includes the following steps:

[0073] S101, count the total number of fluid points and boundary solid points in the porous medium, use the total number as the data space size, allocate storage space for the fluid points and boundary solid points, and construct a storage information array; the storage information array includes the storage positions of each fluid point and each boundary solid point and the mapping relationship with the points in the adjacent directions.

[0074] The size and density distribution of lattice points can be determined according to the requirements for simulation accuracy and the limitations of computing resources. Based on the size and density distribution of the lattice points, the porous medium is divided into lattice points, obtaining a plurality of lattice points. The lattice points include fluid points and boundary solid points, and the boundary solid points are solid points adjacent to at least one fluid point.

[0075] When simulating the flow in a porous medium, the ineffective occupation of memory and computing resources is an urgent problem to be solved. Since most of the space in the porous medium is occupied by solids, these solid points do not contain any flow field information and will not participate in the calculation during the evolution steps of the flow field. However, in the conventional array storage method, such as storing an array in the form of [x, y, z] in a computational domain of size, all the flow field physical information is still attached to the solid points. At this time, although the values of the flow field information on the solid points are all 0, the unified 64-bit precision of the array will cause these flow field information with values of 0 to occupy the same memory as the normal fluid points and introduce the same computational amount in the evolution of the flow field.

[0076] A feasible optimization scheme is to adopt a sparse storage scheme for the porous medium geometry. For example, rewrite the original [x, y, z] array as nx , where nx is the total number of fluid points and boundary solid points remaining after excluding all solid points from the original array. In this way, the storage and calculation are only carried out for the fluid points, which can save storage space and computing resources. As Figure 2 shown, Figure 2 Figure is a comparison of the memory allocation methods for the porous medium geometry between the traditional array and the sparse storage array. In the figure, the left figure is an example of a porous medium geometry, and the right figure is a schematic diagram of the memory allocation of the traditional array and the sparse storage array. Among them, the fluid part represents the part stored in the memory and participating in the calculation, and the solid part represents the part not stored in the memory and not participating in the calculation. However, this operation will lose the original geometric position relationship, and the program calculation of the present invention needs to use these geometric relationships. For example, in the migration step, it involves the information of other points around the fluid point. In the original array storage scheme, the points around a certain point can be expressed by, but in the sparse storage, since the original array is completely expanded and the geometric position relationship is missing, it is impossible to directly find the positions of the points around a certain point in the array. In addition, this method of completely excluding solid points will also make it impossible to implement the wetting boundary condition, resulting in its inapplicability to porous medium simulation.

[0077] In view of the above problems, the present invention innovatively proposes a special sparse storage method, including: counting the total number of fluid points and boundary solid points in the porous medium, using the total number as the data space size, allocating storage space for the fluid points and boundary solid points, and constructing a storage information array. The specific implementation steps for constructing the storage information array include: in the original array [x, y, z] space of the porous medium, in the order of x, y, and z, sequentially allocate the serial numbers of the fluid points and boundary solid points in the storage information array; establish a mapping relationship between each point in the storage space and the points in the adjacent directions, and save the mapping relationship in the storage information array as a 32-bit integer array.

[0078] Specifically, it includes the following steps:

[0079] S1. In the original [x, y, z] array space, count the total number nx of fluid points and boundary solid points.

[0080] S2. Use [nx] as the data space size to allocate storage space for all arrays.

[0081] S3. In the original [x, y, z] array space, in the order of , sequentially allocate the serial numbers of the fluid points and boundary solid points in [nx], as shown in Figure 3 . Figure 3 As shown in the figure, it is a schematic diagram of the positions of lattice points in a traditional array and a sparse storage array in a simple geometry. In the figure, the black part is the area that does not need to participate in the calculation, and the white part is the area that needs to participate in the calculation. The left figure shows the spatial position of the lattice points in the traditional array, and the right figure shows the spatial position of the lattice points in the one-dimensional sparse array.

[0082] S4. Taking D3Q19 as an example, establish a corresponding mapping relationship between the lattice points in 18 directions around each point in nx . For example, the point above in the original array space corresponds to , while in the storage method of the present invention nx , the point above corresponds to . Subsequently, save the mapping relationship as a 32-bit integer array with a size of [18, nx . Thereafter, the specific positions corresponding to the points in each direction around a certain point in nx can be directly obtained.

[0083] Through the above operations, the mapping reduction of the geometric positions in the sparse matrix can be achieved. During this process, the size of the 32-bit integer array established is negligible compared to the original 64-bit floating-point array, significantly reducing the memory occupancy. For example, for a case with a porosity of 10%, compared with the conventional array storage method, this sparse storage format can save more than 85% of the storage and computing resources.

[0084] S102. Divide the porous medium into regions to obtain multiple flow field regions.

[0085] The present invention adopts two different allocation methods to divide the porous medium, namely automatic parallelism and manual parallelism. Under automatic parallelism, the porous medium is divided into regions to obtain multiple flow field regions, including: obtaining the number of graphics processors used for numerical simulation; based on the computing resources, dividing the part of the porous medium that needs to allocate storage space into multiple flow field regions; the number of flow field regions is the same as the number of graphics processors, and the difference between the total numbers of fluid points and boundary solid points in each flow field region is less than a preset threshold. That is, in order to balance the calculations of all graphics processors and achieve load balancing, the total numbers of all fluid points and boundary solid points in each flow field region should be as identical as possible.

[0086] As Figure 4 shown, Figure 4 is a schematic diagram of the division of a flow field region. Specifically, the compiler selects the best computing strategy, and then the regions can be divided according to the computing strategy and computing resources, and the calculation amounts of all flow field regions are the same; it is only necessary to save the large array across multiple devices, and the compiler will partition all the internal contents and compile the communication between devices; while under manual parallelism, precise control over data saving and calculation domain allocation for different devices can be achieved, and there is a higher performance ceiling in some cases.

[0087] S103. Allocate a graphics processor to each flow field region, and respectively perform parallel numerical simulation of the color-gradient lattice Boltzmann method on each flow field region according to the storage information array; during the collision simulation in the numerical simulation process, the multiplication operation adopts the sparse matrix multiplication operation.

[0088] Allocate a graphics processing unit (GPU) to each flow field region to achieve parallel simulation through multiple GPUs; the present invention has completed the program design of this two-phase color-gradient LBM under the JAX framework, which provides a unified NumPy-style interface for running calculations on multiple GPUs in a local or distributed setting, provides built-in just-in-time (JIT) compilation capabilities, and can achieve automatic vectorization of functions.

[0089] The present invention uses code of Single-Program Multiple-Data (SPMD) (which allows different input data to implement the same calculation steps on multiple parallel devices, for example, simultaneously calculating the evolution of different parts in a flow field on multiple devices).

[0090] Optionally, perform parallel numerical simulation of the color gradient lattice Boltzmann method on each flow field region of the porous medium, including: for any flow field region, calculate the wall normal vector in the flow field region, and perform parameter initialization on the flow field region; the initialization parameters include the initial velocity, density, distribution function of each phase, and the color function of the flow field region; iteratively execute the color gradient lattice Boltzmann method on each flow field region according to the initialization parameters until a preset termination condition is reached, to obtain the distribution information of each phase in the porous medium. Among them, the termination condition may be that the number of iterations reaches a preset number threshold.

[0091] As Figure 5 shown, the specific steps of the color gradient lattice Boltzmann method include:

[0092] S501, according to the color function, calculate the unit interface normal vector at the solid boundary using the wetting boundary condition, and calculate the color gradient and local interface curvature according to the unit interface normal vector.

[0093] The implementation of the wetting boundary condition is a key issue in simulating complex geometric two-phase flows. The present invention has developed a wetting boundary condition in three dimensions in the Boltzmann method color model. The basic idea of this method is to adjust the direction of the color gradient at the three-phase contact line so that this direction satisfies the geometric condition of the specified contact angle .

[0094] To better explain the algorithm of the wetting boundary, here the lattice points are divided into two major categories: : Fluid lattice points; : Solid lattice points. At the same time, and can be further divided into two subcategories respectively. For the lattice points belonging to , here they are divided into two subcategories:

[0095] : Boundary fluid points, adjacent to at least one solid lattice point .

[0096] : Interior fluid points, not adjacent to any solid lattice point .

[0097] For those belonging to The grid points are also divided into two sub - categories:

[0098] : Boundary solid points, which are adjacent to at least one fluid grid point proximate.

[0099] : Interior solid points, which are not adjacent to any fluid grid point proximate.

[0100] Figure 6 is a schematic diagram of grid point classification with a solid circle in the flow field as an example. Among them, the black curve represents the solid boundary; the circles on the grid points belong to the grid points of ; the solid dots belong to ; the hollow dots belong to the grid points of ; the solid squares belong to ), this example clearly shows the specific situation of grid point classification near a solid circle. The three - dimensional case is similar to the two - dimensional case and will not be elaborated here.

[0101] Therefore, according to the color function, the unit interface normal vector at the solid boundary is calculated using the wetting boundary condition, including: calculating the estimated color gradient of the fluid grid points in the flow field region according to the color function; determining the estimated unit normal vector perpendicular to the interface according to the estimated color gradient; determining the first unit vector and the second unit vector according to the contact angle, the unit normal vector perpendicular to the wall, and the estimated unit normal vector perpendicular to the interface; obtaining the Euclidean distances between the first unit vector and the second unit vector and the estimated unit normal vector perpendicular to the interface respectively; and determining the unit vector corresponding to the minimum Euclidean distance as the unit interface normal vector at the solid boundary.

[0102] Among them, the fluid grid points include boundary fluid points and interior fluid points; calculating the estimated color gradient of the fluid grid points in the flow field region according to the color function, the specific process includes:

[0103] When calculating physical quantities such as the interface normal vector and the local interface curvature, the problem of finding partial derivatives will be involved. In order to minimize the discrete error generated when calculating partial derivatives, the present invention uses a fourth - order finite - difference algorithm to obtain the partial derivatives of variables. Taking the variable as an example, its partial derivative can be expressed as:

[0104] (1);

[0105] Among them, represents the lattice sound speed, represents along Discrete velocity vectors in the direction indicating the lattice time step indicating the position vector of the fluid point indicating the unit vector in the partial derivative direction indicating the direction corresponding to the partial derivative

[0106] Therefore, the present invention can calculate all the fluid lattice points using formula (1) The color gradient of the lattice points . However, when calculating the at the internal fluid point the color function at the boundary solid point is unknown. Therefore, the present invention obtains the estimated value of the color function at the boundary solid point by weighted averaging of the adjacent values of the color function , as shown in formula (2).

[0107] (2);

[0108] wherein represents the color function of lattice point x represents the weight coefficient along the direction

[0109] After obtaining the estimated value of at the boundary solid point from formula (2), the color gradient at all fluid lattice points can be directly calculated by formula (1). However, since the at the boundary fluid point is calculated from the estimated value of at the boundary solid point , it does not satisfy the specified contact angle condition. Therefore, the present invention refers to this color gradient as the estimated color gradient, denoted as , and the estimated color gradient needs to be corrected to meet the specified contact angle condition

[0110] To satisfy the specified contact angle , the simplest and most effective method to correct the estimated color gradient is to keep the modulus of unchanged and only correct the direction of . At the direction of can be expressed as the estimated unit normal vector :

[0111] (3);

[0112] Among them, represents the modulus of.

[0113] For the grid points belonging to the boundary fluid points two types of unit normal vectors are defined: one is the unit normal vector perpendicular to the wall , and the other is the unit normal vector perpendicular to the fluid interface . The present invention uses an eighth-order isotropic discretization method to obtain the unit normal vector perpendicular to the wall at :

[0114] (4);

[0115] Among them, represents the discrete velocity in the direction, represents the modulus of, represents the weight coefficient of the eighth-order discretization method; represents the fluid-solid indicator function, is 0 when is 1 when. In the three-dimensional case, takes the value of:

[0116] (5).

[0117] Using an eighth-order isotropic discretization algorithm to replace formula (1) to calculate the normal vector perpendicular to the wall can well reduce the spurious velocity, especially reduce the spurious velocity in the case of a curved solid wall.

[0118] After obtaining the normal vector perpendicular to the wall, the unit normal vector perpendicular to the fluid interface is obtained by the following method: As Figure 7 shown, for each contact line therein, there are two unit vectors that satisfy the angle with being the contact angle (that is, the unit vectors and in the figure, respectively obtained by rotating the angle from counterclockwise and clockwise). These two unit vectors (the first unit vector and the second unit vector) can be calculated from and :

[0119] (6);

[0120] Among them, . After obtaining and , take the one with the closest Euclidean distance to as the corrected unit interface normal vector . For example, Figure 7 is a schematic diagram of the implementation of the wetting boundary condition, Figure 7 in which the upper part represents the blue fluid phase, the middle semi-circular part represents the red fluid phase, and the lower part represents the solid phase; is the unit normal vector perpendicular to the wall surface, is the unit normal vector perpendicular to the fluid interface, is the given contact angle, and are respectively obtained by rotating the angle from in the counterclockwise and clockwise directions, is the calculated correction angle. At the contact line on the left side of Figure 7 , take as the corrected unit interface normal vector , while at the contact line on the right side, take as the corrected unit interface normal vector .

[0121] Finally, correct the color gradient at according to the corrected unit interface normal vector to be .

[0122] Optionally, the local interface curvature can be calculated by the following formula:

[0123] (7);

[0124] Among them, represents the unit interface normal vector pointing into the blue fluid and can be written as , represents the curvature magnitude, represents the identity matrix, represents the Hamiltonian operator.

[0125] After some vector calculations, the local interface curvature can be written as:

[0126] (8);

[0127] Among them, , and respectively represent the components of the unit interface normal vector in the x , y , z directions.

[0128] S502. Calculate the force term according to the local interface curvature and velocity, and perform a collision process on each lattice point according to the force term, color gradient, and distribution function to obtain the distribution function of each lattice point after collision.

[0129] The present invention uses a lattice Boltzmann multi-relaxation color model to simulate the immiscible two-phase flow problem, considers a three-dimensional model and adopts a D3Q19 discrete velocity model. Since this three-dimensional model degenerates into a two-dimensional model when the number of lattice points in any dimension is 1, this model can simultaneously achieve applicability to two-dimensional and three-dimensional geometries.

[0130] In the model adopted by the present invention, the two-phase fluids are respectively called the red fluid and the blue fluid, and respectively represent the distribution functions of the red and blue fluids. The total distribution function of the two-phase fluids is defined as , and its collision process is:

[0131] (9);

[0132] wherein, represents the position vector of the fluid point; represents the total distribution function at the position in space after collision at time in the direction of the discrete velocity; represents the total distribution function at the position in space before collision at time in the direction of the discrete velocity; represents the collision operator in the direction; represents the external force term, and the multiplication operation in the force term uses a sparse matrix multiplication operation.

[0133] To improve the stability of the model when simulating two-phase fluids with high viscosity ratios and high elasticity, the present invention uses the Multiple Relaxation Time (MRT) collision model instead of the ordinary Bhatnagar-Gross-Krook (BGK) model. Compared with the collision term using the BGK approximation, the collision term of MRT can reduce non-physical spurious velocities and improve the numerical stability of the model when simulating high viscosity ratios, thus enhancing the ability of the two-phase flow model to simulate high viscosity ratio problems. In the MRT collision model, the collision operator is given by the following formula:

[0134] (10);

[0135] where, represents the multi-relaxation transformation matrix, represents the inverse of; represents the diagonal relaxation matrix; represents the equilibrium distribution function of.

[0136] The multi-relaxation transformation matrix is valued such that its moments with the distribution function contain as many actual physical quantities as possible, such as density, momentum, energy, etc. In the present invention, its value is:

[0137] (11).

[0138] The diagonal relaxation matrix is defined as:

[0139] (12);

[0140] where, is given as:

[0141] (13);

[0142] where, represents the relaxation time. The relaxation time determines the kinematic viscosity of the fluid:

[0143] (14);

[0144] where, represents the time interval.

[0145] is obtained by taking the Maxwell-Boltzmann distribution function with respect to the local velocity of the fluid The second-order Taylor expansion is carried out to obtain the expression as follows:

[0146] (15);

[0147] wherein, represents the total density of the fluid, , is the density of the red fluid, is the density of the blue fluid; represents the lattice sound speed, taking , represents the lattice length, is the lattice time step, where and both take 1; represents the weight coefficient along the direction; represents the discrete velocity vector along the direction.

[0148] In the present invention, for the D3Q19 discrete velocity model adopted in two-dimensional simulation, as shown in Figure 8 , Figure 8 is the schematic diagram of the discrete velocity direction of the D3Q19 discrete velocity model, and its weight coefficient and discrete velocity vector are respectively:

[0149] (16);

[0150] (17).

[0151] In the two-phase flow in porous media, gravity is very small relative to surface tension, and the flow is mainly dominated by viscous force and surface tension. At this time, the influence of the density ratio of the two-phase fluid can be ignored. Therefore, for simplicity, the present invention assumes that the densities of the red fluid and the blue fluid are equal. However, the viscosities of the two-phase fluids have a great influence on the flow in the porous media. In order to simulate the situation where the viscosity differences of the two-phase fluids in the gas-liquid two-phase flow are relatively large, the present invention adopts the harmonic average to determine the viscosity of the two-phase fluid mixture at the interface:

[0152] (18);

[0153] wherein, represents the viscosity of the two-phase fluid mixture, represents the color function; represents the kinematic viscosity of the red fluid; represents the kinematic viscosity of the blue fluid.

[0154] In this color model, the present invention adopts an external body force to achieve the surface tension and the influence on fluid flow. The expression is:

[0155] (19);

[0156] wherein, can be written as:

[0157] (20);

[0158] wherein, represents the weight coefficient along the direction; represents the discrete velocity vector along the direction, represents the local flow velocity, represents the lattice sound speed, taking , represents the surface tension, represents the repulsive force between bubbles when realizing bubble flow, and this term can be ignored when simulating non-bubble flow.

[0159] Using the concept of continuous interface force, the surface tension can be further expressed as:

[0160] (21);

[0161] wherein, represents the surface tension coefficient; represents the color function, and the color function is defined as:

[0162] (22);

[0163] wherein, represents the density of the red fluid; represents the density of the blue fluid. From the above formula, it can be obtained that: when , the area is the red fluid; when , the area is the blue fluid; when , it is the mixed area of the red and blue fluids, that is, the interface.

[0164] S503. According to the distribution function after the collision of each lattice point, re-color each lattice point to obtain the re-colored two-phase distribution function.

[0165] To ensure the separation of the two-phase fluid and generate a clear interface, a recoloring step is performed. The present invention adopts a recoloring algorithm, which can make the red fluid and the blue fluid moderately miscible at the interface and keep the distribution of the two-phase fluid consistent with the color function, so as to reduce the false velocity and overcome the problem of lattice locking. The recolored two-phase distribution function includes the distribution function of the recolored red fluid and the distribution function of the blue fluid. The distribution function of the recolored red fluid and the distribution function of the blue fluid are respectively expressed as:

[0166] (23);

[0167] (24);

[0168] wherein, represents the distribution function of the recolored red fluid, represents the distribution function of the recolored blue fluid, represents the distribution function after lattice collision, represents the total density of the fluid, , is the density of the red fluid, is the density of the blue fluid, represents the weight coefficient along direction, represents the discrete velocity vector along direction, represents the color gradient of the fluid, ( ) represents the phase interface separation parameter related to the interface thickness. To ensure numerical stability and model accuracy, the present invention sets to 0.7.

[0169] S504. According to the recolored two-phase distribution function, each lattice point is migrated to obtain the migrated distribution function.

[0170] After the recoloring step, the distribution function information of the red and blue fluids at each lattice point is migrated so that they are respectively migrated to adjacent lattice points. The expression is:

[0171] (25).

[0172] S505. According to the migrated distribution function, the density and velocity of each phase are updated, and according to the density of each phase, the color function of the flow field region is updated.

[0173] The migrated distribution function obtained by formula (25) is used to calculate the density of the two-phase fluid, that is:

[0174] (26).

[0175] To eliminate the discrete error caused by the acting force when restoring to the N - S equation, the local flow velocity is:

[0176] (27).

[0177] Update the color function of the flow field region according to formula (22).

[0178] S506. Execute the boundary conditions and correct the distribution function of the boundary lattice points.

[0179] Among them, the boundary conditions include inlet - outlet boundary conditions, wall - bounce boundary conditions, etc. According to the boundary conditions, the distribution function of the boundary lattice points is corrected.

[0180] For example, taking the wall - bounce boundary condition as an example, in the LBM simulation, when a fluid particle (the distribution function moving at a discrete velocity can be regarded as a representation of the fluid particle) hits the wall, it will bounce back. Step 1: Determine the wall position and the collision lattice points: First, the position of the wall in the calculation region needs to be determined, which is achieved through the definition of the geometric shape. Then, find out which lattice points are the boundary lattice points that collide with the wall. For example, in a two - dimensional LBM simulation, for a square calculation region, the position of the boundary wall is easy to determine, and the lattice points located on the wall are the collision lattice points. Step 2: Perform distribution function correction - bounce rule: For the correction of the distribution function of the boundary lattice points, according to the bounce rule, when a particle hits the wall, its velocity direction will reverse. For example, for a vertical wall, if a particle hits the wall with a velocity to the right, after bouncing back, its velocity direction becomes to the left, and the corresponding distribution function will be exchanged and corrected according to this rule.

[0181] In an exemplary embodiment, in the D3Q19 MRT collision process of the present invention, M, as shown in formula (10), involves a 19 - order complex matrix multiplication, which occupies a relatively large amount of computing resources, and it is very time - consuming to complete this process using the traditional matrix multiplication method. Therefore, the present invention adopts the following sparse matrix multiplication optimization method:

[0182] Since in is a diagonal matrix with a high degree of sparsity, directly calling the matrix multiplication of the algorithm library will calculate a large number of redundant terms, wasting computing resources. Therefore, the present invention directly expands this matrix multiplication, calculates each term separately, and stores the repeated calculation part in as local variables, further reducing the amount of calculation.

[0183] When expanding this matrix multiplication, each term requires 19 multiplication operations and 18 addition operations. The commonly used method is to perform cumulative assignment in a for loop, which is more concise in program writing. However, the cumulative operation in the for loop reads and writes to memory more frequently, which is not conducive to the calculation speed. Therefore, in the present invention, the calculation of each term is separately expanded, reducing the number of memory reads and writes from the original 57 times and 19 times to 38 times and 1 time respectively, significantly reducing the memory read and write volume.

[0184] After adopting these measures, the present invention can reduce the time-consuming of sparse matrix multiplication to about 20% of the original, greatly optimizing the efficiency of sparse matrix multiplication.

[0185] In an exemplary embodiment, the present invention adopts a single program multiple data (a parallel technology that allows different input data to implement the same calculation steps on multiple parallel devices. For example, the evolution of different parts of the flow field can be calculated simultaneously on multiple devices) code. In terms of the allocation of computing resources for different devices, the present invention can adopt two different allocation methods, namely automatic parallelism and manual parallelism. In automatic parallelism, the compiler selects the best calculation strategy, and the user only needs to slice the data according to the specific example situation and computing resources, and save the large array across multiple devices. The compiler will partition all the internal content and compile the communication between devices. In manual parallelism, precise control over the data storage and calculation domain allocation of different devices can be achieved, with a higher performance ceiling in some examples.

[0186] The specific implementation steps of this program are as follows:

[0187] (1) Import geometric information, classify grid points according to their attributes, and count the number of valid grid points.

[0188] (2) Slice the data of each physical quantity array in the flow field.

[0189] (3) Sparsify the array, obtain the mapping relationship array, and calculate the wall normal vector.

[0190] (4) Initialize the flow field with the required initial conditions.

[0191] (5) Calculate macroscopic physical quantities such as density, velocity, and color function, and execute the inlet and outlet boundary conditions.

[0192] (6) Use the wetting boundary condition to correct the interface normal vector at the solid boundary, and calculate the interface curvature.

[0193] (7) Calculate the force term according to the interface curvature and execute the collision step.

[0194] (8) Re-color.

[0195] (9) Execute the migration step and the wall bounce boundary condition.

[0196] (10) Repeat steps (5)-(9) above until the termination condition is reached.

[0197] The present invention proposes an LBM method that is applicable to both two-dimensional and three-dimensional cases and any arbitrary complex geometry under the color gradient lattice Boltzmann framework. This method is easy to apply, can utilize multi-GPU to accelerate the calculation of super-large grid cases, and has high accuracy, and can accurately achieve the numerical simulation of various complex flows.

[0198] In an exemplary embodiment, in order to verify the accuracy and stability of the parallel acceleration method applicable to two-phase color gradient LBM, the simulation of multi-foam flow in a circular pipe is carried out in this embodiment. This simulation involves the use of sparse storage, the implementation of wetting boundary conditions, the implementation of two-phase foam flow and repulsive force, and is very suitable for verifying the accuracy and stability of the parallel acceleration method applicable to two-phase color gradient LBM proposed by the present invention.

[0199] The specific description of this embodiment is as follows: A uniform-sized foam group with a volume fraction of is initially stationary in the circular pipe flow field and starts to flow along the direction of the circular pipe under the action of a certain pressure gradient, and finally reaches a steady state. Due to the repulsive force between the foams, different foams do not merge but maintain the foam group state, and the existence of the foam group will change the system flow rate at the steady state, which macroscopically manifests as a change in the equivalent viscosity of the two-phase flow system. Here, the ratio of the steady-state flow rate of the system without foam to the steady-state flow rate of the system with foam is called the equivalent viscosity.

[0200] In this embodiment, the factors to be concerned about are the flow pattern of the system and the value of the equivalent viscosity. As Figure 9 shown, Figure 9 is a schematic diagram comparing the simulation results of multi-foam flow in a circular pipe. The left and right figures respectively show the flow field pattern and the equivalent viscosity results when the flow reaches a steady state obtained by the previous single-GPU program and the multi-GPU program of the present invention under the same initial conditions; under the same initial conditions, compared with the previous single-GPU program, the results obtained by the present invention are exactly the same, which proves the stability and accuracy of the method of the present invention.

[0201] In an exemplary embodiment, to verify the efficiency of the parallel acceleration method for two-phase color gradient LBM provided by the present invention, in this embodiment, a comparison of the efficiency with the current mainstream OpenMP-parallelized Central Processing Unit (CPU) program was carried out. The test case was the problem of droplet generation in a microchannel. To ensure the fairness of the comparison, the test hardware was configured with high-end components of the same period: two AMD EPYC 7763 were used on the CPU side, and a single NVIDIA A100 40GB was used on the GPU side. The efficiency test results are as Figure 10 shown, Figure 10 Figure Figure 10 is a schematic diagram of the efficiency comparison between the parallel acceleration method for two-phase color gradient LBM provided by the present invention and the current mainstream OpenMP-parallelized CPU program. The test case was the problem of droplet generation in a microchannel with geometric dimensions of 2475×367×1. When comparing, two AMD EPYC 7763 were used on the CPU side, and a single NVIDIA A100 40GB was used on the GPU side. Figure 10 The meaning of the vertical axis speed in Figure 10 is the time (seconds) required for the program to iterate 10,000 steps. It can be seen that the calculation speed of the parallel acceleration method for two-phase color gradient LBM proposed by the present invention is much higher than that of the current mainstream OpenMP-parallelized CPU program.

[0202] In an exemplary embodiment, to verify the parallelism of the parallel acceleration method for two-phase color gradient LBM, that is, the efficiency when extended to more GPUs, the test case in Embodiment 1 was used in this embodiment, and the acceleration efficiency under different numbers of GPUs was verified by increasing the number of GPUs participating in the calculation. The results are as Figure 11 shown, Figure 11 Figure Figure 11 is the parallel efficiency graph of the parallel acceleration method for two-phase color gradient LBM. The dotted line in the figure is the ideal acceleration efficiency when the parallel efficiency is 100%, and the solid line is the actual acceleration efficiency of the present invention. Due to the existence of a large number of gradient solutions in two-phase flow, involving partial cross-device data transmission, there is a certain gap between the acceleration ratio and the ideal acceleration ratio. Even so, the acceleration ratio of the present invention is still as high as 3 when the total number of GPUs is 4, demonstrating good parallelism.

[0203] In an exemplary embodiment, this embodiment is used to illustrate the application scenario of the present invention. Droplet microfluidics is a technology for precisely manipulating the movement and processing of droplets at the micron scale, which has a wide range of applications in fields such as biomedicine, chemical synthesis, and materials science. Its advantages lie in high efficiency, low consumption, and multi-functional operation. Numerical simulation of droplet microfluidics is a key tool for understanding and optimizing droplet behavior, which can provide detailed information on complex flows and interfacial dynamics within the microfluidic system. However, the main difficulties currently faced by such numerical simulations include: the complexity of interfacial effects and multi-physics field coupling, the accurate capture of non-linear flow behavior, and the high computational cost of implementing numerical simulations under actual microchannel geometries. Since the present invention is applicable to the simulation of complex two-phase flows including flows in complex geometries and foam flows, and can efficiently calculate problems with large grid sizes on a multi-GPU platform, it has become an efficient numerical simulation tool for droplet microfluidics.

[0204] This embodiment demonstrates the numerical simulation of different arrangement patterns spontaneously formed by droplets during the injection process in a three-dimensional microchannel. This behavior of different arrangement patterns spontaneously formed by droplets during the injection process has been observed in earlier experiments, and previous numerical simulation methods were difficult to accurately simulate this behavior due to issues such as computational scale and algorithm accuracy. The comparison between the simulation results of droplet generation in the microchannel by the present invention and previous experiments and numerical simulations is as Figure 12 shown. It can be seen that the simulation accuracy and results of the present invention far exceed those of previous simulations, are basically completely consistent with the experiments, and at the same time provide complete flow field information that cannot be obtained in the experiments, which provides great convenience for further studying the flow mechanism.

[0205] It should be noted that this embodiment is only one of the feasible application scenarios of the present invention. In addition to microfluidic simulation, the present invention can also be applied to the following fields including but not limited to:

[0206] (1) Industrial applications: used to study multiphase flows in oil and gas production, gas-liquid two-phase flows in fuel cells, and droplet dynamics in inkjet printing.

[0207] (2) Biomedicine: Simulate the cell-fluid interaction of blood in biological fluids.

[0208] (3) Environmental engineering: Study the seepage behavior in porous media and evaluate the diffusion and migration of pollutants.

[0209] (4) Materials science: used to design new composite materials and optimize the microstructure of materials by simulating processes such as liquid-phase sintering and phase separation.

[0210] (5) Basic research: Explore the basic mechanisms of fluid interface phenomena such as surface tension and wetting effects, and provide theoretical support for experiments.

[0211] When applying the parallel acceleration method for the two-phase color gradient lattice Boltzmann method provided by the present invention, it is not necessary to execute according to the Figure 1 sequence of each step shown. The specific execution sequence of each step can be determined according to needs, and the present invention does not limit this.

[0212] The above is the parallel acceleration method for the two-phase color gradient lattice Boltzmann method provided by one or more embodiments of the present invention. Based on the same idea, the present invention also provides a corresponding parallel acceleration device for the two-phase color gradient lattice Boltzmann method, as Figure 12 shown.

[0213] Figure 13 FIG. is a schematic diagram of a parallel acceleration device for the two-phase color gradient lattice Boltzmann method provided by the present invention. The device 1300 includes:

[0214] A storage module 1301, configured to count the total number of fluid points and boundary solid points in the porous medium; use the total number as the data space size, allocate storage space for the fluid points and boundary solid points, and construct a storage information array; the storage information array includes the storage positions of each fluid point and each boundary solid point and the mapping relationship with the points in the adjacent directions;

[0215] A slicing module 1302, configured to divide the porous medium into regions to obtain multiple flow field regions;

[0216] A processing module 1303, configured to allocate a graphics processor to each flow field region, and respectively perform parallel numerical simulation of the color gradient lattice Boltzmann method on each flow field region according to the storage information array; in the process of collision simulation during the numerical simulation, the matrix multiplication operation adopts a sparse matrix multiplication operation.

[0217] For the specific limitations on the parallel acceleration device for the two-phase color gradient lattice Boltzmann method, reference can be made to the limitations on the parallel acceleration method for the two-phase color gradient lattice Boltzmann method in the above text, which will not be elaborated here. Each module in the above parallel acceleration device for the two-phase color gradient lattice Boltzmann method can be implemented in whole or in part by software, hardware, and their combination. The above modules can be embedded in the processor in the computer device in hardware form or independent of it, or stored in the memory in the computer device in software form, so that the processor can call and execute the operations corresponding to the above modules.

[0218] The present invention also provides a computer-readable storage medium, which stores a computer program, and the computer program can be used to execute the above Figure 1 provided parallel acceleration method for the two-phase color gradient lattice Boltzmann method.

[0219] The present invention also provides Figure 14 a schematic structural diagram of the computer device shown in the figure, such as Figure 14 shown. At the hardware level, the computer device includes a processor, an internal bus, a network interface, a memory, and a non-volatile memory. Of course, it may also include other hardware required for other services. The processor reads the corresponding computer program from the non-volatile memory into the memory and then runs it to implement the above Figure 1 parallel acceleration method applicable to the two-phase color gradient lattice Boltzmann method provided

[0220] Those of ordinary skill in the art can understand that all or part of the processes in the methods of the above embodiments can be completed by instructing relevant hardware through a computer program. The computer program can be stored in a non-volatile computer-readable storage medium. When the computer program is executed, it may include the processes of the embodiments of the above methods. Among them, any reference to a memory, storage, database, or other medium used in the embodiments provided by the present invention may include at least one of non-volatile and volatile memories. The non-volatile memory may include a read-only memory (ROM), magnetic tape, floppy disk, flash memory, or optical memory, etc. The volatile memory may include a random access memory (RAM) or an external cache memory. By way of illustration and not limitation, RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM), etc.

[0221] The technical features of the above embodiments can be combined arbitrarily. For the sake of brevity of description, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, it should be considered as the scope recorded in the present invention.

Claims

1. A parallel acceleration method for a two-phase color gradient lattice Boltzmann method, characterized in that: include: Count the total number of fluid points and boundary solid points in porous media; Using the total number as the data space size, allocating storage space for the fluid points and the boundary solid points, and constructing a storage information array; the storage information array includes the storage position of each fluid point and each boundary solid point and a mapping relationship with points in adjacent directions; Dividing the porous medium into regions to obtain multiple flow field regions; A graphics processor is allocated to each flow field area, and according to the storage information array, each graphics processor performs parallel numerical simulation of the color gradient lattice Boltzmann method on each flow field area; the matrix multiplication operation during the collision simulation in the numerical simulation process adopts sparse matrix multiplication operation; A parallel numerical simulation of the color gradient lattice Boltzmann method is performed on each of the flow field regions, comprising: For any flow field region, the wall normal vector in the flow field region is calculated, and the parameters of the flow field region are initialized; the initialization parameters include the initial velocity, density, distribution function of each phase and the color function of the flow field region; The color gradient lattice Boltzmann method is iteratively executed on each of the flow field regions according to the initialization parameters until a preset termination condition is reached, thereby obtaining the distribution information of each phase in the porous medium.

2. The method according to claim 1, characterized in that The construction of the storage information array includes: In the original array [x, y, z] space of the porous medium, sequentially assigning serial numbers in the storage information array to the fluid points and the boundary solid points in the order of x, y and z; A mapping relationship is established between each point in the storage space and a point in an adjacent direction, and the mapping relationship is stored in the storage information array as a 32-bit integer array.

3. The method according to claim 1, characterized in that The porous medium is divided into regions to obtain multiple flow field regions, including: Get the number of graphics processors used for numerical simulation; Based on computing resources, the portion of the porous medium to which storage space needs to be allocated is divided into a plurality of flow field regions; the number of flow field regions is the same as the number of graphics processors, and the difference between the total number of fluid points and boundary solid points in each flow field region is less than a preset threshold.

4. The method according to claim 1, characterized in that: The color gradient lattice Boltzmann method comprises: According to the color function, a unit interface normal vector at a solid boundary is calculated using a wetting boundary condition, and a color gradient and a local interface curvature are calculated according to the unit interface normal vector; Calculating a force term according to the local interface curvature and velocity, and performing a collision process on each grid point according to the force term, the color gradient and the distribution function to obtain a distribution function of each grid point after collision; Recoloring each of the grid points according to the distribution function after the collision of each grid point to obtain a recolored two-phase distribution function; Migrating each of the grid points according to the recolored two-phase distribution function to obtain a migrated distribution function; According to the distribution function after the migration, the density and velocity of each phase are updated, and according to the density of each phase, the color function of the flow field area is updated; Boundary conditions are enforced and distribution functions of the boundary grid points are modified.

5. The method according to claim 4, characterized in that The step of calculating the unit interface normal vector at the solid boundary using the wetting boundary condition according to the color function includes: Calculating an estimated color gradient of a fluid grid point in the flow field region according to the color function; Determining an estimated unit normal vector perpendicular to the interface according to the estimated color gradient; Determine a first unit vector and a second unit vector according to the contact angle, a unit normal vector perpendicular to the wall, and the estimated unit normal vector perpendicular to the interface; Obtaining the Euclidean distances between the first unit vector and the second unit vector and the estimated unit normal vector perpendicular to the interface; The unit vector corresponding to the minimum Euclidean distance is determined as the unit interface normal vector at the solid boundary.

6. The method according to claim 4, characterized in that The calculation formula of the local interface curvature is: Where κ represents the local interface curvature, n x 、n y and n z They represent the components of the unit interface normal vector n in the x, y, and z directions respectively.

7. The method according to claim 4, characterized in that The calculation formula of the force term is: in, represents the force term, u represents the velocity, M represents the multi-relaxation transformation matrix, M -1 represents the inverse of M, S represents the diagonal relaxation matrix, ω i represents the weight coefficient along the i direction; e i represents the discrete velocity vector along the i direction, F s represents the surface tension, F rep represents the repulsive force between bubbles when foam flow is realized, c s represents the lattice sound velocity, δ t represents the grid time step, σ represents the surface tension coefficient, κ represents the local interface curvature, Represents a color gradient.

8. The method according to claim 4, characterized in that The collision process is: Where x represents the position vector of the fluid point; represents the total distribution function along the discrete velocity i direction at space x after the collision at time t; f i (x, t) represents the total distribution function along the discrete velocity i direction at space x before the collision at time t; Ω i represents the collision operator in the i direction; Represents the external force term; the matrix multiplication operation in the collision operator adopts sparse matrix multiplication operation.

9. A parallel accelerator for a two-phase color gradient lattice Boltzmann method, characterized in that: include: A storage module is used to count the total number of fluid points and boundary solid points in the porous medium; Using the total number as the data space size, allocating storage space for the fluid points and the boundary solid points, and constructing a storage information array; the storage information array includes the storage position of each fluid point and each boundary solid point and a mapping relationship with points in adjacent directions; A segmentation module, used to divide the porous medium into regions to obtain multiple flow field regions; A processing module is used to allocate a graphics processor to each flow field area, and to perform parallel numerical simulation of the color gradient lattice Boltzmann method on each flow field area through each graphics processor according to the storage information array; the matrix multiplication operation during the collision simulation in the numerical simulation process adopts sparse matrix multiplication operation; A parallel numerical simulation of the color gradient lattice Boltzmann method is performed on each of the flow field areas, including: for any flow field area, calculating the wall normal vector in the flow field area, and initializing the parameters of the flow field area; the initialization parameters include the initial velocity, density, distribution function of each phase and the color function of the flow field area; according to the initialization parameters, the color gradient lattice Boltzmann method is iteratively executed on each of the flow field areas until a preset termination condition is reached, so as to obtain the distribution information of each phase in the porous medium.

Images