Parallel acceleration method suitable for two-phase color gradient lattice Boltzmann method

Abstract

The invention discloses a parallel acceleration method suitable for a two-phase color gradient lattice Boltzmann method, and relates to the technical field of computational fluid mechanics. By taking the total number of fluid points and boundary solid points in the porous medium as the size of a data space, distributing a storage space for the fluid points and the boundary solid points, and constructing a storage information array; the storage information array comprises storage positions of each fluid point and each boundary solid point and a mapping relation between the storage positions and points in the adjacent direction; carrying out regional division on the porous medium to obtain a plurality of flow field regions; a graphics processor is distributed to each flow field area, and according to the storage information array, parallel numerical simulation of a color gradient lattice Boltzmann method is carried out on each flow field area through each graphics processor; in the numerical simulation process, sparse matrix multiplication is adopted as matrix multiplication during collision simulation. The method can improve the efficiency of multiphase flow numerical simulation through the color gradient LBM.

Description

Technical Field

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

[0002] At present, 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 used in many fields, such as aerospace, automotive engineering, chemical processes, weather forecasting, biomedicine, etc.

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

[0004] However, the existing technology of numerical simulation of multiphase flow by color gradient LBM is less efficient. Summary of the invention

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

[0006] The present invention adopts the following technical solutions: The present invention provides a parallel acceleration method applicable to a two-phase color gradient lattice Boltzmann method, comprising: Count the total number of fluid points and boundary solid points in porous media; The total number is used as the data space size, storage space is allocated for fluid points and boundary solid points, and a storage information array is constructed; the storage information array includes the storage position of each fluid point and each boundary solid point and the mapping relationship with the points in the adjacent direction; Divide the porous medium into regions to obtain multiple flow field regions; A graphics processor is allocated to each flow field area. 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.

[0007] Optionally, construct a storage information array, including: In the original array [x, y, z] space of the porous medium, the fluid points and the boundary solid points are assigned serial numbers in the storage information array in the order of x, y and z; A mapping relationship is established between each point in the storage space and the points in the adjacent direction, and the mapping relationship is saved in the storage information array as a 32-bit integer array.

[0008] Optionally, 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 part of the porous medium that requires storage space is divided into multiple flow field areas; the number of flow field areas 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 area is less than a preset threshold.

[0009] Optionally, a parallel numerical simulation of the color gradient lattice Boltzmann method is performed on each flow field region of the porous medium, including: For any flow field area, the wall normal vector in the flow field area is calculated, and the parameters of the flow field area are initialized; the initialization parameters include the initial velocity, density, distribution function of each phase and the color function of the flow field area; The color gradient lattice Boltzmann method is iteratively executed on each flow field region according to the initialization parameters until the preset termination condition is reached, and the distribution information of each phase in the porous medium is obtained.

[0010] Optionally, a color gradient lattice Boltzmann method, comprising: 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 the local interface curvature are calculated according to the unit interface normal vector; According to the local interface curvature and velocity, the force term is calculated, and the collision process is performed on each grid point according to the force term, color gradient and distribution function to obtain the distribution function of each grid point after collision; According to the distribution function of each grid point after collision, each grid point is recolored to obtain the recolored two-phase distribution function; According to the recolored two-phase distribution function, each grid point is migrated to obtain the distribution function after migration; According to the distribution function after 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 the distribution function of the boundary grid points is modified.

[0011] Optionally, a unit interface normal vector at a solid boundary is calculated using a wetting boundary condition according to the color function, including: According to the color function, the color gradient of the fluid grid points in the flow field area is calculated; According to the color gradient, determine the estimated unit normal vector perpendicular to the interface; Determine a first unit vector quantity and a second unit vector according to the contact angle, a unit normal vector perpendicular to the wall, and an estimated unit normal vector perpendicular to the interface; Obtain 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.

[0012] Optionally, the calculation formula of the local interface curvature is: ; in, represents the local interface curvature, , and Denote the unit interface normal vector exist x , y , z Direction component.

[0013] Optionally, the force term is calculated as: ; ; ; in, represents the force term, Indicates speed, represents the multi-relaxation transformation matrix, express The inverse, represents the diagonal relaxation matrix, Indicates along Directional weight coefficient; Indicates along The discrete velocity vector in the direction, represents the surface tension, represents the repulsive force between bubbles when foam flow is realized, represents the lattice sound speed, represents the grid time step, represents the surface tension coefficient, represents the local interface curvature, Represents a color gradient.

[0014] Optionally, the collision process is: ; in, represents the position vector of the fluid point; express After the moment of collision in space The discrete velocity along Total distribution function of direction; express Before the moment of collision in space The discrete velocity along Total distribution function of direction; express Collision operator in direction; Represents the external force term; the matrix multiplication operation in the collision operator adopts sparse matrix multiplication operation.

