Three-phase foam flow injection simulation method for enhanced oil displacement of high-water-content oil reservoir

Abstract

The invention provides a three-phase foam flow injection simulation method for enhanced oil displacement of a high-water-content oil reservoir, and belongs to the field of computational fluid mechanics. An immiscible three-phase flow calculation model is constructed under a lattice Boltzmann framework, and a continuous surface force model and a re-coloring algorithm are introduced into the immiscible three-phase flow calculation model; repulsive force is added between gas-phase foams meeting requirements, so that the three-phase flow model has the capability of simulating non-coalescence foam behaviors; wetting boundary conditions are applied according to a diffusion interface idea, and solid wall surface boundaries with specific wettability under complex geometric conditions are processed; inputting three-phase fluid physical property parameters, providing auxiliary parameters necessary for executing simulation to operate the immiscible three-phase flow calculation model, and obtaining multiple groups of data representing the three-phase flow state in the foam injection oil displacement process through time-marching simulation. According to the invention, the blank of the lattice Boltzmann method in the aspect of simulating the three-phase foam flow of oil, water and gas is filled, and high-precision simulation can be carried out on the three-phase foam injection oil displacement process.

Description

Technical Field

[0001] The invention belongs to the technical field of computational fluid dynamics, and in particular relates to a three-phase foam flow injection simulation method for enhanced oil recovery in high water content oil reservoirs. Background Art

[0002] At present, my country's exploited oil fields have generally entered a high water content stage, the development potential of oil fields has gradually declined, and the traditional water drive oil recovery method has poor recovery effect. In order to maintain crude oil production and extend the development life of old oil fields, it is urgent to innovate the recovery enhancement technology. As a fluid with a high apparent viscosity, foam has the characteristics of being stable when it encounters water, so it can maintain stable properties in high-permeability areas with more water. In addition, foam usually exhibits a high seepage resistance, and the seepage resistance increases with the increase of permeability, so it can effectively block the high permeability layer, improve the water injection profile, and thus improve the water injection sweep efficiency. In addition, the foaming agent is generally a surfactant with good performance, which can reduce the oil-water interfacial tension and play the role of emulsification and carrying oil displacement. Therefore, the enhanced recovery effect of injecting gas in the form of foam is better.

[0003] Three-phase seepage experiment is an important means to study the seepage law of foam flooding, but the cost of experimental research is high. At the same time, it is difficult to accurately obtain important flow field information such as foam morphology and velocity and pressure under three-dimensional conditions due to the limitation of observation methods, resulting in the current understanding of the three-phase seepage mechanism of foam flooding is still incomplete. In order to make up for the shortcomings of experimental research and deepen the understanding of the migration law of foam fluid in high water-content reservoirs in an economical, efficient and intuitive way, the use of high-precision simulation methods to simulate the three-phase seepage process of oil, water and gas in porous media has become an effective alternative. It can truly reproduce the enhanced oil recovery process of foam fluid plugging and profile adjustment, and optimize the optimal operating parameters for enhanced recovery based on the simulation results.

[0004] In the field of numerical simulation, numerical simulation for two-phase flow started early, and relevant literature has widely covered multiple sub-fields such as flow state analysis, heat and mass transfer, pressure drop characteristics and phase distribution, forming a relatively systematic research system. In contrast, the study of three-phase flow faces more complex interactions between fluids and changes in interface topological structure due to the existence of additional phases, making the modeling and simulation of interface dynamics extremely challenging. Especially in immiscible three-phase flow, the singularity of the interface curvature in the coexistence region of the three-phase fluid and the topological changes of different fluid interfaces are the main difficulties in numerical modeling. At present, interface modeling methods are generally divided into sharp interfaces and diffuse interfaces. Among them, the sharp interface method can clearly present the interface structure, while the diffuse interface method is better at capturing the topological changes of the interface. Its implementation strategy is to allow different fluid components to mix at the interface to form a thin diffusion layer, achieve a smooth transition of fluid physical parameters, and introduce a marker function to distinguish different fluids. The main advantage of the diffuse interface is that it can avoid the problem of the inability to model the interfacial tension caused by the singularity of the curvature in the coexistence region of the three-phase fluid, and at the same time has the ability to capture complex interface topological changes, such as the rupture and coalescence of foams.

[0005] When simulating multiphase flow, the way boundary conditions are handled can significantly affect the simulation results. Traditional multiphase flow simulation methods usually strictly apply no-slip boundary conditions at the wall, which can easily lead to the non-physical phenomenon that the three-phase contact line cannot move, or the occurrence of stress singularities near the solid wall, which leads to the divergence of numerical simulation.

[0006] The lattice Boltzmann color model is a numerical model developed by introducing a color indicator function and its gradient under the lattice Boltzmann framework. The color indicator function is used to characterize different phases, while its gradient is intended to represent the interfacial tension phenomenon in multiphase flow. The recoloring operator plays a key role in the color model. It can promote the separation of different fluids along the interface normal at the mesoscopic level, thus ensuring the conservation of mass of different fluid phases during the simulation process, obtaining a sharper interface and taking into account the accuracy and stability of the simulation process. When implementing multiphase flow simulation, the lattice Boltzmann model does not require complex interface tracking, capturing or reconstruction technology, and can automatically describe the dynamic behaviors of interface deformation, rupture and merging. In addition, the color model inherits many advantages of the lattice Boltzmann method, including natural parallel characteristics, easy code writing, and easy processing of complex geometric boundaries. Therefore, it has significant advantages in multiphase seepage simulation at the pore scale and the study of complex interface dynamic behaviors.

[0007] Researchers have developed a variety of three-phase flow solvers under the lattice Boltzmann framework. These solvers focus on interfacial tension modeling and improving the model's ability to handle a wide range of fluid density ratios and viscosity ratios. However, they do not have the ability to simulate "non-coalescing" foam flow, and therefore do not have the ability to simulate three-phase seepage during foam injection enhanced oil recovery. Summary of the invention

[0008] In order to solve the above problems, the present invention provides a three-phase foam flow injection simulation method for enhanced oil recovery in high water-cut oil reservoirs. In particular, the present invention is suitable for foam-enhanced crude oil recovery processes in high water-cut oil reservoirs and multi-scale and multi-physical field simulations of immiscible foam flows of oil, water and gas.

[0009] In order to achieve the above object, the present invention provides the following technical solutions:

[0010] A three-phase foam flow injection simulation method for enhanced oil recovery in high water content reservoirs, comprising:

[0011] An immiscible three-phase flow computational model is constructed under the lattice Boltzmann framework, and a continuous surface force model and a recoloring algorithm are introduced into the immiscible three-phase flow computational model; wherein the immiscible three-phase flow is composed of oil, gas and water; the continuous surface force model is used to generate interfacial tension between different phase interfaces, so that the model has the ability to capture the interface changes of droplets and bubbles in porous media; the recoloring algorithm is used to separate the fluids of each phase along the interface normal and maintain the conservation of mass of each phase to generate an identifiable interface;

