Porous rock multiphase fluid lattice Boltzmann simulation method based on GPU acceleration

Abstract

The invention discloses a porous rock multiphase fluid lattice Boltzmann simulation method based on GPU acceleration, and belongs to the technical field of computational fluid mechanics, a D3019 speed discrete model is utilized to identify porous rock to obtain a digital core model, a speed direction and a weight coefficient are defined, a distribution function is defined by means of an LBM method, and a multi-phase fluid lattice Boltzmann simulation model is obtained. Calculating a macroscopic physical quantity by utilizing a relation between a distribution function and a velocity moment, respectively describing a gas-liquid interface and a complex boundary of the porous rock by virtue of a scalar tracer function and a step function, and constructing a physical model; constructing an equilibrium state distribution function, performing LB equation calculation by adopting a forming schema equation under an external force, constructing a pseudo-potential multiphase model by virtue of a Shang-Chen potential energy formula, and processing the solid-liquid boundary of the porous rock by using a semi-rebound format; and carrying out a numerical solution pseudo code process based on the lattice Boltzmann method. According to the method, the characteristics of multithreading, memory sharing and the like of GPU hardware are fully utilized under the LBM framework, and the steps of the core algorithm are highly parallelized, so that high-performance calculation is realized.

Description

Technical Field

[0001] The invention belongs to the technical field of computational fluid dynamics, and in particular relates to a lattice Boltzmann simulation method of porous rock multiphase fluid based on GPU acceleration. Background Art

[0002] The seepage process of multiphase fluids in porous media often occurs in nature and industrial production practices, such as groundwater, oil and gas reservoir exploitation, and geological storage of carbon dioxide. Accurately and efficiently simulating the migration laws of multiphase fluids in rock pores has become the key to deeply understand the internal mechanism of multiphase fluid permeation process and guide engineering development practices.

[0003] Traditional multiphase fluid numerical simulation methods are mainly divided into two categories: one is the continuity method, such as finite difference, finite element, etc., which obtains flow field information by solving macroscopic control equations (such as Navier-Stokes equations); the other is particle methods, such as molecular dynamics, dissipative particle dynamics, etc., which simulate fluid behavior by tracking the movement of a large number of particles. These two types of methods have their own advantages and disadvantages: the continuity method is suitable for the simulation of macroscopic and engineering scale problems, but it is difficult to accurately describe complex pore structures; the particle method can characterize complex flow phenomena from a microscopic mechanism, but as the number of particles increases sharply, its computational cost also rises sharply, making it difficult to apply to actual engineering problems.

[0004] In recent years, the lattice Boltzmann method (LBM) has attracted more and more attention as a mesoscopic fluid simulation method. LBM solves the Boltzmann equation in discrete velocity space, establishes the connection between macroscopic physical quantities and mesoscopic distribution functions, and describes the evolution of fluids through local collisions and migration between grid points. Compared with traditional methods, LBM has the advantages of clear physical concepts, concise mathematical forms, flexible boundary processing, and easy parallelization. It is especially suitable for flow simulation in complex structures such as porous rocks. At present, LBM has been successfully applied to the simulation of single-phase flow, laminar flow, turbulent flow and other fields in porous rocks. However, there are still many challenges in extending the application of LBM to actual gas-liquid two-phase flow in porous rocks: (1) The pore structure of porous rocks is complex and the pore scale spans a large range, which places high demands on the geometric resolution and flexibility of the calculation method; (2) The interfacial tension effect, contact angle and other physical and chemical effects between the gas and liquid phases have an important influence on the interfacial morphology and migration, and an appropriate multiphase model needs to be introduced to describe them; (3) The actual displacement process often involves millions or even hundreds of millions of pores, and the corresponding calculation scale is huge, which places strict requirements on the execution efficiency of the algorithm and the computing power of the hardware. Summary of the invention

[0005] In view of the shortcomings of the background technology, the purpose of the present invention is to propose a lattice Boltzmann simulation method of porous rock multiphase fluid based on GPU acceleration, comprising the following steps: Step S1, using the D3019 velocity discrete model to identify porous rocks to obtain a digital core model, defining the velocity direction and weight coefficient on the model, defining two distribution functions with the help of the LBM method, calculating macroscopic physical quantities using the connection between the distribution function and the velocity moment, and using the scalar tracer function and the step function to describe the complex boundaries of the gas-liquid interface and porous rocks, respectively, to construct a physical model; Step S2, using the macroscopic physical quantities in the physical model to construct the equilibrium distribution function, using the forcing scheme equation to calculate the LB equation describing the fluid evolution under the action of external force, constructing a pseudo-potential multiphase model with the help of the Shan-Chen potential energy formula, and then using the semi-rebound format to process the solid-liquid boundary of the porous rock, so as to accurately simulate the dynamic fluid behavior in the porous rock; Step S3, a pseudo-code process of numerical solution based on the lattice Boltzmann method.