[0015] The present invention provides a parallel acceleration device suitable for a two-phase color gradient lattice Boltzmann method, comprising: A storage module is 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 location of each fluid point and each boundary solid point and the mapping relationship with the points in the adjacent direction; A slicing module is used to divide the porous medium into regions to obtain multiple flow field regions; The processing module is used to allocate a graphics processor to each flow field area, and 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.

[0016] The present invention provides a computer-readable storage medium, wherein the storage medium stores a computer program, and when the computer program is executed by a processor, the parallel acceleration method applicable to the two-phase color gradient lattice Boltzmann method is implemented.

[0017] The present invention provides a computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the program, the parallel acceleration method applicable to the two-phase color gradient lattice Boltzmann method is implemented.

[0018] At least one of the above technical solutions adopted by the present invention can achieve the following beneficial effects: In the present invention, when performing parallel numerical simulation of the color gradient lattice Boltzmann method on each flow field area of ​​the porous medium, the porous medium is decomposed into multiple flow field areas, multiple graphics processors perform numerical simulation on each flow field area at the same time, and sparse matrix multiplication operation is used when simulating the collision process, which significantly improves the calculation speed; and, the porous medium is stored using a one-dimensional storage information array, that is, only the fluid information is stored in the storage information array, which can avoid calculation of the solid part, greatly reduce the amount of memory required for calculation, and improve the calculation speed, thereby greatly improving the efficiency of numerical simulation of multiphase flow in porous media through color gradient LBM. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] The drawings described herein are used to provide a further understanding of the present invention and constitute a part of the present invention. The exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute an improper limitation of the present invention. In the drawings: Figure 1 A schematic flow chart of a parallel acceleration method for a two-phase color gradient lattice Boltzmann method provided by the present invention; Figure 2 A schematic diagram comparing the memory allocation method of a traditional array and a sparse storage array for porous media geometry provided by the present invention; Figure 3 A schematic diagram of the position of a grid point in a traditional array and a sparse storage array provided by the present invention; Figure 4 A schematic diagram of dividing a flow field area provided by the present invention; Figure 5 A schematic diagram of a color gradient lattice Boltzmann method process provided by the present invention; Figure 6 A schematic diagram of grid point classification provided by the present invention, taking a solid circle in a flow field as an example; Figure 7 A schematic diagram of an implementation of a wetting boundary condition provided by the present invention; Figure 8 A schematic diagram of discrete velocity direction of a D3Q19 discrete velocity model provided by the present invention; Fig. 9 A schematic diagram for comparing simulation results of multi-foam flow in a circular tube provided by the present invention; Fig.10 A schematic diagram showing the efficiency comparison between a parallel acceleration method for two-phase color gradient LBM provided by the present invention and the current mainstream OpenMP parallel CPU program; Fig.11A parallel efficiency diagram of a parallel acceleration method for a two-phase color gradient LBM provided by the present invention; Fig.12 A schematic diagram showing the comparison between the simulation results of the present invention for droplet generation in a microfluidic channel and previous experiments and numerical simulations; Fig.13 A schematic diagram of a parallel acceleration device suitable for a two-phase color gradient lattice Boltzmann method provided by the present invention; Fig.14 A schematic diagram of a computer device for implementing a parallel acceleration method for a two-phase color gradient lattice Boltzmann method provided by the present invention. DETAILED DESCRIPTION

[0020] In order to make the purpose, technical solution and advantages of the present invention clearer, the technical solution 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 part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0021] 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 capture and interaction, resulting in high computational cost. Especially in high-resolution simulations, the amount of computation grows exponentially; in terms of the efficiency of parallel computing, although LBM itself can be well parallelized, data communication and load balancing become bottlenecks when dealing with large-scale grids and multiphase flows. How to optimize parallel algorithms to improve computational efficiency is an important challenge; in terms of memory management, large-scale parallel computing requires processing 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 needs 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. It is necessary to design effective numerical methods to ensure that they remain efficient and accurate 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 simulation of large-scale examples.

[0022] In order to overcome the above problems, the present invention proposes a parallel acceleration method suitable for a two-phase color gradient lattice Boltzmann method, which has the following characteristics: First, it can be applied to both two-dimensional and three-dimensional flows, and can be applied 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, one-dimensional sparse matrix storage is used for complex porous media geometry, combined with optimization methods such as matrix multiplication optimization, which greatly reduces the amount of video memory required for calculation and improves the calculation speed.

[0023] The technical solutions provided by various embodiments of the present invention are described in detail below in conjunction with the accompanying drawings.

