Gpu-accelerated lattice boltzmann method for multiphase flow in porous media

Abstract

The application discloses a GPU-accelerated porous rock multiphase fluid lattice Boltzmann simulation method, and belongs to the technical field of computational fluid dynamics. A digital core model is obtained by using a D3019 velocity discrete model to identify the porous rock, a velocity direction and a weight coefficient are defined, a distribution function is defined by means of an LBM method, macroscopic physical quantities are calculated by using the relationship between the distribution function and the velocity moment, a physical model is constructed by means of a scalar tracer function and a step function to respectively describe a gas-liquid interface and a complex boundary of the porous rock; an equilibrium distribution function is constructed, LB equation calculation is carried out under external force by using a forcing scheme equation, a pseudo-potential multiphase model is constructed by means of a Shan-Chen potential energy formula, and a semi-reflective bounce format is used to process a solid-liquid boundary of the porous rock; and a numerical solution pseudo-code process based on the lattice Boltzmann method is carried out. The application fully utilizes the multi-thread, shared memory and other characteristics of the GPU hardware under the LBM framework, highly parallelizes the core algorithm steps, and thus realizes high-performance calculation.

Description

Technical Field

[0001] This invention belongs to the field of computational fluid dynamics technology, specifically relating to a GPU-accelerated lattice Boltzmann simulation method for porous rock multiphase fluids. Background Technology

[0002] The seepage process of multiphase fluids in porous media is prevalent in nature and industrial production practices, such as groundwater and oil and gas reservoir extraction and geological carbon dioxide sequestration. Accurately and efficiently simulating the migration patterns of multiphase fluids in rock pores is crucial for a deeper understanding of the intrinsic mechanisms of multiphase fluid seepage processes and for guiding engineering development practices.

[0003] Traditional numerical simulation methods for multiphase fluids are mainly divided into two categories: continuous methods, such as finite difference and finite element methods, obtain flow field information by solving macroscopic governing equations (such as the Navier-Stokes equations); and particle methods, such as molecular dynamics and dissipative particle dynamics, simulate fluid behavior by tracking the motion of a large number of particles. Each of these methods has its advantages and disadvantages: continuous methods are suitable for simulating macroscopic, engineering-scale problems, but struggle to finely describe complex pore structures; particle methods can characterize complex flow phenomena from a microscopic mechanism perspective, but their computational cost increases dramatically with the number of particles, making them difficult to apply to practical engineering problems.

[0004] In recent years, the Lattice Boltzmann Method (LBM), as a mesoscale fluid simulation method, has attracted increasing attention. LBM establishes the relationship between macroscopic physical quantities and mesoscopic distribution functions by solving the Boltzmann equations in a discrete velocity space, describing fluid evolution through local collisions and inter-lattice migration. Compared with traditional methods, LBM has advantages such as clear physical concepts, concise mathematical form, flexible boundary treatment, and ease of parallelization, making it particularly suitable for flow simulation in complex structures such as porous rocks. Currently, LBM has been successfully applied to the simulation of single-phase flow, laminar flow, and turbulent flow in porous rocks. However, extending LBM to practical gas-liquid two-phase flow in porous rocks still faces many challenges:

[0005] (1) The pore structure of porous rocks is intricate and the pore scale spans a large range, which places high demands on the geometric resolution and flexibility of the calculation method.

[0006] (2) The interfacial tension effect and contact angle of the gas-liquid two-phase system have an important influence on the phase interface morphology and migration, and a suitable multiphase model needs to be introduced to describe them.

[0007] (3) The actual displacement process often involves millions or even hundreds of millions of pores, which corresponds to a huge computational scale, and puts forward stringent requirements on the execution efficiency of the algorithm and the computing power of the hardware. Summary of the Invention

[0008] To address the shortcomings of the prior art, the present invention aims to propose a GPU-accelerated lattice Boltzmann simulation method for porous rock multiphase fluids, comprising the following steps:

[0009] Step S1: Use the D3019 velocity discretization model to identify porous rocks and obtain a digital core model. Define the velocity direction and weight coefficients on the model. After defining two distribution functions using the LBM method, calculate macroscopic physical quantities using the relationship between the distribution functions and velocity moments. Use scalar tracer functions and step functions to describe the complex boundaries of the gas-liquid interface and porous rocks, respectively, and construct a physical model.

[0010] Step S2: Construct the equilibrium distribution function using macroscopic physical quantities in the physical model. Under the action of external forces, use the forcing scheme equation to calculate the LB equation describing fluid evolution. Construct a pseudo-potential multiphase model with the help of the Shan-Chen potential energy formula. Then use the semi-rebound scheme to process the solid-liquid boundary of porous rocks, thereby accurately simulating the dynamic fluid behavior in porous rocks.