[0006] Preferably, the specific process of step S1 includes: 1) LBM method LBM is achieved by defining two independent distribution functions f i and g i Describe the velocity distribution of the gas phase and the liquid phase respectively, where the subscript i Indicates different discrete speed directions. The D3Q19 speed discrete model defines 19 speed directions e i And the corresponding weight coefficient w i : The calculation formula for a stationary particle with zero velocity in discrete velocity space is:

[0007] The calculation formula for particles in the main axis direction is:

[0008] The calculation formula for diagonal particles is:

[0009] The velocity moment is obtained by weighted summation of the distribution function, and the macroscopic physical quantity can be calculated by the velocity moment. The calculation formula of the macroscopic physical quantity of the relationship between the velocity moment and the distribution function is:

[0010] in, is the gas phase density, is the liquid density, u is the macroscopic velocity field, is the density of the macroscopic velocity field; 2) Describing the location of gas-liquid interfaces and complex boundaries of porous rocks Introducing scalar tracer functions , used to describe the position of the gas-liquid interface, ψ>0 represents the gas phase region, ψ<0 represents the liquid phase region, and ψ=0 is the interface; the calculation formula for the relationship between ψ and the macroscopic density field is:

[0011] in, is the mixed density; The complex boundary of porous rock is described by a step function B(x), which can clearly separate the solid area and the fluid area and be applied to various porous rock models: B=1 represents the solid area, B=0 represents the fluid area, and the wettability is introduced through the dynamic contact angle model.

[0012] Preferably, the specific process in step S2 is as follows: 1) LB equation calculation Constructing equilibrium distribution function through macroscopic physical quantities f i eq , g i eq The calculation formula is:

[0013] In the formula f i eq and g i eq is the equilibrium distribution function of the gas phase and the liquid phase, and u is the fluid velocity; Relaxation time t The kinematic viscosity of the fluid v The relationship is calculated as:

[0014] Among them, c s is the lattice sound speed; External forces are introduced into LBM to simulate the influence of external force effects on gas and liquid phases in fluid dynamics. The forcing scheme equation of Guo is:

[0015] Among them, Ftotal and Gtotal are the total external forces on the gas phase and liquid phase respectively.F i and G i is the external force, including gravity or other body forces; The LB equation is used to simulate the movement and collision of particles on a discrete grid to simulate the evolution of the fluid. The LB equation is:

[0016] in, t f and t g is the corresponding relaxation time, which determines the viscous properties of the fluid; Δt is the time step ,Δx is the spatial step length; 2) Constructing a pseudo-potential multiphase model A stable equilibrium state achieved without external forces ψ eq The calculation formula is:

[0017] In the formula A , B is a constant, controlling the gas-liquid density ratio, Control the surface tension. ψ eq is the equilibrium state of the interface; Surface tension term describing the surface tension effect at the interface Given by Shan-Chen power, the calculation formula of Shan-Chen power is:

[0018] Among them, G controls the strength of surface tension, is the surface tension term; The evolution of the phase interface is described by the improved pseudo-potential multiphase model. The calculation formula of the pseudo-potential multiphase model is:

[0019] in, is the interface relaxation time, with a value between 0.7 and 1.5; 3) Describing the solid-liquid boundary of porous rocks The solid-liquid boundary of porous rock is treated by semi-rebound format; for nodes close to the solid wall, x f , and its corresponding direction i The distribution function formula on is:

[0020] in, jfor i The reverse of e j =- e i , u w is the wall velocity; The format is simple and effective, and can be easily implemented in parallel.

[0021] Preferably, the specific process in step S24 is as follows: Discretizing the above equations, we get an explicit LBM evolution format. The numerical solution process of LBM can be divided into two main steps: collision and migration. With the help of computer control system, the core structure of the whole algorithm is as follows: for n = 0 : nstep do for all lattice x do Calculate the macro quantities p_r^n(x), p_b^n(x), u^n(x), ^n(x) Construct equilibrium states f_i^{eq,n}, g_i^{eq,n}, ^{eq,n} Calculate the external force terms F_i^n, G_i^n, S_ ^n Collision: Update f_i and g_i renew

[0022] Migration: pass the updated f_i and g_i to the adjacent nodes end for end for The main body of the algorithm is a two-layer loop: the outer loop traverses the time step, and the inner loop traverses all grid points; each grid point performs macroscopic quantity calculation, equilibrium state construction, collision relaxation, and migration update in turn, thereby completing the evolution of each step; the solid-liquid boundary is processed using a semi-rebound format.