[0012] When executing the collision step in the moment space, an additional source term is introduced on the right side of the distribution function equation, and the inter-foam repulsion force is introduced into the immiscible three-phase flow calculation model in the form of volume force to correct the immiscible three-phase flow calculation model, wherein the inter-foam repulsion force is calculated by the inter-foam repulsion force model to curb the non-physical aggregation of foams; after executing the collision step, the distribution function in the moment space is first restored by inverse transformation to obtain the total density distribution function, and then the density distribution function is assigned to the three different fluids using the re-coloring algorithm to ensure that different fluids are separated along the interface normal and the mass conservation of each phase of the fluid is maintained; the contact angle geometric formula is used to apply the wetting boundary condition at the fluid-solid boundary to reflect the different wetting characteristics of different fluids on the solid wall;

[0013] The required three-phase fluid physical property parameters are input and the auxiliary parameters necessary for simulation are provided to run the immiscible three-phase flow calculation model. Multiple sets of data characterizing the three-phase flow state during foam injection flooding are obtained through time advancement simulation.

[0014] Preferably, the continuous surface force model is used to generate interfacial tension F between different phase interfaces. s , and its calculation formula is:

[0015]

[0016] in:

[0017] G kl =φ l ▽φ k -φk ▽φ l ;

[0018]

[0019] Where: I is the unit matrix, k, l oil, water, gas represent any two phases in the three-phase fluid; σ kl is the interfacial tension coefficient between fluid k and fluid l; ρ is the density of the three-phase mixed fluid, ρ k is the local density of fluid k, ρ l is the local density of fluid l; φ l is the volume fraction of fluid l, is the density of pure fluid l, φ k is the volume fraction of fluid k, is the density of pure fluid k, n kl is the unit external normal vector of the kl interface, G kl is the color gradient, C kl is the parameter that controls the interfacial tension strength of the kl interface.

[0020] Preferably, the repulsive force between foams is calculated by a repulsive force model between foams. For the repulsive force between foams that is generated between two activated repulsive force points P and Q and that moves away from each other, the specific calculation is as follows:

[0021] The repulsive force generated by point Q on point P is:

[0022]

[0023] Where: K is a constant that controls the strength of the repulsive force; l PQ is the unit vector from P to Q, φ g represents the gas phase volume fraction, x P is the position coordinate of point P, x Q is the position coordinate of point Q, h PQ is the distance between points P and Q;

[0024] The total repulsive force F on point P rep It is expressed as:

[0025]

[0026] Where: n P is the external normal vector of point P;

[0027] The total repulsive force F rep By taking the repulsive force between bubbles into account in the resultant force F, the repulsive force between bubbles can be introduced into the three-phase seepage model.

[0028] Preferably, the re-coloring operation of the re-coloring algorithm is performed on the density distribution function in the velocity space, so before the re-coloring operation, the post-collision velocity moment m′ is first inversely transformed to obtain the post-collision total density distribution function f′, that is, f′(x, t)=M -1 m′(x,t), where M is the orthogonal transformation matrix; then, the total density distribution function after the collision is re-colored:

[0029]

[0030] Where: x represents the spatial position, t is the time, and f i,k ″ represents the distribution function of fluid k after re-coloring, e i is the grid velocity of the D2Q9 model, w i for e i Directional weight coefficient; ρ k is the local density of fluid k, ρ l is the local density of fluid l, ρ is the total density of the three-phase mixed fluid, n kl is the unit external normal vector of the kl interface, β kl is a separation parameter related to the fluid kl interface thickness and satisfies β kl =β lk .

[0031] Preferably, the wetting boundary condition is applied at the fluid-solid boundary by using the contact angle geometric formula, specifically including:

[0032] In any local part of the solid wall, the computational domain is divided into two parts, the fluid domain and the solid domain, by the solid wall boundary. The grid points in the solid domain are recorded as solid grid points, and the remaining grid points are fluid grid points. The solid boundary grid points are defined as follows: for a solid grid point, if any grid point among the 8 grid points around it is a fluid grid point, it is regarded as a solid boundary grid point. Before matching the contact angle, the external normal vector of the solid boundary is calculated first. Assuming that a solid boundary grid point is S, its external normal vector n S The calculation formula is:

[0033]

[0034] Where: c j is the mesoscopic velocity in the eighth-order isotropic discrete space, x S is the position coordinate of point S; s(x) is the indicator function, which takes the value 1 if and only if x is a solid grid point, otherwise it takes the value 0; ω(c j 2 ) is the eighth-order weight coefficient, and its value is as follows:

[0035]

[0036] At any fluid grid point, the sum of the volume fractions of the three fluids is 1, that is, φ o +φ w +φ g =1; therefore, only the weighted contact angles θ of the oil phase and the water phase need to be calculated o and θ w , and its calculation formula is:

[0037] (1-φ o )θ o =φ w θ ow +φ g θ og ;

[0038] (1-φ w )θ w =φ o θ wo +φ g θ wg ;

[0039] Where: k, l represent any two phases of the three-phase fluid o, w, g, φ o is the volume fraction of oil phase, φ w is the volume fraction of water phase, φ g is the gas phase volume fraction; θ kl is the contact angle between the interface of fluids k and l and the solid wall measured from the side of fluid k; contact angle θ kl and θ lk They are complementary angles, i.e. θ lk =π-θ kl ;

[0040] After obtaining the normal vector n outside point S S and the weighted contact angle θ k After that, two lines can be drawn with n S The angles formed are all |π / 2-θ k |Characteristic line l C1 , l C2 , both intersect with the grid lines respectively; the intersection point between the characteristic line and the grid line in the solid domain is an invalid intersection point; the first intersection points of the two characteristic lines in the fluid domain are S1 and S2 respectively, and it is necessary to ensure that the adjacent grid points of the intersection points S1 or S2 on the horizontal or vertical grid lines are all fluid grid points; if the adjacent grid points are not fluid grid points, it is necessary to continue to extend the characteristic line until the above conditions are met;

[0041] Assuming that A and B are the fluid grid points closest to the intersection S1 in the horizontal direction, and C and D are the fluid grid points closest to the intersection S2 in the vertical direction, the fluid volume fractions at points A, B, C, and D are known through flow field calculation, so the fluid volume fractions at the intersections S1 and S2 can be obtained by linear interpolation, that is, the volume fractions of fluid k at the intersections S1 and S2 are:

[0042]

[0043] Where: represents the position coordinates of the intersection point S1, represents the position coordinates of the intersection point S2, x A Represents the position coordinates of point A, x B Indicates the position coordinates of point B, x C Indicates the position coordinates of point C, x D Represents the position coordinates of point D; represents the abscissa of the intersection point S1; represents the ordinate of the intersection point S2; x A represents the horizontal coordinate of point A, x B Indicates the horizontal coordinate of point B; y C Indicates the ordinate of point C, y D represents the ordinate of point D;

[0044] get and After that, the wall wettability is determined. or Assign values ​​to the solid boundary grid points S:

[0045]

[0046] Where: θ k represents the weighted contact angle of fluid k;

[0047] After traversing all boundary solid grid points and executing the above calculation strategy, the volume fractions of oil, water, and gas three-phase fluids matching the required contact angles on all boundary solid grid points can be obtained, thereby affecting the color gradient and interface normal vector on the boundary fluid grid points; in this way, the wetting boundary conditions are implicitly applied.