[0011] Step S3: Pseudocode process for numerical solution based on the lattice Boltzmann method.

[0012] Preferably, the specific process of step S1 includes:

[0013] 1) LBM method

[0014] LBM is defined by two independent distribution functions. f i and g i Describe the velocity distributions in the gas and liquid phases respectively, where the subscripts are... i Representing different discrete velocity directions, the D3Q19 velocity discrete model defines 19 velocity directions. e i and corresponding weighting coefficients w i :

[0015] The formula for calculating a stationary particle with zero velocity in a discrete velocity space is:

[0016]

[0017] The formula for calculating particles along the principal axis is:

[0018]

[0019] The formula for calculating particles in the diagonal direction is:

[0020]

[0021] The velocity moment is obtained by weighted summation of the distribution functions. Macroscopic physical quantities can then be calculated from the velocity moment. The formula for calculating macroscopic physical quantities in relation to the velocity moment and the distribution function is as follows:

[0022]

[0023] in, The density is the gas phase density. The density of the liquid phase is... u For macroscopic velocity fields, The density of the macroscopic velocity field;

[0024] 2) Describe the location of the gas-liquid interface and the complex boundaries of porous rocks.

[0025] Introducing scalar tracer functions ψ is used to describe the location of the gas-liquid interface. ψ>0 represents the gas phase region, ψ<0 represents the liquid phase region, and ψ=0 is the interface. The formula for calculating the relationship between ψ and the macroscopic density field is:

[0026]

[0027] in, For mixed density;

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

[0029] Preferably, the specific process in step S2 is as follows:

[0030] 1) Calculation of LB equation

[0031] Constructing the equilibrium distribution function using macroscopic physical quantities f i eq , g i eq The calculation formula is:

[0032]

[0033] In the formula f i eq and g i eq It is the equilibrium distribution function of the gas and liquid phases, where u is the fluid velocity;

[0034] Relaxation timet With fluid kinematic viscosity v The formula for calculating the relationship is:

[0035]

[0036] Among them, c s For the speed of sound in a grid;

[0037] Introducing external forces into the LBM to simulate the effects of external forces on the gas and liquid phases in fluid dynamics, using Guo's forcing scheme equations:

[0038]

[0039] Where Ftotal and Gtotal are the total external forces acting on the gas phase and liquid phase, respectively. F i and G i These are external forces, including gravity or other volume forces.

[0040] The evolution of fluids is simulated by using the LB equations to model the motion and collisions of particles on a discrete lattice. The LB equations are as follows:

[0041]

[0042] in, t f and t g The corresponding relaxation time determines the viscous properties of the fluid; Δt For time step ,Δx The spatial step size;

[0043] 2) Constructing a pseudo-potential multiphase model

[0044] Equilibrium state that reaches a stable state without external forces ψ eq The calculation formula is:

[0045]

[0046] In the formula A , B The gas-liquid density ratio is a constant. Controlling the magnitude of surface tension, ψ eq This represents the equilibrium state of the interface.

[0047] Surface tension term describing the surface tension effect at the phase interface Provided by the Shan-Chen faction, the formula for calculating the Shan-Chen faction is as follows:

[0048]

[0049] Among them, G controls the surface tension intensity. This is the surface tension term;

[0050] The evolution of the phase interface is described by an improved pseudo-potential multiphase model, the calculation formula of which is:

[0051]

[0052] in, This is the interface relaxation time, with a value between 0.7 and 1.5.

[0053] 3) Describe the solid-liquid boundary of porous rocks

[0054] The solid-liquid boundary of porous rocks is treated using a semi-rebound scheme; for nodes near the solid wall... x f its corresponding direction i The formula for the distribution function on is:

[0055]

[0056] in, j for i The opposite, that is e j =- e i , u w The wall velocity;

[0057] This format is simple, effective, and easy to implement in parallel.

[0058] Preferably, the specific process in step S24 is as follows:

[0059] Discretizing the above equations yields an explicit LBM evolution scheme. The numerical solution process for LBM can be divided into two main steps: collision and migration. With the aid of a computer control system, the core structure of the entire algorithm is as follows:

[0060] for n = 0 : nstep do

[0061] for all lattice x do

[0062] Calculate the macroscopic quantities p_r^n(x), p_b^n(x), u^n(x). ^n(x)

[0063] Construct equilibrium states f_i^{eq,n}, g_i^{eq,n}, ^{eq,n}

[0064] Calculate the external force terms F_i^n, G_i^n, S_ ^n

[0065] Collision: Update f_i and g_i

[0066] renew

[0067] Migration: Pass the updated f_i and g_i to neighboring nodes.

[0068] end for

[0069] end for