[0023] Preferably, the specific steps of the GPU parallel acceleration operation based on building a physical model and LBM method are as follows: Adopt the CUDA programming model: divide the LBM grid based on the D3Q19 model into several CUDA blocks, each block is responsible for the parallel calculation of a subdomain; inside the block, each CUDA thread is responsible for the evolution of a grid point; frequently accessed data is cached through shared memory to reduce time-consuming global memory access; when necessary, use the synchronization instruction __syncthreads() to ensure data consistency.

[0024] In terms of multi-GPU parallelism, a regional decomposition strategy is adopted: the computational domain is divided into several sub-areas and assigned to different GPUs respectively, so as to make full use of the computing resources of multiple GPUs; coupling is achieved through boundary data exchange to ensure data consistency between different GPUs; the communication overhead is concealed through the pipeline strategy, and the transmission and calculation are overlapped, so that the GPU can continue to calculate while waiting for data transmission.

[0025] The present invention also proposes a gas-liquid two-phase flow parallel simulation system applied to the GPU-accelerated porous rock multiphase fluid lattice Boltzmann simulation method according to claim 1, characterized in that the main processes using modular design include: 1) Pre-processing module: import CT images or digital cores of porous rocks to reconstruct their internal pore structure; generate voxelized structured grids and reasonably divide parallel partitions; 2) Numerical solution module, which implements a multi-GPU parallel LBM solver based on CUDA; supports double-precision floating-point operations, sets parallel scale, partition strategy, etc.; 3) Boundary processing module, which realizes parallel processing of multiple boundary conditions such as wall, symmetry, period, pressure, velocity, etc.; supports complex situations such as moving boundary and wetting contact angle; 4) Post-processing module, parallel reading of instantaneous flow field, extraction of key physical quantities; development of parallel visualization tools, rendering of gas-liquid interface, velocity vector field, etc.; The above modules work together through standard interfaces such as MPI communication to realize the whole process of gas-liquid two-phase flow simulation; users input calculation parameters, monitor progress, and analyze results through a graphical interface to achieve human-computer interaction.

[0026] Compared with the prior art, the present invention proposes a lattice Boltzmann simulation method for porous rock multiphase fluid based on GPU acceleration, which has the following advantages: 1. Comprehensively considering complex factors such as porous structure, wettability, surface tension, external force drive, etc., a complete and sophisticated multiphase flow LBM physical model was established; 2. Make full use of the parallel characteristics of GPU hardware to optimize performance from multiple levels such as data structure, storage format, and communication strategy, and achieve efficient large-scale LBM simulation; 3. Developed a complete parallel multiphase flow simulation system, realizing full process integration including pre-processing, solution, and visualization. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] Figure 1 This is a digital core model diagram of porous rock; Figure 2 It is a physical model diagram; Figure 3A diagram showing the water-gas two-phase seepage process of porous rock of the present invention; Figure 4 This is a pore seepage path diagram of porous rocks according to the present invention. DETAILED DESCRIPTION

[0028] The following will be combined with the drawings in the embodiments of the present application to further clearly and completely describe the technical solutions in the embodiments of the present application. It should be noted that the described embodiments are only part of the embodiments of the present application, not all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without making creative work are within the scope of protection of this application.

[0029] In order to make the invention objectives, technical solutions and advantages of the present application clearer, the embodiments of the present application are further described in detail in conjunction with the drawings in the specification: In order to more clearly understand the above-mentioned objects, features and advantages of the present invention, the advantages of the present invention will be further illustrated by comparing the embodiments in conjunction with the drawings and specific implementation methods.

[0030] The present invention proposes a lattice Boltzmann simulation method for porous rock multiphase fluid based on GPU acceleration, and the steps of the method are described in detail: Step S1, using the D3019 velocity discrete model to identify porous rocks to obtain a digital core model (such as Figure 1 As shown in the figure, the velocity direction and weight coefficient are defined on the model. After defining two distribution functions with the help of the LBM method, the macroscopic physical quantities are calculated using the connection between the distribution function and the velocity moment. The complex boundaries of the gas-liquid interface and porous rock are described respectively with the help of scalar tracer function and step function, and then the physical model is constructed (as shown in the figure). Figure 2 shown); Furthermore, the specific process of step 1 is as follows: 1) LBM method LBM (Lattice Boltzmann Method) is a numerical method based on the distribution function of mesoscopic particles. It describes the dynamic behavior of fluids by discretizing velocity directions and positions. Specifically, LBM regards fluids as discrete systems composed of a large number of particles that move on a regular grid and interact with each other through collision rules. In the D3Q19 velocity discretization model (a discretization method for simulating three-dimensional fluid dynamics), there are 19 velocity directions in each grid, including one stationary particle, six particles in the main axis direction, and 12 particles in the diagonal direction. Each velocity direction has a weight coefficient to calculate the contribution of the particle. LBM defines two independent distribution functions f i and g iDescribe the velocity distribution of the gas phase and the liquid phase respectively, where the subscript i Indicates different discrete speed directions. The D3Q19 speed discrete model defines 19 speed directions e i And the corresponding weight coefficient w i : The calculation formula for a stationary particle with zero velocity in three-dimensional space is:

[0031] The particles in the main axis direction refer to the particles moving along the three main axis directions of x, y, and z in the coordinate system. The calculation formula for the particles in the main axis direction is:

[0032] Diagonal particles refer to particles moving along the diagonal directions of the xy, xz, and yz planes. The calculation formula for diagonal particles is:

[0033] The velocity moment is obtained by weighted summation of the distribution function, and the macroscopic physical quantity can be calculated by the velocity moment. The calculation formula of the macroscopic physical quantity of the relationship between the velocity moment and the distribution function is:

[0034] in, is the gas phase density, is the liquid density, u is the macroscopic velocity field, is the density of the macroscopic velocity field; 2) Describing the location of gas-liquid interfaces and complex boundaries of porous rocks Introducing scalar tracer functions , used to describe the position of the gas-liquid interface, ψ>0 represents the gas phase region, ψ<0 represents the liquid phase region, and ψ=0 is the interface; the calculation formula for the relationship between ψ and the macroscopic density field is:

[0035] in, is the mixed density; In order to clearly separate the solid and fluid regions and apply them to various porous rock models, the complex boundary of porous rocks (the interface between the pores inside the rock and the solid part) is described by a step function B(x): B=1 represents the solid region, B=0 represents the fluid region, and the wettability (the degree to which the liquid contacts and spreads on the solid surface) is introduced through a dynamic contact angle model (the smaller the contact angle, the better the wettability).

[0036] Step S2, using the macroscopic physical quantities in the physical model to construct the equilibrium distribution function, using the forcing scheme equation to calculate the LB equation describing the fluid evolution under the action of external force, constructing a pseudo-potential multiphase model with the help of the Shan-Chen potential energy formula, and then using the semi-rebound format to process the solid-liquid boundary of the porous rock, so as to more accurately simulate the dynamic fluid behavior in the porous rock; Furthermore, the specific process of step S2 is as follows: 1) LB equation Constructing equilibrium distribution function through macroscopic physical quantities f i eq , g i eq The calculation formula is:

[0037] In the formula f i eq and g i eq is the equilibrium distribution function of the gas phase and the liquid phase, and u is the fluid velocity; Relaxation time t The kinematic viscosity of the fluid v The relationship is calculated as:

[0038] Among them, c s is the lattice sound speed; In order to simulate the influence of external force effects on gas and liquid phases in fluid dynamics, external forces (such as gravity) are introduced into LBM, and the forcing scheme equation of Guo is adopted as follows:

[0039] Among them, Ftotal and Gtotal are the total external forces on the gas phase and liquid phase respectively. F i and G i is the external force, including gravity or other body forces; To simulate the movement and collision of particles on a discrete grid to simulate the evolution of the fluid, the LB (Boltzmann equation) equation is used:

[0040] in, t f and t gis the corresponding relaxation time, which determines the viscous properties of the fluid; Δt is the time step ,Δx is the spatial step length; 2) Constructing a pseudo-potential multiphase model A stable equilibrium state achieved without external forces ψ eq The calculation formula is:

[0041] In the formula, A and B are constants, controlling the gas-liquid density ratio. Control the surface tension. ψ eq is the equilibrium state of the interface; Surface tension term Given by the Shan-Chen force, the surface tension term The surface tension effect of the phase interface can be described. The calculation formula of the Shan-Chen force is:

[0042] Among them, G controls the strength of surface tension, is the surface tension term; The evolution of the phase interface [the interface between two different phases (such as liquid and gas, solid and liquid, etc.)] is described by the improved pseudo-potential multiphase model. The calculation formula of the pseudo-potential multiphase model is:

[0043] in, is the interface relaxation time, with a value between 0.7 and 1.5; 3) Describing the solid-liquid boundary of porous rocks The solid-liquid boundary of porous rock is processed by half-bounce-back format; for nodes close to the solid wall x f , and its corresponding direction i The distribution function formula on is:

[0044] in, j for i The reverse of e j =- e i , u w is the wall velocity; The format is simple and effective, and can be easily implemented in parallel.