[0024] Figure 1 The figure is a flow chart of a parallel acceleration method for a two-phase color gradient lattice Boltzmann method in the present invention, which specifically includes the following steps: S101, 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 storage information array includes the storage location of each fluid point and each boundary solid point and the mapping relationship with the points in the adjacent direction.

[0025] The size and density distribution of the grid points can be determined according to the requirements for simulation accuracy and the limitation of computing resources. According to the size and density distribution of the grid points, the porous medium is divided into grid points to obtain multiple grid points. The grid points include fluid points and boundary solid points. The boundary solid points are solid points adjacent to at least one fluid point.

[0026] When simulating flow in porous media, the ineffective use of memory and computing resources is an urgent problem to be solved. Since most of the space in porous media is occupied by solids, these solid points do not contain any flow field information and will not participate in the calculation in the flow field evolution step. However, conventional array storage methods, such as storing arrays in the form of [x, y, z] in a computational domain of size, will still attach all the flow field physical information to the solid points. At this time, although the flow field information on the solid points is 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 normal fluid points and introduce the same amount of calculation in the flow field evolution.

[0027] A feasible optimization solution is to use a sparse storage solution for porous media geometry, such as rewriting the original [x, y, z] array into [ nx ], here nx That is, the total number of fluid points and boundary solid points remaining after excluding all solid points in the original array. In this way, storage and calculation will only be performed on fluid points, which can save storage space and computing resources. Figure 2 As shown, Figure 2 This is a comparison of the memory allocation method for porous media geometry between traditional arrays and sparse storage arrays. In the figure, the left figure is an example of porous media geometry, and the right figure is a schematic diagram of the memory allocation method between traditional arrays and sparse storage arrays. The fluid part represents the part stored in the memory and participates in the calculation, and the solid part represents the part that is not stored in the memory and does not participate in the calculation. However, this operation will lose the original geometric position relationship, and the program calculation of the present invention requires these geometric relationships. For example, in the migration step, information about other points around the fluid point is involved. In the original array storage scheme, the points around a point can be expressed by, but in sparse storage, since the original array is fully expanded and the geometric position relationship is missing, it is impossible to directly find the position of the points around a point in the array. In addition, this method of completely excluding solid points will also make the wetting boundary conditions impossible to implement, resulting in its inability to be applied to porous media simulation.

[0028] 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 of 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 assigning serial numbers in the storage information array to the fluid points and boundary solid points; establishing a mapping relationship between each point in the storage space and the points in the adjacent direction, and saving the mapping relationship in the storage information array as a 32-bit integer array.

[0029] Specifically, the steps include: S1, count the total number nx of fluid points and boundary solid points in the original [x, y, z] array space.

[0030] S2, allocates storage space for all arrays with [nx] as the data space size.

[0031] S3, in the original [x, y, z] array space, The order of assigning the fluid points and boundary solid points their serial numbers in [nx] is as follows: Figure 3 As shown, Figure 3 This is a schematic diagram of the location of grid points in traditional arrays and sparse storage arrays 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 location of the grid points in the traditional array, while the right figure shows the spatial location of the grid points in the one-dimensional sparse array.

[0032] S4, taking D3Q19 as an example, nx] establishes corresponding mapping relationships between the grid points in 18 directions around each point in the original array space. The points above correspond to , and in the storage method of the present invention [ nx ]The point above corresponds to . Then, the mapping relationship is converted into a nx ] to save it, and then you can directly get the points in all directions around a certain point in [ nx ] in the corresponding specific location.

[0033] Through the above operations, the mapping and restoration of geometric positions in the sparse matrix can be achieved. The size of the 32-bit integer array established in this process is negligible compared to the original 64-bit floating-point array, which significantly reduces the memory usage. 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 storage and computing resources.

[0034] S102, dividing the porous medium into regions to obtain multiple flow field regions.

[0035] The present invention uses two different allocation methods to divide the porous medium into regions, namely, automatic parallel and manual parallel. Under automatic parallel, 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 computing resources, dividing the part of the porous medium to be allocated with 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 number 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 calculation and load balancing of all graphics processors, the total number of all fluid points and boundary solid points in each flow field region should be as equal as possible.

[0036] like Figure 4 As shown, Figure 4 The figure is a schematic diagram of the division of flow field areas. Specifically, the compiler selects the best computing strategy, and then divides the area according to the computing strategy and computing resources. The computing amount of all flow field areas is the same; large arrays can be saved across multiple devices, and the compiler will partition all internal contents and compile communication between devices; under manual parallelization, precise control of data storage and computing domain allocation of different devices can be achieved, and a higher performance ceiling can be achieved in some examples.

[0037] S103, 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 multiplication operation during the collision simulation in the numerical simulation process adopts sparse matrix multiplication operation.