[0070] The algorithm consists of two loops: the outer loop iterates through the time steps, and the inner loop iterates through all grid points. Each grid point sequentially performs macroscopic quantity calculations, equilibrium state construction, collision relaxation, and migration updates, thereby completing the evolution at each step. The solid-liquid boundary is handled using a semi-bounce scheme.

[0071] Preferably, the specific steps of the GPU parallel acceleration operation based on the physical model construction and LBM method are as follows:

[0072] The CUDA programming model is adopted: the LBM grid based on the D3Q19 model is divided into several CUDA blocks, each block is responsible for the parallel computation of a subdomain; within 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, the synchronization instruction __syncthreads() is used to ensure data consistency.

[0073] In terms of multi-GPU parallelism, a region decomposition strategy is adopted: the computing domain is divided into several sub-regions and assigned to different GPUs, which can make full use of the computing resources of multiple GPUs; coupling is achieved through boundary data exchange, which can ensure data consistency between different GPUs; and the communication overhead is masked by the pipeline strategy, which overlaps transmission and computation, so that the GPU can continue to perform computation while waiting for data transmission.

[0074] This invention also proposes a parallel simulation system for gas-liquid two-phase flow applied to the GPU-accelerated lattice Boltzmann simulation method for porous rock multiphase fluids as described in claim 1, characterized in that the main process of adopting a modular design includes:

[0075] 1) Preprocessing module: Import CT images or digital cores of porous rocks, reconstruct their internal pore structure; generate voxelized structured meshes, and rationally divide them into parallel partitions;

[0076] 2) Numerical solution module, implementing a CUDA-based multi-GPU parallel LBM solver; supports double-precision floating-point operations, and allows setting parallel scale, partitioning strategy, etc.

[0077] 3) Boundary processing module, which enables parallel processing of various boundary conditions such as wall, symmetry, periodicity, pressure, and velocity; supports complex cases such as moving boundaries and wetting contact angle;

[0078] 4) Post-processing module: reads the instantaneous flow field in parallel, extracts key physical quantities; develops parallel visualization tools to render gas-liquid interfaces, velocity vector fields, etc.

[0079] All the above modules work together through standard interfaces such as MPI communication to realize the entire process of gas-liquid two-phase flow simulation; users can input calculation parameters, monitor progress, and analyze results through a graphical interface to achieve human-computer interaction.

[0080] Compared with existing technologies, this invention proposes a GPU-accelerated lattice Boltzmann simulation method for porous rock multiphase fluids. The advantages of this method are:

[0081] 1. Taking into account complex factors such as porous structure, wettability, surface tension, and external force driving, a complete and detailed multiphase flow LBM physical model was established;

[0082] 2. By fully utilizing the parallel characteristics of GPU hardware, performance was optimized from multiple levels, including data structures, storage formats, and communication strategies, achieving efficient large-scale LBM simulation;

[0083] 3. A complete parallel multiphase flow simulation system was developed, realizing the integration of the entire process, including preprocessing, solving, and visualization. Attached Figure Description

[0084] Figure 1 A digital core model of porous rock;

[0085] Figure 2 Physical model diagram;

[0086] Figure 3 This is a diagram illustrating the two-phase flow process of water and air in the porous rock of this invention.

[0087] Figure 4 This is a diagram showing the seepage path in the porous rock pores of the present invention. Detailed Implementation

[0088] The technical solutions of the embodiments of this application will be further described clearly and completely below with reference to the accompanying drawings. It should be noted that the described embodiments are only some embodiments of this application, and not all embodiments. All other embodiments obtained by those skilled in the art based on the embodiments of this application without creative effort are within the scope of protection of this application.

[0089] To make the inventive objectives, technical solutions, and advantages of this application clearer, the embodiments of this application will be further described in detail below with reference to the accompanying drawings: In order to better understand the above-mentioned objectives, features, and advantages of this invention, the advantages of this invention will be further illustrated below by comparing the embodiments with the accompanying drawings and specific implementation methods.

[0090] This invention proposes a GPU-accelerated lattice Boltzmann simulation method for porous rock multiphase fluids, and the steps of this method are described in detail below:

[0091] Step S1: Use the D3019 velocity discretization model to identify porous rocks and obtain a digital core model (e.g., Figure 1 As shown), the model defines the velocity direction and weighting coefficients, defines two distribution functions using the LBM method, calculates macroscopic physical quantities using the relationship between the distribution functions and velocity moments, and describes the complex boundaries of the gas-liquid interface and porous rock using scalar tracer functions and step functions, respectively, thereby constructing a physical model (such as...). Figure 2 (as shown)

[0092] Furthermore, the specific process of step 1 is as follows:

[0093] 1) LBM method