[0045] Step S3, a numerical solution pseudocode process based on the lattice Boltzmann method can efficiently simulate the dynamic behavior of the fluid; Furthermore, the specific process of step S3 is as follows: Discretizing the above equations yields an explicit LBM evolution format. The numerical solution process of LBM can be divided into two main steps: collision and streaming. With the help of a computer control system, the core structure of the entire algorithm is as follows: for n = 0 : nstep do for all lattice x do Calculate the macro quantities p_r^n(x), p_b^n(x), u^n(x), ^n(x) Construct equilibrium states f_i^{eq,n}, g_i^{eq,n}, ^{eq,n} Calculate the external force terms F_i^n, G_i^n, S_ ^n Collision: Update f_i and g_i renew

[0046] Migration: pass the updated f_i and g_i to the adjacent nodes end for end for The main body of the algorithm is a two-layer loop: the outer loop traverses the time step, and the inner loop traverses all grid points; each grid point performs macroscopic quantity calculation, equilibrium state construction, collision relaxation, and migration update in turn, thereby completing the evolution of each step; the solid-liquid boundary is processed using a semi-rebound format.

[0047] Based on the construction of physical models and LBM methods, the CUDA programming model is used to perform GPU parallel acceleration operations to improve computing efficiency; Furthermore, the process of the GPU parallel acceleration operation in the step is as follows: The LBM algorithm is explicit and local, which is very suitable for the single instruction multiple threads (SIMT) parallel mode of GPU. The hardware features of GPU, such as multi-threading and on-chip shared memory, can be fully utilized to highly parallelize the core steps of LBM and achieve high-performance computing. In the specific implementation, the CUDA programming model is adopted; the LBM grid based on the D3Q19 model is divided into several CUDAblocks (basic execution units, a collection of CUDA threads), each block is responsible for the parallel calculation of a subdomain; inside the block, each CUDA thread (minimum execution unit, thread) is responsible for the evolution of a grid point; frequently accessed data is cached through shared memory (cache memory, located inside each block) to reduce time-consuming global memory access; when necessary, the synchronization instruction __syncthreads() (synchronization instruction) is used to ensure data consistency.

[0048] In terms of multi-GPU parallelism, a domain decomposition strategy is adopted; the computational domain is divided into several sub-areas and assigned to different GPUs respectively, so as to make full use of the computing resources of multiple GPUs; coupling is achieved through boundary data exchange to ensure data consistency between different GPUs; the communication overhead is concealed through a pipeline strategy, and transmission and computation are overlapped, so that the GPU can continue computing while waiting for data transmission.

[0049] Based on the above algorithm, a complete gas-liquid two-phase flow parallel simulation system was developed; the system adopts modular design; Furthermore, the main processes of the gas-liquid two-phase flow parallel simulation system include: 1) Pre-processing module: import CT images or digital cores of porous rocks to reconstruct their internal pore structure; generate voxelized structured grids and reasonably divide parallel partitions; 2) Numerical solution module, which implements a multi-GPU parallel LBM solver based on CUDA; supports double-precision floating-point operations, and can set parallel scale, partition strategy, etc.; 3) Boundary processing module, which realizes the parallel processing of multiple boundary conditions such as wall, symmetry, period, pressure, velocity, etc.; supports complex situations such as moving boundary and wetting contact angle; 4) Post-processing module, parallel reading of instantaneous flow field, extraction of key physical quantities; developed parallel visualization tools, rendering gas-liquid interface, velocity vector field, etc.; The above modules work together through standard interfaces such as MPI communication to realize the whole process of gas-liquid two-phase flow simulation; users input calculation parameters, monitor progress, and analyze results through a graphical interface to achieve human-computer interaction.

[0050] like Figure 3 and Figure 4 As shown, the system embodiment of the present invention: 1. Initialization and parameter setting The present invention first imports necessary Python libraries, such as numpy for scientific computing, taichi for GPU acceleration, and pyevtk for VTK format output.

[0051] Define the key parameters of the program: The grid size nx, ny, nz defines the discretization degree of the three-dimensional computational domain and is set to 232^3. A larger grid means a finer but more expensive simulation. The external forces fx, fy, and fz define the volume forces that drive the flow, taking the direction of gravity and adding a small driving force in the x direction; niu_l and niu_g are the kinematic viscosities of the two fluids, which determine the damping characteristics of the flow. Considering common low-viscosity fluids, the value is 0.1 for liquids and a smaller value for gases; psi_solid defines the wettability of the solid, 0.7 is hydrophilic, making the liquid phase tend to adhere to the solid wall; CapA defines the strength of surface tension, with 0.005 being a suitable dimensionless value; The velocity and pressure are applied to the left and right boundaries in the x direction respectively, and the boundary conditions are set as follows: The left boundary fixed injection velocity vx_bcxl=0.001, corresponding to a fixed liquid phase flow rate; The fixed pressure on the right boundary is rho_bcxr=0.995 to ensure that the outlet is in the gas phase; The remaining four boundaries take periodic conditions to simulate infinitely extended porous rocks; The wettability boundary psi_x_left=1.0 makes the left boundary initially filled with liquid phase, driving the liquid phase to seep to the right and expel the gas phase; the right boundary psi_x_right=-1.0 is initially gas phase; the rest are periodic boundaries.