[0038] A graphics processing unit (GPU) is assigned to each flow field area, and parallel simulation is achieved through multiple GPUs; the present invention completes the programming of the two-phase color gradient LBM under the JAX framework, which provides a unified NumPy-style interface for running calculations on multiple GPUs in local or distributed settings, provides a built-in just-in-time (JIT) compilation function, and can achieve automatic vectorization of functions.

[0039] The present invention uses a Single-Program Multiple-Data (SPMD) code (which allows different input data to implement the same computational steps on multiple parallel devices, such as computing the evolution of different parts of a flow field on multiple devices simultaneously).

[0040] Optionally, a parallel numerical simulation of the color gradient lattice Boltzmann method is performed on each flow field region of the porous medium, including: for any flow field region, calculating the wall normal vector in the flow field region, and initializing the parameters of 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; according to the initialization parameters, the color gradient lattice Boltzmann method is iteratively executed on each flow field region until a preset termination condition is reached, and the distribution information of each phase in the porous medium is obtained. The termination condition may be that the number of iterations reaches a preset threshold.

[0041] like Figure 5 As shown in FIG. 1 , the specific steps of the color gradient lattice Boltzmann method include: S501, according to the color function, a unit interface normal vector at the 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.

[0042] The implementation of wetting boundary conditions is a key issue in simulating complex geometric two-phase flows. The present invention develops wetting boundary conditions in three-dimensional cases in the Boltzmann method color model. The basic idea of ​​this method is to adjust the color gradient The direction at the three-phase contact line is such that the direction satisfies the specified contact angle geometric conditions.

[0043] To better explain the algorithm for wetting boundaries, the grid points are divided into two categories: : Fluid grid points; : Solid grid point. At the same time, and They can be divided into two categories: The grid points are divided into two categories: : Boundary fluid point, at least one solid grid point Approaching.

[0044] : Interior fluid point, not connected to any solid grid point Approaching.

[0045] For the The grid points are also divided into two categories: : Boundary solid point, at least one fluid grid point Approaching.

[0046] : Interior solid point, not connected to any fluid grid point Approaching.

[0047] Figure 6 The following is a schematic diagram of grid point classification, taking a solid circle in the flow field as an example. The black curve represents the solid boundary; the circle on the grid point belongs to The grid points of ; The hollow dots belong to ; The square on the grid point belongs to The grid points of ; the solid squares belong to ; The hollow square belongs to ), which clearly shows the specific situation of the classification of grid points near a solid circle. The three-dimensional case is similar to the two-dimensional case and will not be repeated here.

[0048] Therefore, according to the color function, the wetting boundary condition is used to calculate the unit interface normal vector at the solid boundary, including: calculating the estimated color gradient of the fluid grid points in the flow field area 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 quantity 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.

[0049] The fluid grid points include boundary fluid points and internal fluid points. According to the color function, the estimated color gradient of the fluid grid points in the flow field region is calculated. The specific process includes: When calculating physical quantities such as the interface normal vector and the local interface curvature, the problem of partial derivatives is 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 the variables. For example, its partial derivative can be expressed as: (1); in, represents the lattice sound speed, Indicates along The discrete velocity vector in the direction, represents the grid time step, represents the position vector of the fluid point, is the unit vector representing the direction of the partial derivative, Indicates the direction corresponding to the partial derivative.

[0050] Therefore, the present invention can use formula (1) to calculate all the fluid grid points Color gradient of grid points . But when calculating the internal fluid point Where The required boundary solid points The color function at The value is unknown, so the present invention uses the color function The weighted average of neighboring values ​​is used to obtain the boundary solid point Color function The estimated value of is shown in formula (2).

[0051] (2); in, represents the color function of the grid point x, Indicates along Direction weight coefficient.

[0052] The boundary solid point is obtained by formula (2): Department After the estimated value of , the equation (1) can be used to directly calculate the value of all fluid grid points. Color gradient at However, due to the boundary fluid point Department is a solid point through the boundary Department The estimated value is calculated, which cannot meet the specified contact angle condition. Therefore, the present invention calls the color gradient an estimated color gradient, which is expressed as , the estimated color gradient needs to be corrected to satisfy the specified contact angle condition.

[0053] In order to meet the specified contact angle , correct the estimated color gradient The simplest and most effective way is to keep The model remains unchanged, only the correction direction. Department The direction can be expressed as the estimated unit normal vector : (3); in, express Model.

[0054] For the boundary fluid points At the grid 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 adopts the eighth-order isotropic discrete method to obtain The unit normal vector perpendicular to the wall : (4); in, express Discrete speed in direction, express The model, express Model; Represents the weight coefficient of the eighth-order discrete method; represents the fluid-solid indicator function, The time value is 0, When the value is 1. In the three-dimensional case, The value of is: (5).