[0048] The present invention also provides a three-phase foam flow injection simulation system for enhanced oil recovery in high water content oil reservoirs, comprising:

[0049] A model building module is used to build an immiscible three-phase flow calculation model under the lattice Boltzmann framework, and introduce a continuous surface force model and a recoloring algorithm into the immiscible three-phase flow calculation model; wherein the immiscible three-phase flow is composed of oil, gas and water; the continuous surface force model is used to generate interfacial tension between different phase interfaces, so that the model has the ability to capture the interface changes of droplets and bubbles in porous media; the recoloring algorithm is used to separate the phases of fluid along the interface normal and maintain the mass conservation of each phase to generate an identifiable interface;

[0050] The model correction module is used to introduce an additional source term on the right side of the distribution function equation when executing the collision step in the moment space, and introduce the inter-foam repulsion force into the immiscible three-phase flow calculation model in the form of volume force to correct the immiscible three-phase flow calculation model, wherein the inter-foam repulsion force is obtained by calculating the inter-foam repulsion force model, and is used to curb the non-physical aggregation of foams; after executing the collision step, the distribution function in the moment space is first restored by inverse transformation to obtain the total density distribution function, and then the density distribution function is assigned to the three different fluids using the re-coloring algorithm to ensure that different fluids are separated along the interface normal and the mass conservation of each phase fluid is maintained; the contact angle geometric formula is used to apply the wetting boundary condition at the fluid-solid boundary to reflect the different wetting characteristics of different fluids on the solid wall;

[0051] The simulation module is used to input the required three-phase fluid physical properties and provide the auxiliary parameters necessary for executing the simulation to run the immiscible three-phase flow calculation model, and obtain multiple sets of data characterizing the three-phase flow state during the foam injection flooding process through time advancement simulation.

[0052] The present invention also provides a computer device, comprising a memory, a processor and a computer program stored in the memory, wherein the processor executes the computer program to implement any step in the three-phase foam flow injection simulation method for enhanced oil recovery in high water-content reservoirs.

[0053] The present invention also provides a computer-readable storage medium having a computer program stored thereon. When the computer program is loaded by a processor, the computer program can execute any of the steps in the three-phase foam flow injection simulation method for enhanced oil recovery in high water-content reservoirs.

[0054] The three-phase foam flow injection simulation method for enhanced oil recovery in high water content reservoirs provided by the present invention has the following beneficial effects:

[0055] The present invention introduces a continuous surface force model and a recoloring algorithm into the immiscible three-phase flow calculation model, so that the model itself has the ability to simulate three-phase flow under the full range of interfacial tension conditions and maintains the advantages of mass conservation of fluid components; the introduction of the repulsive force model can curb the non-physical aggregation behavior of foam in numerical simulation and correctly reproduce the behavior of foam in porous media; the contact angle geometry formula is used under the diffusion interface framework to apply wetting boundary conditions, avoid stress singularity, flexibly adjust the contact angle, and correctly reflect the wettability of the fluid. Based on the above invention, the constructed immiscible three-phase flow calculation model can accurately simulate the injection of three-phase foam flow. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] In order to more clearly illustrate the embodiment of the present invention and its design scheme, the following briefly introduces the drawings required for this embodiment. The drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0057] Figure 1 This is a flow chart of a three-phase foam flow injection simulation method for enhanced oil recovery in a high water content reservoir according to Example 1 of the present invention;

[0058] Figure 2 Schematic diagram of the range of repulsive force between two nearly contacting bubbles;

[0059] Figure 3 Schematic diagram of the implementation of wetting boundary conditions for complex solid wall surfaces;

[0060] Figure 4 The porous medium geometry and fluid distribution before water flooding; in the figure, black is solid, red is oil phase, and green is water phase;

[0061] Figure 5 This is the curve of residual oil saturation changing with the PV number of injected fluid during water flooding;

[0062] Figure 6 This is the oil-water two-phase distribution diagram in the porous medium after water flooding is stabilized;

[0063] Figure 7 After the foam is injected, the foam first gathers in the throat of the wide channel, showing the blocking effect of blocking the dominant channel and prompting the water phase to open up a new flow channel;

[0064] Figure 8 This is the distribution diagram of the foam being blocked and retained in the narrow channel after the foam is injected;

[0065] Fig. 9 When the plugging and profile control effect of the foam appears, four new displacement fronts are generated;

[0066] Fig.10 Pure water displacement occurs at the top of the porous medium;

[0067] Fig.11 Schematic diagram of the complete destruction of the lamellar residual oil structure at the top of the porous medium. DETAILED DESCRIPTION

[0068] In order to enable those skilled in the art to better understand the technical solution of the present invention and implement it, the present invention is described in detail below in conjunction with the accompanying drawings and specific embodiments. The following embodiments are only used to more clearly illustrate the technical solution of the present invention, and cannot be used to limit the scope of protection of the present invention.

[0069] Example 1

[0070] In the study of three-phase flow, the presence of additional phases will lead to more complex interactions between fluids and changes in the interface topology, which in turn makes the modeling and simulation of the interface dynamics process extremely complicated. Especially in immiscible three-phase flow, the singularity of the interface curvature in the coexistence region of the three-phase fluids and the topological changes of the interfaces of different fluids are the main difficulties in numerical modeling. In order to simulate the enhanced crude oil recovery process of foam injection into high water-content reservoirs in complex porous media, the present invention aims to modify the lattice Boltzmann immiscible three-phase flow calculation model by introducing repulsive forces between foams, and apply wetting boundary conditions using the contact angle geometry formula under the framework of the diffusion interface.

[0071] However, in numerical simulations, foams usually show a tendency to aggregate, which is inconsistent with the separated foam under the combined influence of van der Waals attraction and the repulsive force generated by the surfactant adsorption layer under actual working conditions. Therefore, it is necessary to introduce repulsive force between foams to maintain their non-aggregated state. Based on this, the present invention proposes a simulation algorithm suitable for immiscible three-phase foam flow at the pore scale to fill the research gap in this field and lay a model foundation for exploring and understanding the migration mechanism of foam flow in porous media.

[0072] Specifically, the present invention provides a three-phase foam flow injection simulation method for enhanced oil recovery in high water content reservoirs, specifically a simulation method for pore-scale immiscible oil, water, and gas three-phase foam flow under a lattice Boltzmann framework, such as Figure 1 As shown, the method includes:

[0073] An immiscible three-phase flow computational model is constructed under the lattice Boltzmann framework, and a continuous surface force model and a recoloring algorithm are introduced into the immiscible three-phase flow computational model; wherein the immiscible three-phase flow consists of oil, gas and water; the continuous surface force model is used to generate interfacial tension between different phase interfaces, so that the model has the ability to capture the interface changes of droplets and bubbles in porous media; the recoloring algorithm plays a key role in the color model, which is used to separate the phases of fluid along the interface normal, maintain the mass conservation of each phase, and ensure that a diffusion interface of reasonable thickness is generated under the full range of interfacial tension conditions.