[0052] 2. Field variable storage LBM involves a large number of field variables and needs to utilize limited GPU memory; The distribution function f[i,j,k,l] is a four-dimensional array, (i,j,k) is the spatial index, l is the discrete speed index, ranging from 0 to 18, considering double-precision floating points, 19 speed directions, which is the most memory-consuming variable; F is the buffer of the distribution function f, which is used to store the intermediate state after collision and before propagation, and is of the same type as f; Macroscopic physical quantities such as density rho[i,j,k] and velocity v[i,j,k] are three-dimensional arrays with fewer components than f and less overhead; The density of the gas and liquid phases, rho_r, rho_b, and the region identifier psi are three-dimensional arrays, occupying only 1 double space; To save memory, the solid distribution solid[i,j,k] is marked as a Boolean array, which only takes up 1 byte; Next, define some basic variables of the LBM model: e[l]=[(0,0,0),(±1,0,0)...] is the discrete velocity direction, taking the D3Q19 model in 19 directions; e_f is the floating point version of the velocity direction, used in macro velocity calculations; w[l]=[1 / 3, 1 / 18, ...] is the isotropic Gaussian integral weight to ensure isotropy; LR[l] records the reverse mapping of 18 velocity directions for the bounce-back format; Finally, we define the matrix M[19,19] and its inverse inv_M[19,19] in the MRT collision operator to link the macroscopic quantity with the moment space. The derivation of these matrices requires mathematical techniques such as the Chapman-Enskog expansion, which is the theoretical basis of LBM and will not be expanded here. At the end of the parameter setting, the relaxation coefficient of the two-phase fluid is calculated: wl, wg are calculated based on the kinematic viscosities niu_l and niu_g of the liquid and gas phases, respectively, to ensure accurate macroscopic fluid dynamics properties.

[0053] 3. Initialization After the field variable is defined, it needs to be assigned a reasonable initial value, which can be done through the init() function: Traverse the non-solid area (solid==0) and initialize its physical quantities: psi is initialized to -1, representing the gas phase, which matches the wettability of psi_solid=0.7; The initial values ​​of density rho, velocity v, and other macroscopic quantities are determined based on the hydrostatic relations; The density distributions of the two phases, rho_r and rho_b, are assigned values ​​according to the positive and negative values ​​of psi; Finally, the initial value of the distribution function f is calculated using the corresponding equilibrium state function feq to ensure that the initial field matches the macroscopic quantity; Buffer F copies the value of f as the initial state; The microstructure of porous rock, such as the distribution of solid, is loaded from an external file through init_geo(); in this example, the pre-generated geometry is read from a *.txt file in the form of 0 / 1 and assigned to solid; similarly, the initial distribution of the two phases, psi, can also be imported from an external file; this enables flexible initialization of complex porous structures.