[0055] The eighth-order isotropic discrete algorithm is used to replace formula (1) to calculate the normal vector of the vertical wall. It can effectively reduce the false velocity, especially the false velocity in the case of curved surface and solid wall.

[0056] Get the normal vector of the vertical wall Then, the unit normal vector perpendicular to the fluid interface It is obtained by the following method: Figure 7 As shown, for each contact line, there are two unit vectors that satisfy The angle between the two is the contact angle (i.e., the unit vector in the figure and , respectively Rotate the angle counterclockwise and clockwise These two unit vectors (the first unit vector and the second unit vector) can be expressed by and The calculation results in: (6); in, . In obtaining and Then, take the middle distance between the two The closest Euclidean distance is taken as the corrected unit interface normal vector .For example, Figure 7 Schematic diagram of the implementation of wetting boundary conditions. Figure 7 In the figure, the upper half represents the blue fluid phase, the middle semicircle represents the red fluid phase, and the lower half represents the solid phase; is the unit normal vector perpendicular to the wall, is the unit normal vector perpendicular to the fluid interface, For a given contact angle, and Respectively by Rotate the angle counterclockwise and clockwise get, To calculate the correction angle, Figure 7 At the left contact line, take As the modified unit interface normal vector , and in the right contact line As the modified unit interface normal vector .

[0057] Finally, according to the corrected unit interface normal vector Color gradient at for .

[0058] Optionally, the local interface curvature It can be calculated by the following formula: (7); in, represents the unit interface normal vector pointing to the interior of the blue fluid and can be written as , represents the curvature size, represents the identity matrix, represents the Hamiltonian operator.

[0059] After some vector calculations, the local interface curvature It can be written as: (8); in, , and denote the unit interface normal vector exist x , y , z The direction component.

[0060] S502, calculating the 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 the distribution function of each grid point after the collision.

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

[0062] In the model adopted by the present invention, the two-phase fluids are respectively called red fluid and blue fluid. and They are represented as the distribution functions of the red and blue fluids, respectively. The total distribution function of the two-phase fluid is defined as , and its collision process is: (9); in, represents the position vector of the fluid point; express After the moment of collision in space The discrete velocity along Total distribution function of direction; express Before the moment of collision in space The discrete velocity along Total distribution function of direction; express Collision operator in direction; Represents the external force term, and the multiplication operation in the force term adopts sparse matrix multiplication operation.

[0063] In order to improve the stability of the model when simulating two-phase fluids with high viscosity ratio and high elasticity, the present invention uses a multiple relaxation time (MRT) collision model instead of the ordinary BGK (Bhatnagar-Gross-Krook, BGK) model. Compared with the collision term using the BGK approximation, the MRT collision term can reduce non-physical false velocity and improve the numerical stability of the model when simulating high viscosity ratios, thereby improving the ability of the two-phase flow model to simulate high viscosity ratio problems. In the MRT collision model, the collision operator It is given by the following formula: (10); in, represents the multi-relaxation transformation matrix, express The inverse of represents the diagonal relaxation matrix; express The equilibrium distribution function of .

[0064] Multi-relaxation transformation matrix The value of needs to make it and the moment of the distribution function contain as many actual physical quantities as possible, such as density, momentum, energy, etc. In the present invention, its value is: (11).

[0065] Diagonal relaxation matrix is defined as: (12); in, Given as: (13); in, Indicates relaxation time. The relaxation time determines the kinematic viscosity of the fluid. : (14); in, Indicates a time interval.

[0066] The Maxwell-Boltzmann distribution function is applied to the local velocity of the fluid. Perform a second-order Taylor expansion to obtain the expression: (15); in, 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 velocity, take , represents the grid length, is the grid time step, where and All take 1; Indicates along Directional weight coefficient; Indicates along The discrete velocity vector in the direction.

[0067] In the present invention, the D3Q19 discrete velocity model used in the two-dimensional simulation is as follows: Figure 8 As shown, Figure 8 The discrete velocity direction diagram of the D3Q19 discrete velocity model, and its weight coefficient and the discrete velocity vector They are: (16); (17).

[0068] In the two-phase flow of porous media, gravity is very small relative to surface tension, and the flow is mainly dominated by viscosity and surface tension. At this time, the density ratio of the two-phase fluid can be ignored. Therefore, for the sake of simplicity, the present invention assumes that the density of the red fluid and the blue fluid is equal. The viscosity of the two-phase fluid has a great influence on the flow in the porous medium. In order to simulate the situation where the viscosity of the two-phase fluid in the gas-liquid two-phase flow is greatly different, the present invention uses the harmonic mean to determine the viscosity of the two-phase fluid mixture at the interface: (18); in, 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.

[0069] In this color model, the present invention adopts an external volume force To realize the surface tension and its effect on fluid flow. The expression is: (19); in, It can be written as: (20); in, Indicates along Directional weight coefficient; Indicates along The discrete velocity vector in the direction, represents the local flow velocity, represents the lattice sound velocity, take , represents the surface tension, It represents the repulsive force between bubbles when realizing foam flow. This item can be ignored when simulating non-foam flow.

[0070] Using the concept of continuous interfacial forces, surface tension It can be further expressed as: (twenty one); in, represents the surface tension coefficient; Represents color function, color function It is defined as: (twenty two); in, represents the density of the red fluid; represents the density of the blue fluid. From the above formula, we can get: When , the area is red fluid; when When , the area is blue fluid; when , it is the mixing area of ​​the red and blue fluids, that is, the interface.

[0071] S503, recoloring each grid point according to the distribution function after the collision of each grid point, to obtain a recolored two-phase distribution function.

[0072] In order to ensure the separation of the two-phase fluid and produce 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 grid lock problem. 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 Respectively expressed as: (twenty three); (twenty four); in, represents the distribution function of the recolored red fluid, represents the distribution function of the recolored blue fluid, represents the distribution function after the grid point collision, represents the total density of the fluid, , is the density of the red fluid, is the density of the blue fluid, Indicates along The weight coefficient of the direction, Indicates along The discrete velocity vector in the 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 will Set to 0.7.

[0073] S504, according to the recolored two-phase distribution function, migrate each grid point to obtain a migrated distribution function.

[0074] After the recoloring step, the distribution function information of the red and blue fluids at each grid point is migrated to the adjacent grid points respectively. The expression is: (25).

[0075] S505, updating the density and velocity of each phase according to the distribution function after migration, and updating the color function of the flow field area according to the density of each phase.

[0076] The distribution function after migration obtained by formula (25) is used to calculate the density of the two-phase fluid, that is: (26).

[0077] In order to eliminate the discretization error caused by the force when restoring the NS equation, the local flow velocity for: (27).

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

[0079] S506, executing boundary conditions and modifying the distribution function of the boundary grid points.

[0080] Among them, the boundary conditions include inlet and outlet boundary conditions, wall rebound boundary conditions, etc. According to the boundary conditions, the distribution function of the boundary grid points is corrected.

[0081] For example, taking the wall rebound 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, a rebound occurs. Step 1: Determine the wall position and collision grid: First, determine the position of the wall in the calculation area, which is achieved by defining the geometric shape. Then, find out which grid points are the boundary grid points that collide with the wall. For example, in a two-dimensional LBM simulation, for a square calculation area, the position of the boundary wall is easy to determine, and the grid points located on the wall are the collision grid points. Step 2: Perform distribution function correction-bounce rule: For the distribution function correction of the boundary grid point, according to the rebound rule, when the particle hits the wall, its velocity direction will be reversed. For example, for a vertical wall, if the particle hits the wall at a speed to the right, its velocity direction will change to the left after rebounding, and the corresponding distribution function will also be exchanged and corrected according to this rule.

[0082] In an exemplary embodiment, the D3Q19 MRT collision process of the present invention, as shown in formula (10), involves The 19th-order complex matrix multiplication of occupies a lot of computing resources, and it is very time-consuming to complete the process using the traditional matrix multiplication method. To this end, the present invention adopts the following sparse matrix multiplication optimization method: because In is a diagonal matrix with high sparsity. Directly calling the matrix multiplication of the algorithm library will calculate a large number of redundant items, wasting computing resources. Therefore, the present invention directly expands the matrix multiplication, calculates each item separately, and converts The repeated calculation parts in are stored as local variables, which further reduces the amount of calculation.

[0083] When the matrix multiplication is expanded, each item needs to perform 19 multiplication calculations and 18 addition calculations. The commonly used method is to perform cumulative assignment in a for loop, which is more concise when programming, but the cumulative operation in the for loop reads and writes the memory more frequently, which is not conducive to the calculation speed. Therefore, the present invention expands the calculation of each item separately, so that the number of memory reads and writes is reduced from the original 57 times and 19 times to 38 times and 1 time, respectively, significantly reducing the memory read and write amount.

[0084] After taking these measures, the present invention can reduce the time consumption of sparse matrix multiplication to about 20% of the original time, greatly optimizing the efficiency of sparse matrix multiplication.

[0085] In an exemplary embodiment, the present invention uses 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, to simultaneously calculate the evolution of different parts of the flow field on multiple devices) code. In terms of the allocation of computing resources to different devices, the present invention can adopt two different allocation methods, namely automatic parallelization and manual parallelization. In automatic parallelization, the compiler selects the best computing strategy, and the user only needs to perform data slicing according to the specific case and computing resources, and save large arrays across multiple devices. The compiler will partition all internal contents and compile inter-device communication. In manual parallelization, precise control of data storage and computational domain allocation for different devices can be achieved, with a higher performance ceiling in some cases.