[0094] LBM (Lattice Boltzmann Method) is a numerical method based on mesoscopic particle distribution functions. It describes the dynamic behavior of fluids by discretizing velocity directions and positions. Specifically, LBM treats the fluid as a discrete system composed of a large number of particles that move on a regular grid and interact through collision rules. In the D3Q19 velocity discretization model (a discretization method for simulating three-dimensional fluid dynamics), each cell contains 19 velocity directions, including one stationary particle, six particles along the principal axes, and twelve particles along the diagonals. Each velocity direction has a weighting coefficient used to calculate the particle's contribution.

[0095] LBM defines two independent distribution functions. f i and g i Describe the velocity distribution in the gas phase and liquid phase respectively, where the subscripts are... iRepresenting different discrete velocity directions, the D3Q19 velocity discretization model defines 19 velocity directions. e i and corresponding weighting coefficients w i :

[0096] The formula for calculating the velocity of a stationary particle in three-dimensional space is:

[0097]

[0098] Principal axis particles refer to particles that move along the three principal axes (x, y, z) of the coordinate system. The formula for calculating principal axis particles is:

[0099]

[0100] Diagonal particles refer to particles that move along the diagonal directions of the xy, xz, and yz planes. The formula for calculating diagonal particles is:

[0101]

[0102] The velocity moment is obtained by weighted summation of the distribution functions. Macroscopic physical quantities can then be calculated from the velocity moment. The formula for calculating macroscopic physical quantities in relation to the velocity moment and the distribution function is as follows:

[0103]

[0104] in, The density is the gas phase density. The density of the liquid phase is... u For macroscopic velocity fields, The density of the macroscopic velocity field;

[0105] 2) Describe the location of the gas-liquid interface and the complex boundaries of porous rocks.

[0106] Introducing scalar tracer functions ψ is used to describe the location of the gas-liquid interface. ψ>0 represents the gas phase region, ψ<0 represents the liquid phase region, and ψ=0 is the interface. The formula for calculating the relationship between ψ and the macroscopic density field is:

[0107]

[0108] in, For mixed density;

[0109] 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 internal pores and the solid portion of the rock) is described by a step function B(x): B=1 represents the solid region, B=0 represents the fluid region, and 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).

[0110] Step S2: Construct the equilibrium distribution function using macroscopic physical quantities in the physical model. Under the action of external forces, use the forcing scheme equation to calculate the LB equation describing fluid evolution. Construct a pseudo-potential multiphase model with the help of the Shan-Chen potential energy formula. Then use the semi-rebound scheme to process the solid-liquid boundary of porous rocks, thereby simulating the dynamic fluid behavior in porous rocks more accurately.

[0111] Furthermore, the specific process of step S2 is as follows:

[0112] 1) LB equation

[0113] Constructing the equilibrium distribution function using macroscopic physical quantities f i eq , g i eq The calculation formula is:

[0114]

[0115] In the formula f i eq and g i eq It is the equilibrium distribution function of the gas and liquid phases, where u is the fluid velocity;

[0116] Relaxation time t With fluid kinematic viscosity v The formula for calculating the relationship is:

[0117]

[0118] Among them, c s For the speed of sound in a grid;

[0119] To simulate the effects of external forces on the gas and liquid phases in fluid dynamics, external forces (such as gravity) are introduced into the LBM, using Guo's forcing scheme equations:

[0120]

[0121] Where Ftotal and Gtotal are the total external forces acting on the gas phase and liquid phase, respectively. F i and G i These are external forces, including gravity or other volume forces.

[0122] To simulate the evolution of fluids by modeling the motion and collisions of particles on a discrete lattice, the Boltzmann equation (LB) is used:

[0123]

[0124] in, t f and t g The corresponding relaxation time determines the viscous properties of the fluid; Δt For time step ,Δx The spatial step size;

[0125] 2) Constructing a pseudo-potential multiphase model

[0126] Equilibrium state that reaches a stable state without external forces ψ eq The calculation formula is:

[0127]

[0128] In the formula, A and B are constants that control the gas-liquid density ratio. Controlling the magnitude of surface tension, ψ eq This represents the equilibrium state of the interface.

[0129] Surface tension term The surface tension term is given by the Shan-Chen force. The surface tension effect at the phase interface can be described by the formula for calculating the Shan-Chen force:

[0130]

[0131] Among them, G controls the surface tension intensity. This is the surface tension term;

[0132] The evolution of the phase interface [the interface between two different phases (such as liquid and gas, solid and liquid, etc.)] is described by an improved pseudo-potential multiphase model, the calculation formula of which is:

[0133]

[0134] in, This is the interface relaxation time, with a value between 0.7 and 1.5.

[0135] 3) Describe the solid-liquid boundary of porous rocks

[0136] The solid-liquid boundary of porous rocks is treated using a half-bounce-back scheme; for nodes near the solid wall... x f its corresponding direction i The formula for the distribution function on is:

[0137]

[0138] in, j for i The opposite, that is e j =- e i , u w The wall velocity;

[0139] This format is simple, effective, and easy to implement in parallel.

[0140] Step S3, the numerical solution pseudocode process based on the lattice Boltzmann method, can efficiently simulate the dynamic behavior of fluids;

[0141] Furthermore, the specific process of step S3 is as follows:

[0142] Discretizing the above equations yields the explicit LBM evolution scheme. The numerical solution process for LBM can be divided into two main steps: Collision and Streaming. With the aid of a computer control system, the core structure of the entire algorithm is as follows:

[0143] for n = 0 : nstep do

[0144] for all lattice x do

[0145] Calculate the macroscopic quantities p_r^n(x), p_b^n(x), u^n(x). ^n(x)

[0146] Construct equilibrium states f_i^{eq,n}, g_i^{eq,n}, ^{eq,n}

[0147] Calculate the external force terms F_i^n, G_i^n, S_ ^n

[0148] Collision: Update f_i and g_i

[0149] renew

[0150] Migration: Pass the updated f_i and g_i to neighboring nodes.

[0151] end for

[0152] end for

[0153] The algorithm consists of two loops: the outer loop iterates through the time steps, and the inner loop iterates through all grid points. Each grid point sequentially performs macroscopic quantity calculations, equilibrium state construction, collision relaxation, and migration updates, thereby completing the evolution at each step. The solid-liquid boundary is handled using a semi-bounce scheme.

[0154] Based on the construction of physical models and the LBM method, the CUDA programming model is used to perform GPU parallel acceleration operations, thereby improving computational efficiency.

[0155] Furthermore, the GPU parallel acceleration process in the aforementioned steps is as follows:

[0156] The LBM algorithm is characterized by its explicit and local nature, making it very suitable for the single instruction multithreading (SIMT) parallel mode of GPUs. The hardware features of GPUs, such as multithreading and on-chip shared memory, can be fully utilized to highly parallelize the core steps of LBM and achieve high-performance computing.

[0157] In the specific implementation, the CUDA programming model is adopted; the LBM grid based on the D3Q19 model is divided into several CUDA blocks (basic execution units, a collection of CUDA threads), each block is responsible for the parallel computation of a subdomain; within the block, each CUDA thread (the smallest execution unit, a thread) is responsible for the evolution of a grid point; frequently accessed data is cached through shared memory (cached memory, located inside each block) to reduce time-consuming global memory access; when necessary, the synchronization instruction __syncthreads() is used to ensure data consistency.

[0158] In terms of multi-GPU parallelism, a domain decomposition strategy is adopted; the computation domain is divided into several sub-regions and assigned to different GPUs, which can make full use of the computing resources of multiple GPUs; coupling is achieved through boundary data exchange, which can ensure data consistency between different GPUs; and communication overhead is masked through a pipeline strategy, which overlaps transmission and computation, allowing the GPU to continue to perform computation while waiting for data transmission.

[0159] Based on the above algorithm, a complete parallel simulation system for gas-liquid two-phase flow was developed; the system adopts a modular design.

[0160] Furthermore, the main processes of the parallel simulation system for gas-liquid two-phase flow include:

[0161] 1) Preprocessing module: Import CT images or digital cores of porous rocks, reconstruct their internal pore structure; generate voxelized structured meshes, and rationally divide them into parallel partitions;

[0162] 2) Numerical solution module, which implements a CUDA-based multi-GPU parallel LBM solver; supports double-precision floating-point operations and allows setting parallel scale, partitioning strategy, etc.

[0163] 3) Boundary processing module, which realizes parallel processing of various boundary conditions such as wall, symmetry, periodicity, pressure, velocity, etc.; supports complex cases such as moving boundary and wetting contact angle;

[0164] 4) Post-processing module: reads instantaneous flow field in parallel and extracts key physical quantities; developed parallel visualization tools to render gas-liquid interface, velocity vector field, etc.

[0165] All the above modules work together through standard interfaces such as MPI communication to realize the entire process of gas-liquid two-phase flow simulation; users can input calculation parameters, monitor progress, and analyze results through a graphical interface to achieve human-computer interaction.

[0166] like Figure 3 and Figure 4 As shown, an embodiment of the system of the present invention is as follows:

[0167] 1. Initialization and parameter settings

[0168] This invention first imports the necessary Python libraries, such as numpy for scientific computing, taichi for GPU acceleration, and pyevtk for VTK format output.

[0169] Define the key parameters of the program:

[0170] The mesh sizes nx, ny, nz define the degree of discretization of the three-dimensional computational domain. Taking 232^3, a larger mesh means a more refined but more expensive simulation.

[0171] 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;