[0054] 4. Main loop After initialization, the program enters the main loop and repeats the classic LBM steps of "collision-propagation"; At each time step: ① colission() performs "collision", that is, relaxes the distribution function in momentum space: Traverse each fluid node (i, j, k) and complete the following steps: The macroscopic density rho, velocity u, etc. at the node (i,j,k) are calculated by f[i,j,k,l]; From the macroscopic quantity, calculate the equilibrium distribution function f_eq[l; Introduce external forces, such as gravity and buoyancy, and store them as F[l]; By using the MRT matrix M[l,m], f[l] and f_eq[l] are transformed into moment space; For the non-conserved quantities in the moment space, a BGK type relaxation is performed with the relaxation coefficient w; Transform back to velocity space and superimpose the external force term to obtain the distribution function after the collision; Anisotropic relaxation ensures the physical properties of the macroscopic fluid, such as incompressibility; Store the distribution function after the collision to F[i,j,k,l] to complete the "collision"; ② streaming1() performs "propagation", that is, translating the distribution function in physical space: Traverse all nodes (i, j, k) and translate F according to e[l]; Specifically, in the degree direction l, F[i,j,k,l] propagates to its neighboring point F[i+ex,j+ey,k+ez,l]; This step is performed synchronously throughout the flow field, simulating the free flight of fluid particles; After propagation, the distribution function at the solid boundary is updated to reflect the no-slip condition; ③ Boundary_condition() handles boundary conditions: For periodic boundaries, the propagation of the distribution function f is continuous and no special treatment is required; For velocity boundaries, such as the left boundary, reset f by feq_bcxl; this equilibrium distribution corresponds to the injection velocity vx_bcxl; For pressure boundaries, such as the right boundary, it is also reset with feq_bcxr, corresponding to the outlet density rho_bcxr; For the solid-liquid boundary, a semi-rebound format is used to simulate the viscous no-slip condition, which is achieved by finding the reverse velocity through LR[l]; Wettability is introduced by the contact angle dynamics equation as an additional term superimposed to control the dynamic contact line; ④ streaming3() copies the propagated distribution function F back to f and calculates the macro quantity: First, copy the entire flow field F to f to prepare for the next step of "collision"; At the same time, the macroscopic density rho[i,j,k] and velocity u[i,j,k] at the node (i,j,k) are calculated by using the moment relationship based on the updated f; Update the interface scalar psi by the density fields rho_r, rho_b; At this point, one time step of fluid iteration is completed; repeat multiple time steps to simulate the dynamic evolution of the flow; ⑤ Output the physical quantities of the flow field to a file for post-processing analysis and visualization.

[0055] At certain intervals, such as 100 steps, the entire flow field is outputted through I / O; Using the pyevtk library, scalar and vector fields are organized into *.vtr format, which can be read by software such as Paraview; The entire process is encapsulated in a high-level Python script, which can flexibly set parameters and control the start and stop of the simulation. It is compiled and executed on the GPU through CUDA, giving full play to the advantages of parallel computing and completing simulations of millions of grids in a few hours.

[0056] Although the preferred embodiments of the present application have been described, those skilled in the art may make additional changes and modifications to these embodiments once they have learned the basic creative concept. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments and all changes and modifications falling within the scope of the present application.

[0057] Obviously, those skilled in the art can make various changes and modifications to the present application without departing from the spirit and scope of the present application. Thus, if these modifications and variations of the present application fall within the scope of the claims of the present application and their equivalents, the present application is also intended to include these modifications and variations.

Claims

1. A GPU-accelerated lattice Boltzmann simulation method for porous rock multiphase fluid, characterized in that: include: Step S1, using the D3019 velocity discrete model to identify porous rocks to obtain a digital core model, defining the velocity direction and weight coefficient on the model, defining two distribution functions with the help of the LBM method, calculating macroscopic physical quantities using the connection between the distribution function and the velocity moment, and using the scalar tracer function and the step function to describe the complex boundaries of the gas-liquid interface and porous rocks, respectively, to construct a physical model; Step S2, using the macroscopic physical quantities in the physical model to construct the equilibrium distribution function, using the forcing scheme equation to calculate the LB equation describing the fluid evolution under the action of external force, constructing a pseudo-potential multiphase model with the help of the Shan-Chen potential energy formula, and then using the semi-rebound format to process the solid-liquid boundary of the porous rock; Step S3, a pseudo-code process of numerical solution based on the lattice Boltzmann method.

2. The method for lattice Boltzmann simulation of porous rock multiphase fluid based on GPU acceleration according to claim 1 is characterized in that: The specific process of step S1 is as follows: 1) LBM method LBM defines two independent distribution functions f i and g i Describe the velocity distribution of the gas phase and the liquid phase respectively, where the subscript i Indicates different discrete speed directions. The D3Q19 speed discrete model defines 19 speed directions. e i And the corresponding weight coefficient w i : The calculation formula for a stationary particle with zero velocity in discrete velocity space is: ； The calculation formula for particles in the main axis direction is: ； The calculation formula for diagonal particles is: ； The velocity moment is obtained by weighted summation of the distribution function, and the macroscopic physical quantity can be calculated by the velocity moment. The calculation formula of the macroscopic physical quantity of the relationship between the velocity moment and the distribution function is: ； in, is the gas phase density, is the liquid density, u is the macroscopic velocity field, is the density of the macroscopic velocity field; 2) Describe the location of the gas-liquid interface and the boundaries of porous rocks The position of the gas-liquid interface is traced by a scalar tracer function Description, where ψ>0 represents the gas phase region, ψ<0 represents the liquid phase region, and ψ=0 is the interface; the calculation formula for the relationship between ψ and the macroscopic density field is: ； in, is the mixed density; The boundary of porous rock is described by a step function B(x), which separates the solid region and the fluid region and applies it to various porous rock models: B=1 represents the solid region, B=0 represents the fluid region, and wettability is introduced through the dynamic contact angle model.