[0086] The specific implementation steps of the program are as follows: (1) Import geometric information and classify grid points according to their attributes, and count the number of valid grid points.

[0087] (2) Data slicing is performed on the arrays of various physical quantities in the flow field.

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

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

[0090] (5) Calculate macroscopic physical quantities such as density, velocity, and color functions and implement import and export boundary conditions.

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

[0092] (7) Based on the interface curvature, calculate the force term and execute the collision step.

[0093] (8) Recoloring.

[0094] (9) Execute the migration step and wall rebound boundary conditions.

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

[0096] In the framework of color gradient lattice Boltzmann, the present invention proposes an LBM method that is applicable to both two-dimensional and three-dimensional and arbitrary complex geometry. The method is simple to apply, can utilize multiple GPUs to accelerate the calculation of ultra-large grid examples with high accuracy, and can accurately realize numerical simulation of various complex flows.

[0097] In an exemplary embodiment, in order to verify the accuracy and stability of the parallel acceleration method for two-phase color gradient LBM, a simulation of multi-foam flow in a circular tube was performed in this embodiment. The simulation also involves the use of sparse storage, the implementation of wetting boundary conditions, the implementation of two-phase foam flow and repulsive force, which is very suitable for verifying the accuracy and stability of the parallel acceleration method for two-phase color gradient LBM proposed in the present invention.