[0172] niu_l and niu_g are the kinematic viscosity 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.

[0173] psi_solid defines the wettability of a solid; 0.7 indicates hydrophilicity, which makes the liquid phase tend to adhere to the solid wall.

[0174] CapA defines the strength of surface tension, with 0.005 being a suitable dimensionless value;

[0175] Velocity and pressure are applied to the left and right boundaries in the x-direction, respectively. Boundary conditions are set as follows:

[0176] The left boundary has a fixed injection velocity vx_bcxl=0.001, corresponding to a fixed liquid phase flow rate;

[0177] The right boundary is fixed at a pressure of rho_bcxr=0.995 to ensure that the outlet is in the gas phase.

[0178] The remaining four boundaries are set to periodic conditions to simulate infinitely extending porous rock;

[0179] The wettability boundary psi_x_left=1.0 initially fills the left boundary 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.

[0180] 2. Field variable storage

[0181] LBMs involve a large number of field variables and require the use of limited GPU memory;

[0182] The distribution function f[i,j,k,l] is a four-dimensional array, where (i,j,k) is the spatial index and l is the discrete velocity index, taking values ​​from 0 to 18. Considering double-precision floating-point, there are 19 velocity directions, which is the most memory-intensive variable.

[0183] F is a buffer for the distribution function f, used to store intermediate states after collision and before propagation, and has the same type as f;

[0184] 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, resulting in lower overhead.

[0185] The densities rho_r and rho_b of the gas and liquid phases, and the region identifier psi, are three-dimensional arrays that occupy only one double space.

[0186] To save memory, the solid distribution solid[i,j,k] is marked as a boolean array, occupying only 1 byte;

[0187] Next, we define some basic variables for the LBM model:

[0188] e[l]=[(0,0,0),(±1,0,0)...] represents the discrete velocity directions, taking the D3Q19 model with 19 directions;

[0189] e_f represents the velocity direction in the floating-point version, used in macroscopic velocity calculations;

[0190] w[l]=[1 / 3, 1 / 18, ...] represents the isotropic Gaussian integral weights, ensuring isotropy;

[0191] LR[l] records the reverse mapping of 18 velocity directions for bounce-back format;

[0192] Finally, we define the matrix M[19,19] and its inverse inv_M[19,19] in the MRT collision operator to link macroscopic quantities with the moment space; the derivation of these matrices requires mathematical techniques such as Chapman-Enskog expansion, which are the theoretical foundation of LBM, and will not be elaborated here.

[0193] Finally, calculate the relaxation coefficient of the two-phase fluid:

[0194] wl and wg are calculated based on the kinematic viscosity niu_l and niu_g of the liquid and gas phases, respectively, to ensure accurate macroscopic fluid dynamic properties.

[0195] 3. Initialization

[0196] After the field variables are defined, they need to be assigned appropriate initial values, which can be done using the init() function:

[0197] Traverse the non-solid region (solid==0) and initialize its physical quantities:

[0198] The psi is initialized to -1, representing the gas phase, which matches the wettability of psi_solid=0.7;

[0199] The initial values ​​of density rho, velocity v, and other macroscopic quantities are determined based on hydrostatic relations.

[0200] The density distributions rho_r and rho_b of the two phases are assigned values ​​according to the sign of psi;

[0201] Finally, the initial value of the distribution function f is calculated using the corresponding equilibrium function feq to ensure that the initial field matches the macroscopic quantity.

[0202] The buffer F copies the value of f as the initial state;

[0203] The microstructure of porous rocks, such as the solid distribution solid, is loaded from an external file via init_geo(); in this example, the pre-generated geometry is read from a *.txt file in 0 / 1 form and assigned to solid; similarly, the initial distribution psi of the two phases can also be imported from an external file; this enables flexible initialization of complex porous structures.

[0204] 4. Main loop

[0205] After initialization, the program enters its main loop and repeatedly performs the classic LBM steps of "collision-propagation".

[0206] At each time step:

[0207] ① colission() performs a "collision", that is, relaxes the distribution function in momentum space:

[0208] Iterate through each fluid node (i,j,k) and complete the following steps:

[0209] Calculate the macroscopic density rho, velocity u, etc. at node (i,j,k) using f[i,j,k,l].

[0210] Calculate the equilibrium distribution function f_eq[l] from macroscopic quantities;

[0211] Introduce external force terms, such as gravity and buoyancy, and store them as F[l].

[0212] Using the MRT matrix M[l,m], transform f[l] and f_eq[l] to the moment space;

[0213] For non-conserved quantities in the moment space, perform BGK-type relaxation with relaxation coefficient w;

[0214] Transform back to velocity space and simultaneously add external force terms to obtain the distribution function after the collision;

[0215] Anisotropic relaxation ensures the physical properties of macroscopic fluids, such as incompressibility;