3. The method for lattice Boltzmann simulation of porous rock multiphase fluid based on GPU acceleration according to claim 1 is characterized in that: The specific process of step S2 is as follows: 1) LB equation calculation Constructing equilibrium distribution function through macroscopic physical quantities f i eq , g i eq The calculation formula is: ； In the formula f i eq and g i eq is the equilibrium distribution function of the gas phase and the liquid phase, and u is the fluid velocity; Relaxation time τ The kinematic viscosity of the fluid v The relationship is calculated as: ； Among them, c s is the lattice sound speed; Introducing external forces into LBM to simulate the influence of external forces in fluid dynamics on the gas and liquid phases, the forcing scheme equation of Guo is: ； Among them, Ftotal and Gtotal are the total external forces on the gas phase and liquid phase respectively. F i and G i is the external force, including gravity or other body forces; The LB equation is used to simulate the movement and collision of particles on a discrete grid to simulate the evolution of the fluid. The LB equation is: ； in, τ f and τ g is the corresponding relaxation time, which determines the viscous properties of the fluid; Δt is the time step , Δx is the spatial step length; 2) Constructing a pseudo-potential multiphase model A stable equilibrium state achieved without external forces ψ eq The calculation formula is: ； In the formula, A and B are constants, which control the gas-liquid density ratio. Control the surface tension. ψ eq is the equilibrium state of the interface; Surface tension term describing the surface tension effect at the interface Given by Shan-Chen power, the calculation formula of Shan-Chen power is: ； Among them, G controls the strength of surface tension, is the surface tension term; The evolution of the phase interface is described by the improved pseudo-potential multiphase model. The calculation formula of the pseudo-potential multiphase model is: ； in, is the interface relaxation time, with a value between 0.7 and 1.5; 3) Describing the solid-liquid boundary of porous rocks The solid-liquid boundary of porous rock is treated by semi-rebound format; for nodes close to the solid wall, x f , and its corresponding direction i The distribution function formula on is: ； in, j for i The reverse of e j =- e i , u w is the wall velocity; The format is simple and effective, and can be easily implemented in parallel.

4. The method for lattice Boltzmann simulation of porous rock multiphase fluid based on GPU acceleration according to claim 1 is characterized in that: The specific process in step S3 is as follows: Discretizing the above equations yields an explicit LBM evolution format. The numerical solution process of LBM can be divided into two main steps: collision and migration. Its core structure is as follows: for n = 0 : nstep do for all lattice x do Calculate the macro quantities p_r^n(x), p_b^n(x), u^n(x), ^n(x) Construct equilibrium states f_i^{eq,n}, g_i^{eq,n}, ^{eq,n} Calculate the external force terms F_i^n, G_i^n, S_ ^n Collision: Update f_i and g_i renew , Migration: pass the updated f_i and g_i to the adjacent nodes end for end for The main body of this core structure is a two-layer loop: the outer loop traverses the time step, and the inner loop traverses all grid points; each grid point performs macroscopic quantity calculation, equilibrium state construction, collision relaxation, and migration update in turn, thereby completing the evolution of each step; the solid-liquid boundary is processed using a semi-rebound format.

5. The method for lattice Boltzmann simulation of porous rock multiphase fluid based on GPU acceleration according to claim 1, characterized in that: The physical model and LBM method are constructed using GPU parallel acceleration operations, and the specific steps are as follows: Adopt CUDA programming model: divide the LBM grid based on D3Q19 model into several CUDA blocks, each block is responsible for parallel calculation of a subdomain; inside the block, each CUDA thread is responsible for the evolution of a grid point; frequently accessed data is cached through shared memory to reduce time-consuming global memory access; use synchronization instruction __syncthreads() to ensure data consistency when necessary; In terms of multi-GPU parallelism, a regional decomposition strategy is adopted: the computational domain is divided into several sub-areas, which are assigned to different GPUs respectively to utilize the computing resources of multiple GPUs; coupling is achieved through boundary data exchange to ensure data consistency between different GPUs; the communication overhead is concealed through a pipeline strategy, and transmission and computation are overlapped so that the GPU can continue computing while waiting for data transmission.

6. A gas-liquid two-phase flow parallel simulation system applied to the GPU-accelerated porous rock multiphase fluid lattice Boltzmann simulation method according to claim 1, characterized in that: The main processes of adopting modular design include: 1) Pre-processing module: import CT images or digital cores of porous rocks to reconstruct their internal pore structure; generate voxelized structured grids and reasonably divide parallel partitions; 2) Numerical solution module, which implements a multi-GPU parallel LBM solver based on CUDA; supports double-precision floating-point operations, sets parallel scale, partition strategy, etc.; 3) Boundary processing module, which realizes parallel processing of multiple boundary conditions such as wall, symmetry, period, pressure, velocity, etc.; supports complex situations such as moving boundary and wetting contact angle; 4) Post-processing module, parallel reading of instantaneous flow field, extraction of key physical quantities; development of parallel visualization tools, rendering of gas-liquid interface, velocity vector field, etc.; The above modules work together through standard interfaces such as MPI communication to realize the whole process of gas-liquid two-phase flow simulation; users input calculation parameters, monitor progress, and analyze results through a graphical interface to achieve human-computer interaction.

Images