[0098] The specific description of this embodiment is as follows: the volume fraction is The uniform size foam group initially remains still in the circular tube flow field, begins to flow along the circular tube under the action of a certain pressure gradient, and finally reaches a steady state. Due to the repulsive force between the foams, different foams will not merge but remain in the state of foam groups, and the existence of the foam group will change the system flow rate in the steady state, which is reflected in the change of the equivalent viscosity of the two-phase flow system on a macro scale. 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.

[0099] In this embodiment, the factors that need to be paid attention to are the flow pattern of the system and the value of the equivalent viscosity. Fig. 9 As shown, Fig. 9 It is a schematic diagram for comparing simulation results of multi-foam flow in a circular tube. The left and right figures are respectively the flow field morphology and 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, the results obtained by the present invention are completely consistent with those of the previous single GPU program, which proves the stability and accuracy of the method of the present invention.

[0100] In an exemplary embodiment, in order to verify the efficiency of the parallel acceleration method for two-phase color gradient LBM provided by the present invention, this embodiment compares the efficiency of the central processing unit (CPU) program with the current mainstream OpenMP parallel program. The test case is the droplet generation problem in the microchannel. In order to ensure the fairness of the comparison, the test hardware adopts the high-end configuration of the same period: two AMD EPYC 7763s are used on the CPU side, and a single NVIDIA A100 40GB is used on the GPU side. The efficiency test results are as follows: Fig.10 As shown, Fig.10 The figure 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 CPU program parallel to the current mainstream OpenMP. The test case is the droplet generation problem in a microchannel with a geometric size of 2475×367×1. When comparing, two AMD EPYC 7763s are used on the CPU side, and a single NVIDIA A100 40GB is used on the GPU side. Fig.10The speed on the vertical axis means the time (in seconds) required for each 10,000 iterations of the program. 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 parallel CPU program.

[0101] In an exemplary embodiment, in order to verify the parallelism of the parallel acceleration method applicable to the two-phase color gradient LBM, that is, the efficiency when expanded to more GPUs, this embodiment adopts the example in the first embodiment, and verifies the acceleration efficiency under different numbers of GPUs by increasing the number of GPUs involved in the calculation. The results are as follows: Fig.11 As shown, Fig.11 This is a parallel efficiency diagram of the parallel acceleration method for the 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 large number of gradient solutions in the two-phase flow, the cross-device transmission of some data is involved, resulting in a certain distance 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, reflecting good parallelism.

[0102] In an exemplary embodiment, this embodiment is used to illustrate the application scenario of the present invention. Droplet microfluidics is a technology for manipulating the precise movement and processing of droplets at the micrometer scale. It has a wide range of applications in the fields of biomedicine, chemical synthesis, and materials science, and its advantages are high efficiency, low consumption, and multifunctional operation. Numerical simulation of droplet microfluidics is a key tool for understanding and optimizing droplet behavior, and can provide detailed information on complex flows and interface dynamics in microfluidic systems. However, the main difficulties currently faced by such numerical simulations include: the complexity of interface effects and multi-physical field coupling, accurate capture of nonlinear flow behavior, and high computational costs of implementing numerical simulations under actual microchannel geometry. Since the present invention is applicable to complex two-phase flow simulations including flow in complex geometries and foam flow, and can achieve efficient calculation of large grid size problems on multi-GPU platforms, it has become an efficient numerical simulation tool for droplet microfluidics.