[0216] Store the distribution function after the collision into F[i,j,k,l] to complete the "collision";

[0217] ② streaming1() performs "propagation", that is, it translates the distribution function in physical space:

[0218] Iterate through all nodes (i,j,k) and shift by F according to e[l];

[0219] Specifically, in the degree direction l, F[i,j,k,l] propagates to its neighbor F[i+ex,j+ey,k+ez,l];

[0220] This step is performed synchronously throughout the entire flow field, simulating the free flight of fluid particles;

[0221] After propagation, the distribution function at the solid boundary is updated to reflect the no-slip condition;

[0222] ③ Boundary_condition() handles boundary conditions:

[0223] For periodic boundaries, the propagation of the distribution function f is continuous and requires no special treatment;

[0224] For velocity boundaries, such as the left boundary, f is reset with feq_bcxl; this balanced distribution corresponds to the injection velocity vx_bcxl.

[0225] For pressure boundaries, such as the right boundary, they are also reset with feq_bcxr, corresponding to the outlet density rho_bcxr;

[0226] For the solid-liquid boundary, a semi-rebound scheme is used to simulate the viscous no-slip condition, and the reverse velocity is found through LR[l].

[0227] Wetting properties are introduced by the contact angle dynamics equation and superimposed as an additional term to control the dynamic contact line;

[0228] ④ streaming3() copies the propagated distribution function F back to f and calculates macroscopic quantities:

[0229] First, the F of the entire flow field is copied to f to prepare for the next step of "collision";

[0230] Meanwhile, using the updated f and moment relations, the macroscopic density rho[i,j,k] and velocity u[i,j,k] at node (i,j,k) are calculated;

[0231] The interface scalar psi is updated by the density fields rho_r and rho_b;

[0232] At this point, a fluid iteration of one time step has been completed; repeating this process multiple times simulates the dynamic evolution of the flow.

[0233] ⑤ Output the physical quantities of the flow field to a file for post-processing analysis and visualization.

[0234] Every certain number of steps, such as 100 steps, I / O output is performed on the entire flow field;

[0235] Using the pyevtk library, scalar and vector fields can be organized into *.vtr format, which can be read by software such as Paraview;

[0236] The entire process is encapsulated in a high-level Python script, which allows for flexible parameter settings and control over the start and stop of the simulation. It is compiled and executed on the GPU via CUDA, giving full play to the advantages of parallel computing and enabling simulations of millions of grids to be completed in a few hours.

[0237] Although preferred embodiments of this application have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments as well as all changes and modifications falling within the scope of this application.

[0238] Obviously, those skilled in the art can make various modifications and variations to this application without departing from the spirit and scope of this application. Therefore, if such modifications and variations fall within the scope of the claims of this application and their equivalents, this application also intends to include such modifications and variations.

Claims

1. A GPU-accelerated porous rock multiphase fluid lattice Boltzmann simulation method, characterized in that, The application relates to a method for simulating the evolution of a fluid in a porous rock, comprising the following steps: Step S1, a digital core model is obtained by identifying the porous rock by using a D3Q19 velocity discrete model, velocity directions and weight coefficients are defined on the model, two distribution functions are defined by using the LBM method, macroscopic physical quantities are calculated by using the relationship between the distribution functions and the velocity moments, a scalar tracer function and a step function are used to respectively describe a gas-liquid interface and a complex boundary of the porous rock, and a physical model is constructed; Step S2, the macroscopic physical quantities in the physical model are used to construct equilibrium distribution functions, LB equation calculation of fluid evolution is carried out by using a forcing scheme equation under the action of an external force, a pseudo-potential multiphase model is constructed by using a Shan-Chen potential formula, and a semi-rebound format is used to process the solid-liquid boundary of the porous rock; 1) LB equation calculation Formulas for constructing equilibrium distribution functions from macroscopic physical quantities The formula is: ; In the formula It is the equilibrium distribution function of the gas phase. It is the equilibrium distribution function of the liquid phase. For fluid velocity; Nineteen velocity directions are defined for the D3Q19 velocity discrete model; To and The corresponding weighting coefficients; This refers to the gas phase density. The density of the liquid phase; Relaxation time The relationship between the fluid kinematic viscosity The calculation formula is: ； wherein, is the lattice sound speed; is the time step; is the spatial step; An external force is introduced in the LBM to simulate the influence of the external force effect on the gas phase and the liquid phase in fluid dynamics, and a forcing scheme equation of Guo is used: ; wherein, is the total external force on the gas phase; is the total external force on the liquid phase, and is an external force term, including gravity or other body forces; is the relaxation time of the gas phase; is the relaxation time of the liquid phase; The motion and collision of particles on a discrete lattice are simulated by using the LB equation to simulate the evolution of the fluid, and the LB equation is: ; wherein, is the relaxation time of the gas phase; is the relaxation time of the liquid phase, determining the viscous properties of the fluid; is the time step; is the spatial step; 2) Construction of a pseudo-potential multiphase model Equilibrium state in which the system reaches a stable state without external force The calculation formula is: ; Where A, B are constants, control the gas-liquid density ratio, Control the surface tension size, Is the equilibrium state of the interface; Is a scalar tracer function; Surface tension term describing the effect of surface tension at phase interfaces Given by the Shan-Chen potential, the formula for calculating the Shan-Chen potential is: ; where G controls the surface tension strength, is the surface tension term; is the lattice sound speed; The evolution of the phase interface is described by an improved pseudo-potential multiphase model, and the calculation formula of the pseudo-potential multiphase model is: ; wherein is the interface relaxation time, which is in the range 0.7-1.5; 3) Description of the solid-liquid boundary of the porous rock The solid-liquid boundary of porous rock is treated by a half-reflecting format; for the nodes close to the solid wall , the distribution function formula in the corresponding direction is ； wherein is the inverse of , , is the wall velocity; The format is simple and effective, and is easy to realize in parallel; Step S3, a numerical solution pseudo-code process based on the lattice Boltzmann method.