[0074] When executing the collision step in the moment space, an additional source term is introduced on the right side of the distribution function equation, and the inter-foam repulsion force is introduced into the immiscible three-phase flow calculation model in the form of volume force to correct the immiscible three-phase flow calculation model, wherein the inter-foam repulsion force is obtained by calculating the inter-foam repulsion force model, which is used to curb the non-physical aggregation of foams; after executing the collision step, the distribution function in the moment space is first restored through an inverse transformation to obtain the total density distribution function, also called the total density distribution function in the velocity space; and then the re-coloring algorithm is used to assign density distribution functions to different fluids to ensure that different phases are separated along the interface normal and maintain mass conservation.

[0075] In numerical simulations, foams usually show a tendency to aggregate, which is inconsistent with the separated state of foams under actual working conditions under the combined influence of van der Waals attraction and the repulsive force generated by the surfactant adsorption layer. Therefore, it is necessary to introduce a repulsive force model between foams to maintain their non-aggregate state and correctly reproduce the foam behavior.

[0076] In the framework of diffusion interface, the contact angle geometric formula is used to impose wetting boundary conditions at the fluid-solid boundary to reflect the different wetting characteristics of different fluids on the solid wall.

[0077] The main advantages of the diffusion interface are: avoiding the interface tension that cannot be modeled due to the singular curvature of the phase interface in the coexistence zone of three-phase fluids; and having the ability to capture complex interface topological changes, such as droplet breakup and coalescence. The diffusion interface method allows the contact line to move by introducing a diffusion mechanism, so no-slip boundary conditions can still be used when dealing with wetting boundary conditions. This not only greatly simplifies the implementation of wetting boundary conditions under irregular wall conditions, but also avoids the introduction of slip models and slip lengths, which are dependent on the microscopic properties of materials and fluids and are difficult to measure.

[0078] Input the required three-phase fluid physical parameters and provide the auxiliary parameters necessary for simulation to run the immiscible three-phase flow calculation model, and obtain multiple sets of data characterizing the three-phase flow state during foam injection and oil recovery through time advancement simulation. Among them, the physical parameters include fluid viscosity, interfacial tension coefficient of the two fluids, three-phase contact angle, etc.; auxiliary parameters include D2Q9 velocity set, transformation matrix, diagonal relaxation matrix, etc.

[0079] Furthermore, the construction of the immiscible three-phase flow calculation model includes the following contents:

[0080] First, three different immiscible fluids are defined in the model. In the foam flooding system, the three-phase fluid consists of oil, water, and gas, and their density distribution functions are f o,i 、f w,i 、f g,i , where subscripts o, w, g represent the corresponding fluids (same as below), and i is the lattice velocity subscript. In the embodiment of the present invention, a two-dimensional nine-velocity (D2Q9) lattice model is used, and the value range of i is 0-8.

[0081] Total density distribution function f i It is defined as the sum of the density distribution functions of the three-phase fluid, that is:

[0082] f i =f o,i +f w,i +f g,i (1)

[0083] In the process of performing collision operation on the total density distribution function, in order to improve the numerical stability of the simulation, a multi-relaxation collision model is proposed, so the collision process is carried out in the moment space. Let f = (f0, f1, ..., f8) T , the corresponding velocity moment is m=Mf, where M is the orthogonal transformation matrix, and the explicit expression is as follows:

[0084]

[0085] Therefore, the collision operation in moment space can be expressed as:

[0086]

[0087] Where: m′(x,t) is the velocity moment at spatial position x and time t after the collision, m eq is the equilibrium state function of the moment space, and D is the relaxation matrix; is the force term in moment space, I is the unit matrix, δ t is the time step.

[0088] The relaxation matrix D is a non-negative diagonal matrix with diagonal elements d i (i = 0 ~ 8) are 9 relaxation parameters, each of which corresponds to a moment component, and the relaxation time corresponding to various moment components can be adjusted independently, thereby improving the numerical stability of the lattice Boltzmann simulation under high viscosity ratio conditions. The relaxation parameters of D are:

[0089] D=diag[1,d1,d2,1,d4,1,d6,d7,d8] (4)

[0090] Where: d7 and d8 are terms related to the dynamic viscosity μ of the three-phase mixed fluid, and d7 = d8; for symmetry considerations, d4 = d6 needs to be satisfied. In this model, in order to obtain the viscosity-independent wall position while taking into account numerical stability, the relaxation parameters are selected as follows:

[0091]

[0092] In formula (5): ρ = ρ o +ρ w +ρ g is the density of the three-phase mixed fluid. The dynamic viscosity μ can be obtained by harmonic averaging:

[0093]

[0094] The equilibrium function m in moment space eq By transforming the equilibrium distribution function in velocity space Obtained, its expression is as follows:

[0095]

[0096] Where: u is the velocity of the fluid, u x and u y are the components of u in the horizontal and vertical directions respectively.

[0097] Similarly, the force term in moment space It can also be obtained through similar transformations, as shown below:

[0098]

[0099] Where: F x 、F y are the horizontal and vertical components of force F, respectively. It is particularly noted that in the embodiment of the present invention, force F is the interfacial tension F s and the repulsive force between bubbles F rep For oil, water, and gas three-phase fluids, the interfacial tension F s The calculation formula is:

[0100]

[0101] in:

[0102] G kl =φ l ▽φ k -φ k ▽φ l (11)

[0103]

[0104] In the above formula: k and l can traverse o, w, and g respectively; I is the unit matrix, k and l represent any two phases in the three-phase fluid; σ kl is the interfacial tension coefficient between fluid k and fluid l; ρ is the density of the three-phase mixed fluid, ρ k is the local density of fluid k, ρ l is the local density of fluid l; φ l is the volume fraction of fluid l, is the density of pure fluid l, φ k is the volume fraction of fluid k, is the density of pure fluid k, n kl is the unit external normal vector of the kl interface, G kl is the color gradient, C kl is the parameter that controls the interfacial tension strength of the kl interface.

[0105] Although F s It can produce interfacial tension effect, but it cannot ensure phase separation, so it is necessary to introduce a recoloring operator to ensure that oil, water and gas are immiscible. The recoloring operation is performed on the density distribution function in the velocity space, so before the recoloring operation, m′ is first transformed by inverse transformation to obtain the total density distribution function f′ after the collision, that is, f′(x, t) = M -1 m′(x,t), where M is the orthogonal transformation matrix. Then, the total density distribution function after the collision is re-colored:

[0106]

[0107] Where: x represents the spatial position, t is the time, and f i,k ″ represents the distribution function of fluid k after re-coloring, e i is the grid velocity of the D2Q9 model, w i for e i Directional weight coefficient; ρ k is the local density of fluid k, ρ l is the local density of fluid l, ρ is the total density of the three-phase mixed fluid, n kl is the unit external normal vector of the kl interface. kl is a separation parameter related to the fluid kl interface thickness and satisfies β kl =β lk , and its calculation formula is:

[0108]

[0109] Where: g(x) is a piecewise function about x, which is defined as:

[0110]

[0111] Finally, the re-colored oil, water, and gas distribution functions are migrated to adjacent grid points, namely:

[0112] f i,k (x+e i δ t ,t+δ t )=f i,k ″(x,t) (18)

[0113] After the migration is completed, the new distribution function is used to calculate the fluid macroscopic physical quantities such as density and velocity, as shown in formula (19), and enter the next cycle.