[0103] This embodiment demonstrates the numerical simulation of different arrangements of droplets spontaneously formed during the injection process in a three-dimensional microfluidic channel. This behavior of droplets spontaneously forming different arrangements during the injection process has been observed in earlier experiments, but previous numerical simulation methods have been difficult to accurately simulate this behavior due to problems with calculation scale and algorithm accuracy. The simulation results of droplet generation in the microfluidic channel of the present invention are compared with previous experiments and numerical simulations. Fig.12 It can be seen that the simulation accuracy and results of the present invention far exceed the previous simulation results, and are basically consistent with the experiment. At the same time, it also provides complete flow field information that cannot be obtained in the experiment, which greatly facilitates further research on the flow mechanism.

[0104] 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: (1) Industrial applications: used to study multiphase flow in oil and gas production, gas-liquid two-phase flow in fuel cells, and droplet dynamics in inkjet printing.

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

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

[0107] (4) Materials Science: Used to design new composite materials and optimize the material microstructure by simulating processes such as liquid phase sintering and phase separation.

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

[0109] 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 Figure 1 The steps are executed in the order shown. The specific execution order of the steps can be determined according to needs, and the present invention does not limit this.

[0110] The above is a parallel acceleration method for a 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 a two-phase color gradient lattice Boltzmann method, such as Fig.12 shown.

[0111] Fig.13 A schematic diagram of a parallel acceleration device for a two-phase color gradient lattice Boltzmann method provided by the present invention, the device 1300 comprises: The storage module 1301 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 boundary solid points, and constructing a storage information array; the storage information array includes the storage location of each fluid point and each boundary solid point and the mapping relationship with the points in the adjacent direction; A segmentation module 1302 is used to divide the porous medium into regions to obtain multiple flow field regions; The processing module 1303 is used to allocate a graphics processor to each flow field area, and 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.

[0112] For the specific definition of the parallel acceleration device applicable to the two-phase color gradient lattice Boltzmann method, please refer to the definition of the parallel acceleration method applicable to the two-phase color gradient lattice Boltzmann method above, which will not be repeated here. Each module in the above-mentioned parallel acceleration device applicable to the two-phase color gradient lattice Boltzmann method can be implemented in whole or in part by software, hardware and a combination thereof. The above-mentioned modules can be embedded in or independent of the processor in the computer device in the form of hardware, or can be stored in the memory of the computer device in the form of software, so that the processor can call and execute the operations corresponding to the above modules.

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

[0114] The present invention also provides Fig.14 The structural diagram of the computer device shown in FIG. Fig.14 As shown in the figure, at the hardware level, the computer device includes a processor, an internal bus, a network interface, a memory, and a non-volatile memory, and may also include other hardware required for the business. The processor reads the corresponding computer program from the non-volatile memory into the memory and then runs it to achieve the above Figure 1 A parallel acceleration method for the two-phase color gradient lattice Boltzmann method is provided.

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

[0116] The technical features of the above embodiments may be arbitrarily combined. To make the description concise, 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, they should be considered to be within the scope of 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.

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 parallel numerical simulation of each flow field region of the porous medium using the color gradient lattice Boltzmann method comprises: 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.

5. The method according to claim 4, 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.

6. The method according to claim 5, 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.

7. The method according to claim 5, characterized in that The calculation formula of the local interface curvature is: ； in, represents the local interface curvature, , and Denote the unit interface normal vector exist x , y , z Direction component.

8. The method according to claim 5, characterized in that The calculation formula of the force term is: ； ； ； in, represents the force term, Indicates speed, represents the multi-relaxation transformation matrix, express The inverse, represents the diagonal relaxation matrix, Indicates along Directional weight coefficient; Indicates along The discrete velocity vector in the direction, represents the surface tension, represents the repulsive force between bubbles when foam flow is realized, represents the lattice sound speed, represents the grid time step, represents the surface tension coefficient, represents the local interface curvature, Represents a color gradient.

9. The method according to claim 5, characterized in that The collision process is: ； in, represents the position vector of the fluid point; express After the moment of collision in space The discrete velocity along Total distribution function of direction; express Before the moment of collision in space The discrete velocity along Total distribution function of direction; express Collision operator in direction; Represents the external force term; the matrix multiplication operation in the collision operator adopts sparse matrix multiplication operation.

10. 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; The processing module is used to allocate a graphics processor to each flow field area, and 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.

Images