2. The GPU-accelerated multi-phase fluid lattice Boltzmann simulation method for porous rocks according to claim 1, wherein, The specific process of the step S1 is as follows: 1) LBM method LBM defines two independent distribution functions and respectively depict the velocity distribution of gas and liquid phase, where the subscript denotes different discrete velocity directions, and the D3Q19 velocity discrete model defines 19 velocity directions and corresponding weight coefficients : The calculation formula of a stationary particle with zero velocity in a discrete velocity space is: ; The calculation formula of a particle in the principal axis direction is: ； ; The calculation formula of a particle in the diagonal direction is: ; ； The velocity moment is obtained by weighted summation of the distribution functions, and the macroscopic physical quantity can be calculated through the velocity moment, and the calculation formula of the macroscopic physical quantity related to the distribution function is: ； wherein, is the gas phase density, is the liquid phase density, is the density of the macroscopic velocity field; 2) Description of the position of the gas-liquid interface and the boundary of the porous rock The position of the gas-liquid interface is given by the scalar tracer function The description, wherein, represents the gas phase region, represents the liquid phase region, at which the interface is located; The calculation formula related to the macroscopic density field is: ; wherein, is the mixed density; The boundary of the porous rock is described by using a step function B(x), which separates the solid region and the fluid region and is applied to various porous rock models: B=1 represents the solid region, B=0 represents the fluid region, and wettability is introduced by a dynamic contact angle model.

3. The GPU-accelerated multi-phase fluid lattice Boltzmann simulation method for porous rocks according to claim 1, wherein, The physical model and the LBM method are accelerated by using GPU parallel operation, and the specific steps are as follows: A CUDA programming model is used: the LBM grid based on the D3Q19 model is divided into a plurality of CUDA blocks, each block is responsible for parallel calculation of a subdomain; in the block, each CUDA thread is responsible for the evolution of a lattice point; frequently accessed data is cached through a shared memory, reducing time-consuming global memory access; data consistency is ensured by using a synchronization instruction __syncthreads (). In terms of multi-GPU parallelism, a domain decomposition strategy is used: the calculation domain is divided into a plurality of sub-regions, which are respectively assigned to different GPUs, and the calculation resources of multiple GPUs are used; the coupling is realized through boundary data exchange to ensure the data consistency between different GPUs; the communication overhead is covered by using a pipeline strategy, and the transmission and calculation are overlapped, so that the GPU continues to calculate while waiting for data transmission.

4. A gas-liquid two-phase flow parallel simulation system applied to the GPU-accelerated multi-phase flow lattice Boltzmann simulation method for porous rocks according to claim 1, characterized in that, The process of modular design includes: 1) Pre-processing module, import CT images or digital cores of porous rock, reconstruct its internal pore structure; generate voxelized structured grid, reasonable division of parallel partition; 2) Numerical solution module, realize multi-GPU parallel LBM solver based on CUDA; Support double-precision floating-point operation, set parallel scale, partition strategy; 3) Boundary processing module, realize parallel processing of wall, symmetry, periodic, pressure, velocity boundary conditions; support dynamic boundary, wet contact angle complex situation; 4) Post-processing module, parallel read transient flow field, extract key physical quantities; develop parallel visualization tools, render gas-liquid interface, velocity vector field; All the above modules work together through MPI communication to realize the whole process of gas-liquid two-phase flow simulation; users input calculation parameters, monitor progress, analyze results through graphical interface to realize human-computer interaction.

Images