[0114]

[0115] In the three-phase foam flow demonstration example based on the present invention, the introduction of the repulsive force between foams includes the following contents:

[0116] As mentioned above, the introduction of the repulsive force term is intended to curb the non-physical aggregation between bubbles. In a three-phase foam flow system, a surfactant is introduced as a foaming agent and adsorbed on the foam surface. When two independent bubbles approach each other, complex microscopic and mesoscopic interactions will be generated locally in the contact area, such as van der Waals forces, electrostatic forces, hydrophobic / hydrophilic interactions, etc. Under the combined influence of close contact effects such as intermolecular van der Waals forces and electrostatic repulsion generated by foam surfactant molecules, the foam can maintain its topological structure for a long time without aggregation. For this reason, the present invention introduces repulsive forces between bubbles to simulate this phenomenon, and the effectiveness of the repulsive force term needs to meet certain conditions.

[0117] like Figure 2 As shown in the figure, the dotted circle relatively located in the upper left corner and the dotted circle relatively located in the lower right corner represent the influence range of points P and Q respectively, and the radius of the circle is the critical value for generating repulsive force. In addition to point P, point Q can also affect points 1, 2, and 3; similarly, in addition to point Q, point P can also affect points 4, 5, 6, and 7. The black area in the upper left corner and the gray area in the lower right corner represent two independent bubbles that are in close contact with each other. The two activated repulsive force points P and Q will generate repulsive forces pointing to each other. To activate the repulsive force, three conditions must be met: ①P and Q are both located in the diffusion interface area of ​​the bubble to which they belong; ②The distance between P and Q is less than a certain critical value; ③P and Q belong to grid points in different bubbles. If the above three conditions cannot be met at the same time, the repulsive force cannot be activated and the default value 0 is taken.

[0118] To satisfy condition ①, it is necessary to examine the gas volume fraction φ at points P and Q respectively. g , if 0<φ g (x P )<1、0<φg (x Q )<1, it is considered that both points P and Q are located within the diffusion interface.

[0119] To meet condition ②, it is necessary to manually set the maximum threshold distance L that can activate the repulsive force. max , and space it from P and Q For comparison, if h PQ ≤L max , a repulsive force can be activated between two points.

[0120] To meet condition ③, it is necessary to calculate the external normal vector n of point P based on the above two conditions. P and the external normal vector n of point Q Q Taking point P as an example, its external normal vector n P The calculation formula is as follows:

[0121]

[0122] Similarly, the external normal vector n of point Q can be calculated Q After obtaining the external normal vectors at points P and Q, if they satisfy the following conditions at the same time:

[0123] n P ·n Q <0 (22)

[0124] n P ·l PQ >0 (23)

[0125] It is considered that points P and Q are not in the same bubble. In formula (23), l PQ is the unit vector pointing from P to Q.

[0126] After the above three conditions are met, the repulsive force generated by point Q on point P is:

[0127]

[0128] Where: K is a constant that controls the strength of the repulsive force; l PQ is the unit vector from P to Q, φ g represents the gas phase volume fraction, x P is the position coordinate of point P, x Q is the position coordinate of point Q, h PQ is the distance between points P and Q. Figure 1 As shown in the figure, the point Q that can produce a repulsive force on point P is usually not unique, and the total repulsive force F on point P is rep It can be expressed as:

[0129]

[0130] Where: n P is the external normal vector of the interface at point P.

[0131] Finally, the total repulsive force can be included in the force F to introduce the inter-foam repulsive force into the three-phase percolation model.

[0132] In the three-phase foam flow demonstration example based on the present invention, applying the wetting boundary includes the following:

[0133] In multiphase flow problems, the wettability of different fluids on solid walls is usually different. The contact angle formed between the interface between different fluids and the solid wall at the wall will strongly affect the dynamic behavior of the interface, resulting in significant differences in the interface morphology and the speed of the moving contact line. For multiphase flow simulations in confined spaces, accurate contact angle modeling can not only help capture the evolution of the phase interface near the solid wall over time, but also reflect the different wetting effects of different fluids on the solid wall. With the idea of ​​diffusion interface, the present invention intends to match the required contact angle size by assigning a virtual fluid volume fraction to the solid grid points close to the wall boundary (hereinafter referred to as solid boundary grid points).

[0134] like Figure 3 As shown in the figure, in any complex solid wall, the entire flow field is divided into two parts, the fluid domain and the solid domain, by the solid wall boundary. The grid points in the solid domain are recorded as solid grid points, and the remaining grid points are fluid grid points. The solid boundary grid point is defined as: for a certain solid grid point, if any grid point among the 8 grid points around it is a fluid grid point, it is regarded as a solid boundary grid point. Before matching the contact angle, the external normal vector of the solid boundary is calculated first. Taking a certain solid boundary grid point S as an example, its external normal vector n S The calculation formula is:

[0135]

[0136] Where: c j is the mesoscopic velocity in the eighth-order isotropic discrete space, x S is the position coordinate of point S; s(x) is the indicator function, which takes the value 1 if and only if x is a solid grid point, otherwise it takes the value 0; ω(c j 2 ) is the eighth-order weight coefficient, and its value is as follows:

[0137]

[0138] For any fluid grid point, theoretically, each grid point contains three fluids. Therefore, for any phase of fluid, it can form two contact angles with the other two phases of fluid at the solid wall. However, the implementation of the wetting boundary condition requires a unique contact angle, so it needs to be obtained in a weighted manner. At any fluid grid point, the sum of the volume fractions of the three fluids is 1, that is, φ o +φw +φ g = 1. Therefore, only the weighted contact angles θ of the oil phase and the water phase need to be calculated o and θ w , and its calculation formula is:

[0139] (1-φ o )θ o =φ w θ ow +φ g θ og (28)

[0140] (1-φ w )θ w =φ o θ wo +φ g θ wg (29)

[0141] Where: φ o is the volume fraction of oil phase, φ w is the volume fraction of water phase, φ g is the gas phase volume fraction, and the contact angle between the interface of fluids k and l and the solid wall measured from the fluid k side; where k and l represent any two phases of the three-phase fluids o, w, and g. Obviously, the contact angle θ kl and θ lk They are complementary angles, i.e. θ lk =π-θ kl .

[0142] like Figure 3 As shown, after obtaining the normal vector n outside point S S and the weighted contact angle θ k (k=o and w) then two lines can be drawn with n S The angles formed are all |π / 2-θ k |Characteristic line l C1 , l C2 , both of which intersect with the grid lines. The intersection of the characteristic line and the grid lines in the solid domain, such as point The first intersection points of the two characteristic lines in the fluid domain are S1 and S2, respectively, and it is necessary to ensure that the adjacent grid points of the intersection points S1 or S2 on the horizontal or vertical grid lines are all fluid grid points, such as points A, B or points C, D. If the adjacent grid points are not fluid grid points, the characteristic lines need to be extended until the above conditions are met.

[0143] Through flow field calculation, the fluid volume fractions at points A, B, C, and D are all known, so the fluid volume fractions at points S1 and S2 can be obtained by linear interpolation, that is, the volume fractions of fluid k at the intersection points S1 and S2 are:

[0144]

[0145] Where: represents the position coordinates of the intersection point S1, represents the position coordinates of the intersection point S2, x A Represents the position coordinates of point A, x B Indicates the position coordinates of point B, x C Indicates the position coordinates of point C, x D Represents the position coordinates of point D; represents the abscissa of the intersection point S1; represents the ordinate of the intersection point S2; x A represents the horizontal coordinate of point A, x B Indicates the horizontal coordinate of point B; y C Indicates the ordinate of point C, y D Represents the vertical coordinate of point D.

[0146] get and After that, it is necessary to judge the wall wettability. or Assign values ​​to the solid boundary grid points S:

[0147]

[0148] Where: θ k represents the weighted contact angle of fluid k.

[0149] After traversing all boundary solid grid points and executing the above calculation strategy on them, the volume fractions of the three-phase fluids of oil, water, and gas on all boundary solid grid points can be obtained and used to calculate the color gradient in equation (11), thereby affecting the interfacial tension F s In this way, the influence of the contact angle is accounted for.

[0150] In the process of oil displacement, water phase is usually injected first to continue displacement until no more oil can be produced, and then plugging and profile control technical measures are taken in the oil layer to further improve the recovery rate. Among them, foam has the characteristics of blocking large but not small, blocking water but not oil, and has the ability to reduce oil-water interfacial tension and change the wetting properties of rock wall, so it can well expand the affected area and improve oil displacement efficiency.

[0151] Figure 4 The porous media geometry and fluid distribution before water flooding are shown. The porous media calculation domain is obtained by filling isotropic, staggered cylinders, and three coarse, medium and fine crack channels are artificially constructed. The right side of the crack channel is the dense area; the red is the oil phase to be displaced, the green is the water phase, and the blue is the foam phase (not yet appeared). The displacement fluid is injected from the uniformly distributed holes on the left at a constant flow rate and flows out from the right boundary.

[0152] Figure 5 The curve of residual oil saturation in porous media changing with the increase of injected pore volume (PV) during the first step of water flooding to steady state (i.e., further flooding can not produce more crude oil) is shown. The main criterion for flooding to stability is that the decline rate of the residual oil saturation curve is less than a given small amount.

[0153] Figure 6 It shows the distribution patterns of oil and water phases in the porous medium after the first stage of displacement; this is also the starting point of the next stage of foam flooding.

[0154] Figures 7 to 11 Several key nodes after the second stage of foam injection are shown, and how the residual oil is displaced and produced. In particular, the foam fluid contains both gas and liquid phases. In order to inject uniform and size-controlled foam, a buffer layer needs to be set at the left inlet section and the porous medium area, and the foam needs to be injected in an alternating manner of water and gas injection.

[0155] Figure 7 It shows that in the initial stage of foam injection, the migration rate of foam along the wide channel is the fastest, while the migration along the narrow channel is blocked; in addition, the foam in the wide channel first enters the dense area on the right and gathers at the conical expansion throat, producing a blocking effect, prompting the water phase to open up new flow channels along the cone surface (see the area marked by the white rectangular box).

[0156] Figure 8 It shows that after a period of foam injection, the wide channels and medium channels are filled with foam; while the foam in the narrow channels is difficult to migrate, so it is retained in the narrow channels and produces a large flow resistance, which prompts the water phase to penetrate into the uppermost flaky residual oil and open up new flow channels (see the area marked by the white rectangular frame).

[0157] Fig. 9 It shows that after a period of continuous foam injection, the four flaky residual oils on both sides of the three dominant channels all showed different degrees of foam flow intrusion, and four new displacement fronts were generated (see the white line highlighted area). Among them, the upper side of the top flaky residual oil is a rigid wall, which produces the largest longitudinal flow resistance, while the foam migration in the lower narrow channel is the most difficult (that is, the best blocking effect), so the water phase displacement front moves the fastest, and the water phase and foam phase front separation phenomenon occurs locally; while the flaky residual oil in the lower middle part is blocked by foam on both sides, resulting in insufficient flow resistance, slow forward breakthrough rate of the water phase, and foam enrichment at the displacement front; in addition, the oblique front edge also proves that narrow channels have the ability to produce stronger foam flow resistance than wide channels.

[0158] Fig.10 The key node where pure water first breaks through the top flaky residual oil is shown (see the area marked by the white rectangular box).

[0159] Fig.11 The three-phase fluid morphology in the porous medium after pure water displacement at the top is shown. Among them, the uppermost sheet-like residual oil structure is completely destroyed, and the remaining drop-like residual oil is retained in the blind end or pore throat because it is close to the wall or wrapped by the newly generated branch-type dominant channel; similarly, the left side of the middle and upper sheet-like residual oil is also well displaced, thereby significantly increasing the affected area; while the residual oil in the middle and lower part cannot be effectively displaced because the foam on both sides cannot generate enough water flow resistance; finally, only pure water displacement channels appear in the bottom residual oil, but there is no foam infiltration, which is caused by the narrow new channel and the existence of a wide channel with low flow resistance next to it.

[0160] In general, Figures 6 to 11 The demonstration example shown proves that the three-phase foam flow simulation method described in the present invention can better simulate the whole process of injecting foam into water-containing reservoirs to improve crude oil recovery, and at the same time embodies the following three characteristics: ① It can effectively simulate the non-miscible three-phase flow behavior of oil, water and gas, and the interfacial tension effect in the foam generation and migration process is effectively captured; ② The introduction of repulsive force effectively curbs the non-physical aggregation of foam, so that aggregation will not occur in various situations such as foam aggregation and deformation; ③ The wetting boundary condition can accurately model the wettability of the three-phase fluid on the solid surface and reasonably characterize the movement of the contact line on the solid surface. In addition, this demonstration example well demonstrates the internal mechanism of foam injection "plugging and profile adjustment" to increase the swept area and enhance oil recovery, and for the first time reveals the migration law of the three-phase fluid in the foam displacement process from the perspective of numerical simulation. In summary, the method described in the present invention can effectively simulate the foam injection oil displacement process in high water-containing reservoirs.

[0161] The present invention organically combines a variety of independently developed algorithms and models, and can simulate the immiscible three-phase foam flow under the full range of interfacial tension conditions, reflecting many characteristics: conservation of mass of each phase, small pseudo-velocity, prevention of non-physical aggregation of foam, correct description of fluid wetting characteristics and near-wall solid-liquid flow characteristics, etc. Finally, based on these characteristics, this model can reasonably reveal the migration law of three-phase foam flow in non-uniform porous media, as well as the internal mechanism of foam-enhanced oil recovery, plugging and profile control.

[0162] The three-phase foam flow injection simulation method for enhanced oil recovery in high water content reservoirs provided by the present invention has the following advantages:

[0163] 1. The model itself has the ability to simulate three-phase flow under the full range of interfacial tension conditions and maintain advantages such as mass conservation; 2. The introduction of the repulsive force model can curb the non-physical aggregation behavior of foam in numerical simulation and correctly reproduce the behavior of foam in porous media; 3. In the framework of the diffusion interface, the wetting boundary conditions are imposed by assigning virtual fluid volume fractions to the solid boundary grid points to match the required contact angle size, which not only avoids the stress singularity caused by the conflict between the contact line movement and the no-slip boundary condition, but also can flexibly adjust the contact angle to accurately reflect the wetting characteristics of the fluid.

[0164] Based on the same inventive concept, an embodiment of the present invention further provides a three-phase foam flow injection simulation system for enhanced oil recovery in high water content reservoirs, comprising:

[0165] A model building module is used to build an immiscible three-phase flow calculation model under the lattice Boltzmann framework, and introduce a continuous surface force model and a recoloring algorithm into the immiscible three-phase flow calculation model; wherein the immiscible three-phase flow is composed of oil, gas and water; the continuous surface force model is used to generate interfacial tension between different phase interfaces, so that the model has the ability to capture the interface changes of droplets and bubbles in porous media; the recoloring algorithm is used to separate the phases of fluid along the interface normal and maintain the mass conservation of each phase to generate an identifiable interface;

[0166] The model correction module is used to introduce an additional source term on the right side of the distribution function equation when executing the collision step in the moment space, and introduce the inter-foam repulsion force into the immiscible three-phase flow calculation model in the form of volume force to correct the immiscible three-phase flow calculation model, wherein the inter-foam repulsion force is obtained by calculating the inter-foam repulsion force model, and is used to curb the non-physical aggregation of foams; after executing the collision step, the distribution function in the moment space is first restored by inverse transformation to obtain the total density distribution function, and then the density distribution function is assigned to the three different fluids using the re-coloring algorithm to ensure that different fluids are separated along the interface normal and the mass conservation of each phase fluid is maintained; the contact angle geometric formula is used to apply the wetting boundary condition at the fluid-solid boundary to reflect the different wetting characteristics of different fluids on the solid wall;

[0167] The simulation module is used to input the required three-phase fluid physical properties and provide the auxiliary parameters necessary for executing the simulation to run the immiscible three-phase flow calculation model, and obtain multiple sets of data characterizing the three-phase flow state during the foam injection flooding process through time advancement simulation.

[0168] Each module in the above-mentioned three-phase foam flow injection simulation system for enhanced oil recovery in high water-cut reservoirs can be implemented in whole or in part by software, hardware and a combination thereof. Each of the above-mentioned modules can be embedded in or independent of a processor in a computer device in the form of hardware, or can be stored in a memory in a computer device in the form of software, so that the processor can call and execute operations corresponding to each of the above modules.

[0169] The present invention also provides a computer device, including a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps in the embodiment of the three-phase foam flow injection simulation method for enhanced oil recovery in high water content reservoirs. The specific implementation method can be found in the method embodiment, which will not be repeated here.

[0170] Furthermore, the present invention also provides a non-temporary computer-readable storage medium containing instructions, and a computer program is stored on the storage medium. For example, a memory containing instructions, the above instructions can be executed by a processor of a computer device to complete the above method. For example, the non-temporary computer-readable storage medium can be a ROM, a random access memory (RAM), a CD-ROM, a magnetic tape, a floppy disk, and an optical data storage device. When the computer program is executed by the processor, the steps in the embodiment of the three-phase foam flow injection simulation method for enhanced oil recovery in high water-cut oil reservoirs can be implemented. The specific implementation method can be found in the method embodiment, which will not be repeated here.

[0171] It will be appreciated by those skilled in the art that embodiments of the present invention may provide methods, systems or computer program products. Therefore, the present invention may take the form of a complete hardware embodiment, a complete software embodiment or an embodiment combining software and hardware. Moreover, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0172] The present invention is described with reference to flowcharts and / or block diagrams of methods, apparatus (systems) and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowchart and / or block diagram, as well as a combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0173] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture including an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1A function specified in one or more boxes.

[0174] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process in the computer or other programmable device. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0175] It should be pointed out that the specific implementation methods described above can enable those skilled in the art to understand the invention more comprehensively, but do not limit the invention in any way. Therefore, although the invention has been described in detail in this specification and embodiments, those skilled in the art should understand that the invention can still be modified or replaced by equivalents; and all technical solutions and improvements that do not deviate from the spirit and scope of the invention are included in the protection scope of the patent for the invention. Any figure mark in the claims should not be regarded as limiting the claims involved. Any simple change or equivalent replacement of the technical solution that can be obviously obtained by any technician familiar with the field within the technical scope disclosed in the present invention belongs to the protection scope of the present invention.

Claims

1. A three-phase foam flow injection simulation method for enhanced oil recovery in high water content reservoirs, characterized in that: include: An immiscible three-phase flow computational model is constructed under the lattice Boltzmann framework, and a continuous surface force model and a recoloring algorithm are introduced into the immiscible three-phase flow computational model; wherein the immiscible three-phase flow is composed of oil, gas and water; the continuous surface force model is used to generate interfacial tension between different phase interfaces, so that the model has the ability to capture the interface changes of droplets and bubbles in porous media; the recoloring algorithm is used to separate the fluids of each phase along the interface normal and maintain the conservation of mass of each phase to generate an identifiable interface; When executing the collision step in the moment space, an additional source term is introduced on the right side of the distribution function equation, and the inter-foam repulsion force is introduced into the immiscible three-phase flow calculation model in the form of volume force to correct the immiscible three-phase flow calculation model, wherein the inter-foam repulsion force is calculated by the inter-foam repulsion force model to curb the non-physical aggregation of foams; after executing the collision step, the distribution function in the moment space is first restored by inverse transformation to obtain the total density distribution function, and then the density distribution function is assigned to the three different fluids using the re-coloring algorithm to ensure that different fluids are separated along the interface normal and the mass conservation of each phase of the fluid is maintained; the contact angle geometric formula is used to apply the wetting boundary condition at the fluid-solid boundary to reflect the different wetting characteristics of different fluids on the solid wall; The required three-phase fluid physical property parameters are input and the auxiliary parameters necessary for simulation are provided to run the immiscible three-phase flow calculation model. Multiple sets of data characterizing the three-phase flow state during foam injection flooding are obtained through time advancement simulation.

2. The three-phase foam flow injection simulation method for enhanced oil recovery in high water content reservoirs according to claim 1, characterized in that: The continuous surface force model is used to generate interfacial tension F between different phase interfaces. s , and its calculation formula is: in: Where: I is the unit matrix, k and l represent any two phases of the three-phase fluid of oil, water and gas; σ kl is the interfacial tension coefficient between fluid k and fluid l; ρ is the density of the three-phase mixed fluid, ρ k is the local density of fluid k, ρ l is the local density of fluid l; φ l is the volume fraction of fluid l, is the density of pure fluid l, φ k is the volume fraction of fluid k, is the density of pure fluid k, n kl is the unit external normal vector of the kl interface, G kl is the color gradient, C kl is the parameter that controls the interfacial tension strength of the kl interface.

3. The three-phase foam flow injection simulation method for enhanced oil recovery in high water content reservoirs according to claim 1, characterized in that: The inter-foam repulsive force is calculated by the inter-foam repulsive force model. For the inter-foam repulsive force between two activated repulsive force points P and Q that are moving away from each other, the specific calculation is as follows: The repulsive force generated by point Q on point P is: Where: K is a constant that controls the strength of the repulsive force; l PQ is the unit vector from P to Q, φ g represents the gas phase volume fraction, x P is the position coordinate of point P, x Q is the position coordinate of point Q, h PQ is the distance between points P and Q; The total repulsive force F on point P rep It is expressed as: Where: n P is the external normal vector of point P; The total repulsive force F rep By taking the repulsive force between bubbles into account in the resultant force F, the repulsive force between bubbles can be introduced into the three-phase seepage model.

4. The three-phase foam flow injection simulation method for enhanced oil recovery in high water content reservoirs according to claim 1, characterized in that: The recoloring operation of the recoloring algorithm is performed on the density distribution function in the velocity space. Therefore, before the recoloring operation, the post-collision velocity moment m′ is first inversely transformed to obtain the post-collision total density distribution function f′, that is, f′(x, t)=M -1 m′(x,t), where M is the orthogonal transformation matrix; then, the total density distribution function after the collision is re-colored: Where: x represents the spatial position, t is the time, and f i,k ″ represents the distribution function of fluid k after re-coloring, e i is the grid velocity of the D2Q9 model, w i for e i Directional weight coefficient; ρ k is the local density of fluid k, ρ l is the local density of fluid l, ρ is the total density of the three-phase mixed fluid, n kl is the unit external normal vector of the kl interface, β kl is a separation parameter related to the fluid kl interface thickness and satisfies β kl =β lk .

5. The three-phase foam flow injection simulation method for enhanced oil recovery in high water content reservoirs according to claim 1, characterized in that: The contact angle geometric formula is used to apply the wetting boundary condition at the fluid-solid boundary, specifically including: In any local part of the solid wall, the computational domain is divided into two parts, the fluid domain and the solid domain, by the solid wall boundary. The grid points in the solid domain are recorded as solid grid points, and the remaining grid points are fluid grid points. The solid boundary grid points are defined as follows: for a solid grid point, if any grid point among the 8 grid points around it is a fluid grid point, it is regarded as a solid boundary grid point. Before matching the contact angle, the external normal vector of the solid boundary is calculated first. Assuming that a solid boundary grid point is S, its external normal vector n S The calculation formula is: Where: c j is the mesoscopic velocity in the eighth-order isotropic discrete space, x S is the position coordinate of point S; s(x) is the indicator function, which takes the value 1 if and only if x is a solid grid point, otherwise it takes the value 0; ω(c j 2 ) is the eighth-order weight coefficient, and its value is as follows: At any fluid grid point, the sum of the volume fractions of the three fluids is 1, that is, φ o +φ w +φ g =1; therefore, only the weighted contact angles θ of the oil phase and the water phase need to be calculated o and θ w , and its calculation formula is: (1-φ o )i o =φ w i ow +φ g i og ; (1-φ w )i w =φ o i wo +φ g i wg ; Where: k, l represent any two phases of the three-phase fluid o, w, g, φ o is the volume fraction of oil phase, φ w is the volume fraction of water phase, φ g is the gas phase volume fraction; θ kl is the contact angle between the interface of fluids k and l and the solid wall measured from the side of fluid k; contact angle θ kl and θ lk They are complementary angles, i.e. θ lk =π-θ kl ; After obtaining the normal vector n outside point S S and the weighted contact angle θ k After that, two lines can be drawn with n S The angles formed are all |π / 2-θ k |Characteristic line l C1 , l C2 , both intersect with the grid lines respectively; the intersection point between the characteristic line and the grid line in the solid domain is an invalid intersection point; the first intersection points of the two characteristic lines in the fluid domain are S1 and S2 respectively, and it is necessary to ensure that the adjacent grid points of the intersection points S1 or S2 on the horizontal or vertical grid lines are all fluid grid points; if the adjacent grid points are not fluid grid points, it is necessary to continue to extend the characteristic line until the above conditions are met; Assuming that A and B are the fluid grid points closest to the intersection S1 in the horizontal direction, and C and D are the fluid grid points closest to the intersection S2 in the vertical direction, the fluid volume fractions at points A, B, C, and D are known through flow field calculation, so the fluid volume fractions at the intersections S1 and S2 can be obtained by linear interpolation. Then the volume fractions of fluid k at the intersections S1 and S2 are: Where: represents the position coordinates of the intersection point S1, Indicates the position coordinates of the intersection point S2, x A Represents the position coordinates of point A, x B Indicates the position coordinates of point B, x C Indicates the position coordinates of point C, x D Represents the position coordinates of point D; represents the horizontal coordinate of the intersection S1; x A Indicates the horizontal coordinate of point A; x B represents the horizontal coordinate of point B, represents the ordinate of the intersection point S2; y C Indicates the ordinate of point C, y D represents the ordinate of point D; get and After that, the wall wettability is determined. or Assign values ​​to the solid boundary grid points S: Where: θ k represents the weighted contact angle of fluid k; After traversing all boundary solid grid points and executing the above calculation strategy, the volume fractions of oil, water, and gas three-phase fluids matching the required contact angles on all boundary solid grid points can be obtained, thereby affecting the color gradient and interface normal vector on the boundary fluid grid points; in this way, the wetting boundary conditions are implicitly applied.

6. A three-phase foam flow injection simulation system for enhanced oil recovery in high water content reservoirs, characterized in that: include: A model building module is used to build an immiscible three-phase flow calculation model under the lattice Boltzmann framework, and introduce a continuous surface force model and a recoloring algorithm into the immiscible three-phase flow calculation model; wherein the immiscible three-phase flow is composed of oil, gas and water; the continuous surface force model is used to generate interfacial tension between different phase interfaces, so that the model has the ability to capture the interface changes of droplets and bubbles in porous media; the recoloring algorithm is used to separate the phases of fluid along the interface normal and maintain the mass conservation of each phase to generate an identifiable interface; The model correction module is used to introduce an additional source term on the right side of the distribution function equation when executing the collision step in the moment space, and introduce the inter-foam repulsion force into the immiscible three-phase flow calculation model in the form of volume force to correct the immiscible three-phase flow calculation model, wherein the inter-foam repulsion force is obtained by calculating the inter-foam repulsion force model, and is used to curb the non-physical aggregation of foams; after executing the collision step, the distribution function in the moment space is first restored by inverse transformation to obtain the total density distribution function, and then the density distribution function is assigned to the three different fluids using the re-coloring algorithm to ensure that different fluids are separated along the interface normal and the mass conservation of each phase fluid is maintained; the contact angle geometric formula is used to apply the wetting boundary condition at the fluid-solid boundary to reflect the different wetting characteristics of different fluids on the solid wall; The simulation module is used to input the required three-phase fluid physical properties and provide the auxiliary parameters necessary for executing the simulation to run the immiscible three-phase flow calculation model, and obtain multiple sets of data characterizing the three-phase flow state during the foam injection flooding process through time advancement simulation.

7. A computer device comprising a memory, a processor and a computer program stored in the memory, characterized in that: The processor executes the computer program to implement the steps of the method according to any one of claims 1 to 5.

8. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is loaded into a processor, it can execute the steps of the method according to any one of claims 1 to 5.

Images