A three-phase foam flow injection simulation method for enhanced oil recovery in high water cut reservoirs

Abstract

The application provides a three-phase foam flow injection simulation method for enhanced oil displacement in high water cut reservoirs, belongs to the field of computational fluid dynamics, constructs a non-miscible three-phase flow calculation model under a lattice Boltzmann framework, and introduces a continuous surface force model and a re-labeling algorithm in the non-miscible three-phase flow calculation model; repulsive forces are added between required gas-phase foams, so that the three-phase flow model has the ability to simulate non-coalesced foam behavior; wetting boundary conditions are applied according to the diffusion interface idea to process solid wall boundaries with specific wettability under complex geometric conditions; non-miscible three-phase flow calculation models are run by inputting three-phase fluid physical property parameters and providing necessary auxiliary parameters for simulation, and a plurality of groups of data representing three-phase flow states in the foam injection oil displacement process are obtained through time advancement simulation. The application fills the gap of the lattice Boltzmann method in simulating oil, water and gas three-phase foam flow, and can perform high-precision simulation on the three-phase foam injection oil displacement process.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of computational fluid dynamics, and particularly relates to a three-phase foam flow injection simulation method for enhanced oil displacement in high water cut reservoirs. BACKGROUND

[0002] At present, the oil fields that have been exploited in China have generally entered the high water cut stage, and the development potential of oil fields is gradually declining. The recovery effect of the traditional water flooding oil recovery method is poor. In order to maintain the crude oil production and prolong the development life of old oil fields, it is urgent to innovate and improve the recovery technology. Foam, as a fluid with high apparent viscosity, has the characteristics of stability when encountering water, and thus can maintain a stable state in high permeability parts with more water. In addition, foam generally exhibits high seepage resistance, and the seepage resistance increases with the increase of permeability, so it can effectively plug high permeability layers, improve the 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 at the same time play the role of emulsification and oil displacement, so 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 experimental research cost is high, and at the same time, it is difficult to accurately obtain important flow field information such as foam shape, velocity and pressure under three-dimensional conditions due to the limitation of observation means, which makes the current understanding of the three-phase seepage mechanism of foam flooding still not comprehensive. In order to make up for the deficiency of experimental research, and to deepen the understanding of the migration law of foam fluid in high water cut reservoirs in an economical, efficient and intuitive way, using high-precision simulation method 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 process of foam fluid plugging profile control for enhanced oil displacement, and according to the simulation results, the best working condition parameters for enhanced recovery are obtained.

[0004] In the field of numerical simulation, the numerical simulation of two-phase flow has been studied for a long time, and the related literatures have covered many subfields, such as flow regime analysis, heat and mass transfer, pressure drop characteristics, and phase distribution. In contrast, the study of three-phase flow is more challenging due to the complex interactions between the fluids and the changes in the interface topology. In particular, in immiscible three-phase flow, the singularity of the interface curvature in the coexistence region and the topological changes of the different fluid interfaces are the main difficulties in numerical modeling. Currently, the interface modeling methods are generally divided into sharp interface and diffuse interface methods. 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. The implementation strategy of the diffuse interface method is to allow the mixing of different fluid components at the interface to form a thin diffusion layer, which realizes the smooth transition of the fluid property parameters, and introduces a marker function to distinguish different fluids. The main advantage of the diffuse interface method is that it can avoid the difficulty of modeling the interfacial tension caused by the singularity of the interface curvature in the coexistence region, and has the ability to capture complex interface topological changes, such as bubble breakup and coalescence.

[0005] When simulating multiphase flow, the way of handling boundary conditions can significantly affect the simulation results. Traditional multiphase flow simulation methods usually strictly apply no-slip boundary conditions at the wall surface, which can easily lead to non-physical phenomena such as the inability of the three-phase contact line to move, or stress singularities near solid walls, which can cause numerical simulation divergence.

[0006] The lattice Boltzmann color model is a numerical model developed by introducing a color indicator function and its gradient in the lattice Boltzmann framework. The color indicator function is used to represent different phases, and its gradient is used to represent the interfacial tension phenomenon in multiphase flow. The re-coloring operator plays a key role in the color model, which can promote the separation of different fluids along the interface normal at the mesoscopic level, thus ensuring the conservation of mass of different fluids in the simulation process, obtaining a relatively sharp interface, and ensuring 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 techniques, and can automatically describe the dynamic behaviors of interface deformation, breakup, and coalescence. In addition, the color model inherits many advantages of the lattice Boltzmann method, including natural parallelism, easy code writing, and easy handling of complex geometric boundaries, so it has significant advantages in porous-scale multiphase flow simulation and complex interface dynamic behavior research.

[0007] Researchers have developed various three-phase flow solvers under the lattice Boltzmann framework, which focus on interface tension modeling and improving the model's performance in handling large ranges of fluid density ratios and viscosity ratios. However, these solvers do not have the ability to simulate "non-coalescing" foam flow, and therefore do not have the ability to simulate three-phase flow in the process of foam injection enhanced oil recovery. SUMMARY

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

[0009] To achieve the above purpose, the application provides the following technical scheme:

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

[0011] A non-miscible three-phase flow calculation model is constructed under the lattice Boltzmann framework, and a continuous surface force model and a re-coloring algorithm are introduced into the non-miscible three-phase flow calculation model; wherein the non-miscible 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 interfacial changes of droplets and bubbles in porous media; the re-coloring algorithm is used to separate each phase fluid along the interface normal and maintain the mass conservation of each phase, and an identifiable interface is generated;

[0012] When performing the collision step in the moment space, an additional source term is introduced at the right end of the distribution function equation, and the repulsion force between foams is introduced into the non-miscible three-phase flow calculation model in the form of volume force to modify the non-miscible three-phase flow calculation model, wherein the repulsion force between foams is obtained by calculating the repulsion force between foams, which is used to suppress the non-physical coalescence of foams; after the collision step is performed, the total density distribution function is restored by inverse transformation in the moment space, and then the density distribution function is allocated to the three different fluids by using the re-coloring algorithm, so as to ensure that different fluids are separated along the interface normal and maintain the mass conservation of each phase fluid; the wetting boundary condition is applied at the fluid-solid boundary by using the contact angle geometric formula, so as to reflect the different wetting characteristics of different fluids on the solid wall surface;

[0013] The required three-phase fluid physical property parameters are input, and the necessary auxiliary parameters for performing simulation are provided to run the non-miscible three-phase flow calculation model, and a plurality of groups of data representing the state of three-phase flow in the foam injection oil displacement process are obtained by time advancement simulation.

[0014] Preferably, the continuous surface force model is used to generate interfacial tension F s , and the calculation formula is:

[0015]

[0016] Wherein:

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

[0018]

[0019] where I is the identity matrix, k, l, o, w, g denote any two phases in the three-phase fluid; σ kl is the interfacial tension coefficient between fluid k and fluid l; p is the density of the three-phase mixture, p k is the local density of fluid k, p 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 outer normal vector of the k-l interface, G kl is the color gradient, C kl is a parameter controlling the interfacial tension strength of the k-l interface.

[0020] Preferably, the inter-bubble repulsive force is obtained by an inter-bubble repulsive force model, for two activated repulsive force points P and Q, the inter-bubble repulsive force generated by each other is calculated as follows:

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

[0022]

[0023] where K is a constant controlling the strength of the repulsive force; l PQ is the unit vector from P to Q, φ g denotes the volume fraction of gas phase, 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 rep experienced by point P is:

[0025]

[0026] where n P is the outer normal vector of point P;

[0027] The total repulsive force F rep is taken into the resultant force F, which can introduce the inter-bubble repulsive force into the three-phase percolation model.

[0028] Preferably, the re-coloring operation of the re-coloring algorithm is performed in the velocity space for the density distribution function, so before the re-coloring operation, the post-collision momentum m' is transformed into the post-collision total density distribution function f' through an inverse transformation, i.e. f'(x, t) = M -1 m'(x, t), wherein M is an orthogonal transformation matrix; subsequently, the re-coloring operation is performed on the post-collision total density distribution function:

[0029]

[0030] wherein: x represents a spatial position, t is time, f i,k represents the distribution function of fluid k after re-coloring, e i is the lattice velocity of the D2Q9 model, w i is the weight coefficient in the e i direction; p k is the local density of fluid k, p l is the local density of fluid l, p is the total density of the three-phase mixed fluid, n kl is the unit outer normal vector of the k-l interface, b kl is a separation parameter related to the thickness of the fluid k-l interface, and satisfies b kl = b lk .

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

[0032] At any local solid wall surface, the calculation domain is divided into a fluid domain and a solid domain by the solid wall boundary, the lattice points in the solid domain are recorded as solid lattice points, and the remaining lattice points are fluid lattice points; the solid boundary lattice point is defined as follows: for a certain solid lattice point, if any of the eight lattice points around it is a fluid lattice point, it is regarded as a solid boundary lattice point; before matching the contact angle, the outer normal vector of the solid boundary is first calculated; assuming that a certain solid boundary lattice point is S, the calculation formula of the outer normal vector n S of S is as follows:

[0033]

[0034] wherein: c j is the eight-order isotropic mesoscopic velocity in discrete space, x S is the position coordinate of S; s(x) is an indicator function, which takes 1 only when x is a solid lattice point, and otherwise takes 0; w(c j 2 is the eight-order weight coefficient, and the values are as follows:

[0035]

[0036] The sum of the volume fractions of the three fluids at any fluid grid point is 1, i.e. φ o + φ w + φ g = 1; therefore, only the weighted contact angles θ o and θ w of the oil phase and the water phase need to be calculated, and the calculation formulae are:

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

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

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

[0040] After obtaining the outward normal vector n S of the S point and the weighted contact angle θ k , two characteristic lines l C1 and l C2 with an included angle of | π / 2- θ k | with n S can be drawn, and the two characteristic lines intersect with the grid lines, respectively; the intersection point of 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 required to ensure that the adjacent grid points of the intersection point S1 or S2 located on the horizontal or vertical grid line are fluid grid points; if the adjacent grid points are not fluid grid points, the characteristic line needs to be continued to be extended until the above condition is met;

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

[0042]

[0043] In the formula: x represents the position coordinate of intersection S1, x represents the position coordinate of intersection S2, A x represents the position coordinate of point A, B x represents the position coordinate of point B, C x represents the position coordinate of point C, D x represents the position coordinate of point D; x represents the horizontal coordinate of intersection S1; x represents the vertical coordinate of intersection S2; A x represents the horizontal coordinate of point A, B x represents the horizontal coordinate of point B; C x represents the vertical coordinate of point C, D x represents the vertical coordinate of point D;

[0044] obtained and After that, according to the wall wetting judgment, the value of or is assigned to the solid boundary grid point S:

[0045]

[0046] In the formula: θ k x represents the weighted contact angle of fluid k;

[0047] After traversing all boundary solid grid points and performing the above calculation strategy, the oil, water and gas three-phase fluid volume fractions matching the required contact angle on all boundary solid grid points can be obtained, and then the color gradient and interface normal vector on the boundary fluid grid points are affected; in this way, the wetting boundary condition is implicitly applied.

[0048] The application also provides a three-phase foam flow injection simulation system for enhanced oil displacement in a high-water-content oil reservoir, comprising:

[0049] The model construction module is configured to construct a non-miscible three-phase flow calculation model under a lattice Boltzmann framework, and introduce a continuous surface force model and a re-coloring algorithm into the non-miscible three-phase flow calculation model; wherein the non-miscible three-phase flow is composed of oil, gas and water; the continuous surface force model is configured to generate an interfacial tension between different phase interfaces, so that the model has the ability to capture the interfacial changes of liquid droplets and gas bubbles in a porous medium; and the re-coloring algorithm is configured to separate each phase fluid along the normal direction of the interface, maintain the mass conservation of each phase fluid, and generate an identifiable interface.

[0050] The model correction module is configured to introduce an additional source term into the right end of a distribution function equation when performing a collision step in a moment space, introduce an inter-foam repulsion force into the non-miscible three-phase flow calculation model in the form of a volume force to correct the non-miscible three-phase flow calculation model, wherein the inter-foam repulsion force is obtained by calculating an inter-foam repulsion force model, and is configured to suppress the non-physical coalescence of foams; after the collision step is performed, the total density distribution function is restored from the distribution function in the moment space through inverse transformation, and the density distribution function is allocated to three different fluids by using the re-coloring algorithm, so as to ensure that different fluids are separated along the normal direction of the interface, and the mass conservation of each phase fluid is maintained; and the wetting boundary condition is applied at the fluid-solid boundary by using a contact angle geometric formula, so as to reflect the different wetting characteristics of different fluids on the solid wall surface.

[0051] The simulation module is configured to input required three-phase fluid physical property parameters and provide auxiliary parameters necessary for performing simulation to run the non-miscible three-phase flow calculation model, and obtain a plurality of groups of data representing the state of the three-phase flow in the process of foam injection and oil displacement by time advancement simulation.

[0052] The application further 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 the steps of any one of the three-phase foam flow injection simulation methods for enhanced oil displacement in high water cut reservoirs.

[0053] The application further provides a computer readable storage medium, wherein the storage medium stores a computer program, and the computer program can execute the steps of any one of the three-phase foam flow injection simulation methods for enhanced oil displacement in high water cut reservoirs when loaded by a processor.

[0054] The three-phase foam flow injection simulation method for enhanced oil displacement in high water cut reservoirs provided by the application has the following beneficial effects:

[0055] The application introduces a continuous surface force model and a re-coloring algorithm in a non-miscible three-phase flow calculation model, so that the model itself has the ability to simulate three-phase flow under full-range interfacial tension conditions, and maintains the advantages of fluid component mass conservation; the introduction of the repulsive force model can suppress the non-physical coalescence behavior of foam in numerical simulation, correctly reproduce the behavior of foam in porous media; under the diffusion interface framework, the contact angle geometry formula is used to apply the wetting boundary condition, avoiding stress singularity, being able to flexibly control the contact angle, and correctly reflecting the fluid wettability. Based on the above application, the non-miscible three-phase flow calculation model can accurately simulate the injection of three-phase foam flow. BRIEF DESCRIPTION OF DRAWINGS

[0056] In order to more clearly illustrate the embodiments of the application and the design scheme thereof, the following will briefly introduce the drawings required by the embodiments. The drawings in the following description are only part of the embodiments of the application, and other drawings can be obtained by those skilled in the art without creative labor on the basis of these drawings.

[0057] Figure 1 The flow chart of the three-phase foam flow injection simulation method for enhanced oil displacement in high water cut reservoirs of embodiment 1 of the application;

[0058] Figure 2 The schematic diagram of the repulsive force range between two close contact foams;

[0059] Figure 3 The schematic diagram of the implementation of the wetting boundary condition of the complex solid wall surface;

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

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

[0062] Figure 6 The distribution diagram of oil and water phases in the porous medium after water flooding is stable;

[0063] Figure 7 After injecting the foam, the foam aggregation phenomenon first appears in the wide channel expansion throat, demonstrating the plugging profile control effect of blocking the dominant channel and prompting the water phase to open a new flow channel;

[0064] Figure 8 After injecting the foam, the foam is blocked and stays in the narrow channel;

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

[0066] Figure 10 The pure water displacement phenomenon occurs at the top of the porous medium;

[0067] Figure 11 The schematic diagram for the complete destruction of the sheet residual oil structure at the top of the porous medium. DETAILED DESCRIPTION

[0068] In order to make the technical personnel of the present application better understand the technical solutions and can be implemented, the present application is described in detail below in conjunction with the drawings and specific examples. The following examples are only used to more clearly illustrate the technical solutions of the present application, and cannot limit the protection scope of the present application.

[0069] Example 1

[0070] In the study of three-phase flow, the presence of additional phases will lead to more complex fluid-fluid interaction and interface topological structure change, which makes the modeling and simulation of the interface dynamics process extremely complex. Especially in the immiscible three-phase flow, the singularity of the interface curvature in the three-phase coexistence region and the topological change of different fluid interfaces are the main difficulties of numerical modeling. In order to realize the simulation of the process of strengthening oil recovery by injecting foam into high water cut oil reservoirs in complex porous media, the present application aims to correct the lattice Boltzmann immiscible three-phase flow calculation model by introducing the repulsive force between the foams, and to apply the wetting boundary condition by using the contact angle geometric formula under the diffusion interface framework.

[0071] However, in numerical simulation, the foams usually show a tendency to coalesce, which is not consistent with the separated state of the foams under the comprehensive influence of van der Waals attraction and repulsive force generated by the surfactant adsorption layer in the actual working condition, so it is necessary to introduce repulsive force between the foams to maintain their non-coalescence state. Based on this, the present application proposes a simulation algorithm suitable for immiscible three-phase foam flow at pore scale to fill the research gap in the field, and lays a model foundation for mining and understanding the migration mechanism of foam flow in porous media.

[0072] Specifically, the present application provides a three-phase foam flow injection simulation method for strengthening oil displacement in high water cut oil reservoirs, specifically a simulation method suitable for immiscible oil, water and gas three-phase foam flow at pore scale under the lattice Boltzmann framework, as shown in the figure, the method comprises: Figure 1

[0073] ​A non-miscible three-phase flow calculation model is constructed under the lattice Boltzmann framework, and a continuous surface force model and a re-coloring algorithm are introduced into the non-miscible three-phase flow calculation model; wherein the non-miscible 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 interfacial changes of droplets and bubbles in a porous medium; the re-coloring algorithm plays a key role in the color model, and is used to separate each phase fluid along the interface normal, maintain the mass conservation of each phase, and ensure that a reasonable thickness of the diffuse interface is generated under the condition of full-range interfacial tension.

[0074] When performing a collision step in the moment space, an additional source term is introduced at the right end of the distribution function equation, and the repulsive force between foams is introduced into the non-miscible three-phase flow calculation model in the form of a volume force to modify the non-miscible three-phase flow calculation model, wherein the repulsive force between foams is obtained by calculating a repulsive force model between foams, and is used to suppress the non-physical coalescence of foams; after the collision step is performed, the total density distribution function in the moment space is first restored by inverse transformation to obtain the total density distribution function, also known as the velocity space total density distribution function; then the re-coloring algorithm is used to allocate the density distribution function to different fluids, to ensure that different phases separate along the interface normal, and to maintain mass conservation.

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

[0076] Under the framework of the diffuse interface, a contact angle geometric formula is used to apply a wetting boundary condition at the fluid-solid boundary to reflect the different wetting characteristics of different fluids on the solid wall surface.

[0077] The main advantages of the diffuse interface are: avoiding the modeling of interfacial tension caused by the singular interfacial curvature of the three-phase fluid coexistence region; having the ability to capture the complex interfacial topological changes, such as the breakup and coalescence of droplets. The diffuse interface method allows the contact line to move by introducing a diffusion mechanism, so the no-slip boundary condition can still be used when dealing with wetting boundary conditions, which not only greatly simplifies the implementation of wetting boundary conditions under irregular wall conditions, but also avoids the introduction of parameters such as slip model and slip length, which depend on the micro properties of materials and fluids and are difficult to measure.

[0078] The required three-phase fluid physical property parameters are input, and auxiliary parameters necessary for executing the simulation are provided to run the non-miscible three-phase flow calculation model, and a plurality of groups of data characterizing the state of the three-phase flow in the process of foam injection and oil displacement are obtained through time advancement simulation. Among them, the physical property parameters include fluid viscosity, interfacial tension coefficient of two fluids, three-phase contact angle, etc.; the auxiliary parameters include a D2Q9 velocity set, a transformation matrix, a diagonal relaxation matrix, etc.

[0079] Further, the non-miscible three-phase flow calculation model comprises the following contents:

[0080] First, three different non-miscible fluids are defined in the model. In the foam oil displacement system, the three-phase fluid is composed of oil, water and gas, and their density distribution functions are f o,i , f w,i , f g,i , respectively, where the subscripts o, w and g represent the corresponding fluids (as follows), and i is the lattice velocity index. In the embodiment of the present application, a two-dimensional nine-velocity (D2Q9) lattice model is used, so the value range of i is 0-8.

[0081] The total density distribution function f i 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 the collision operation on the total density distribution function, in order to improve the numerical stability of the simulation, a multi-relaxation collision model is used, so the collision process is performed in the moment space. Denote f = (f0, f1, …, f8) T , and the corresponding velocity moment is m = Mf, where M is an orthogonal transformation matrix, and the explicit expression is as follows:

[0084]

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

[0086]

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

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

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

[0090] where d7, d8 are terms related to the dynamic viscosity μ of the three-phase mixture, and d7 = d8; for symmetry reasons, it is necessary to satisfy d4 = d6. In the present model, in order to obtain a viscosity independent wall position while taking into account numerical stability, the relaxation parameters are chosen as follows:

[0091]

[0092] where ρ = ρ o + ρ w + ρ g is the density of the three-phase mixture. The dynamic viscosity μ can be obtained by harmonic averaging:

[0093]

[0094] The equilibrium function m eq in the moment space is obtained by transforming the equilibrium distribution function f in the velocity space, which is given by:

[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 F in the moment space can also be obtained by a similar transformation, which is given by:

[0098]

[0099] where F x and F y are the components of the force F in the horizontal and vertical directions, respectively. It is particularly pointed out that the force F in the embodiments of the present application is the resultant force of the interfacial tension F s and the repulsive force F rep between the bubbles. For a three-phase fluid of oil, water, and gas, the formula for calculating the interfacial tension F s is:

[0100]

[0101] where:

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

[0103]

[0104] In the above formula: k and l can traverse o, w, g respectively; I is a unit matrix, k, l represent any two phases in the three-phase fluid; σ kl is the interfacial tension coefficient between fluid k and fluid l; p is the density of the three-phase mixed fluid, p k is the local density of fluid k, p 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 outer normal vector of the k-l interface, G kl is the color gradient, C kl is a parameter for controlling the interfacial tension strength of the k-l interface.

[0105] Although F s can generate interfacial tension effect, but cannot ensure phase separation, so it is necessary to introduce a re-coloring operator to ensure that oil, water and gas are mutually insoluble. The re-coloring operation is performed on the density distribution function in the velocity space, so before the re-coloring operation, the total density distribution function f′ after collision is obtained by inverse transformation of m′, that is, f′(x, t) = M -1 m′(x, t), where M is an orthogonal transformation matrix. Subsequently, the re-coloring operation is performed on the total density distribution function after collision:

[0106]

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

[0108]

[0109] In the formula: 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 Equation (19), and enters 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:

[0116] As mentioned above, the introduction of the repulsive force term is intended to curb the non-physical aggregation between foams. In a three-phase foam flow system, a surfactant is introduced as a foaming agent and adsorbed on the foam surface. When two independent foams 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 foams 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 dashed circles in the upper left and lower right corners represent the influence ranges of points P and Q, respectively. The radius of the circles is the critical value for generating repulsive forces. In addition to point P, point Q can also influence points 1, 2, and 3; similarly, in addition to point Q, point P can also influence 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 in close contact. Two activated repulsive points P and Q will generate a repulsive force directed toward each other. Activating the repulsive force requires three conditions: ① Both points P and Q are within the diffusion interface region of their respective bubbles; ② The distance between points P and Q is less than a certain critical value; and ③ Points P and Q belong to different grid points within different bubbles. If all three conditions are not met, the repulsive force is not activated and takes the default value of 0.

[0118] In order to satisfy condition ①, it is necessary to examine the gas phase 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 apart 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 both of them are satisfied:

[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 the 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, 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 normal vector outside the interface at point P.

[0131] Finally, the inter-foam repulsive force can be introduced into the three-phase percolation model by including the total repulsive force into the force F.

[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 lattice points near the wall boundary (hereinafter referred to as solid boundary lattice points).

[0134] like Figure 3 As shown in the figure, at 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 follows: 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 a weighted approach is needed to obtain it. 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 is 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) after that, two lines can be made with n S The angles formed are all |π / 2-θ k |Characteristic line l C1 、l C2 , and they intersect with the grid lines respectively. 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. 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 and B or points C and 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] wherein: represents the position coordinate of intersection point S1, represents the position coordinate of intersection point S2, x A represents the position coordinate of point A, x B represents the position coordinate of point B, x C represents the position coordinate of point C, x D represents the position coordinate of point D; represents the abscissa of intersection point S1; represents the ordinate of intersection point S2; x A represents the abscissa of point A, x B represents the abscissa of point B; y C represents the ordinate of point C, y D represents the ordinate of point D.

[0146] obtained and after that, it is necessary to judge the wall wetting property to assign or to the solid boundary grid point S:

[0147]

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

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

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

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

[0152] Figure 5 The first step of water flooding to the steady state (i.e. further displacement of crude oil cannot be produced) process, the residual oil saturation in the porous medium with the increase of injection pore volume (PV) change curve, the main criterion for displacement to stable is the residual oil saturation curve of the decline rate is less than a given small amount.

[0153] Figure 6 The first step of water flooding to the steady state (i.e. further displacement of crude oil cannot be produced) process, the residual oil saturation in the porous medium with the increase of injection pore volume (PV) change curve, the main criterion for displacement to stable is the residual oil saturation curve of the decline rate is less than a given small amount.

[0154] Figure 7 to Figure 11 The second stage of injection of foam after a number of key nodes, and show how the residual oil is displaced. It is pointed out that the foam fluid contains gas and liquid phases, in order to inject uniform, size controllable foam, need to set the buffer layer in the left side of the entrance and porous media area, and with water gas injection method of injection of foam.

[0155] Figure 7 The initial stage of injection of foam, the foam along the wide channel migration rate is the fastest, and along the narrow channel migration is blocked; in addition, the foam in the wide channel first into the right side of the dense area, and gather in the tapered expansion throat, produce plugging effect, make water phase along the cone to open up new flow channel (see white rectangular frame marked area).

[0156] Figure 8 The injection of foam for a period of time, the wide channel and the channel are filled with foam; and the narrow channel foam migration difficult, thus retained in the larger flow resistance, prompting the water phase into the uppermost sheet residual oil, open up new flow channel (see white rectangular frame marked area).

[0157] Figure 9 The injection of foam for a period of time, the three dominant channel two sides of the four sheet residual oil appear different degree of foam flow invasion phenomenon, and new four displacement front (see white line highlighted area). Among them, the top sheet residual oil on the side of the rigid wall, the largest longitudinal flow resistance, and the lower part of the narrow channel foam migration is the most difficult (i.e. the best plugging effect), so the water phase displacement front moves the fastest, and local water phase, foam phase front separation phenomenon; and the lower sheet residual oil due to both sides of the foam plugging, resulting in flow resistance, water phase breakthrough rate is slow, foam enrichment in the displacement front; in addition, the diagonal front also proved that the narrow channel has stronger foam flow resistance than the wide channel.

[0158] Figure 10 The key node of the top sheet residual oil first appears pure water breakthrough (see white rectangular frame marked area).

[0159] Figure 11 The three-phase fluid morphology in the porous medium after further displacement evolution after the appearance of pure water displacement at the top is shown. Among them, the uppermost sheet-shaped residual oil structure is completely destroyed, and the remaining drop-shaped residual oil is retained in the blind end or pore throat due to being close to the wall surface or being wrapped by newly generated branched dominant channels; similarly, the sheet-shaped residual oil on the left side of the upper part is also well displaced, thereby significantly increasing the swept area; the middle and lower part of the residual oil cannot be effectively displaced due to the inability of the foam on both sides to generate sufficient water flow resistance; finally, only pure water displacement channels appear in the bottom residual oil, but no foam penetrates, which is caused by the narrower new channels and the wide channels with low flow resistance next to them.

[0160] Overall, Figure 6 to Figure 11 The demonstration example shown demonstrates that the three-phase foam flow simulation method described in the present application can better simulate the whole process of injecting foam into a water-containing reservoir to improve oil recovery, while embodying the following three characteristics: ① it can effectively simulate the non-miscible three-phase percolation 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 prevents non-physical coalescence of foam, so that coalescence does not occur under various conditions such as foam aggregation and deformation; ③ the wet boundary condition can accurately model the wettability of three-phase fluid on the solid surface, and reasonably represent the movement of the contact line on the solid surface. In addition, the demonstration example better demonstrates the internal mechanism of foam injection "plugging profile control" to improve sweep area and enhance oil recovery, and for the first time reveals the migration law of three-phase fluid in the foam displacement process from the perspective of numerical simulation. In summary, the method described in the present application can effectively simulate the foam injection and oil displacement process in a high water-cut reservoir.

[0161] The present application combines a variety of independently developed algorithms and models together, and can simulate non-miscible three-phase foam flow under full-range interfacial tension conditions, and embodies many characteristics: mass conservation of each phase, small pseudo-velocity, prevention of non-physical coalescence of foam, correct description of fluid wettability and near-wall solid-liquid flow characteristics, etc. Based on these characteristics, the model can reasonably reveal the migration law of three-phase foam flow in a non-uniform porous medium, and the internal mechanism of foam enhanced oil displacement and plugging profile control.

[0162] The three-phase foam flow injection simulation method for enhanced oil displacement in a high water-cut reservoir provided by the present application 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 simulations and accurately reproduce the behavior of foam in porous media; 3. In the diffusion interface framework, the wetting boundary condition is imposed by assigning a virtual fluid volume fraction to the solid boundary grid to match the required contact angle size. This not only avoids stress singularities caused by the conflict between the contact line movement and the no-slip boundary condition, but also allows for flexible adjustment of 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-cut reservoirs, comprising:

[0165] A model construction module is used to construct an immiscible three-phase flow computational model within the lattice Boltzmann framework and introduce a continuous surface force model and a recoloring algorithm into the immiscible three-phase flow computational model. 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, enabling the model to capture changes in the interfaces of droplets and bubbles in porous media. The recoloring algorithm is used to separate the phases of fluid along the interface normal while maintaining the conservation of mass of each phase, thereby generating a recognizable 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. 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 recoloring algorithm to ensure that the different fluids are separated along the interface normal and maintain the conservation of mass of each phase fluid. The contact angle geometric formula is used to apply wetting boundary conditions 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. Through time-marching simulation, multiple sets of data characterizing the three-phase flow state during the foam injection flooding process are obtained.

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

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

[0170] Further, the present application also provides a non-transitory computer readable storage medium containing instructions, and the storage medium stores a computer program. For example, the storage medium containing instructions can be executed by the processor of the computer device to complete the above method. For example, the non-transitory 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, etc. When the computer program is executed by the processor, the steps in the method embodiment of the three-phase foam flow injection simulation method for enhanced oil displacement in high water cut reservoirs can be implemented. The specific implementation method can be referred to the method embodiment, which will not be repeated here.

[0171] Those skilled in the art should understand that the embodiments of the present application can provide a method, a system or a computer program product. Therefore, the present application can take the form of a complete hardware embodiment, a complete software embodiment or an embodiment combining software and hardware aspects. Moreover, the present application can 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 application is described with reference to flowcharts and / or block diagrams of the method, device (system) and computer program product according to the embodiments of the present application. It should be understood that each flow and / or block in the flowcharts and / or block diagrams can be implemented by computer program instructions, and the combination of the flows and / or blocks in the flowcharts and / or block diagrams. These computer program instructions can be provided to the 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 produce a device that implements the functions specified in the flowcharts and / or block diagrams. Figure 1 The device that implements the functions specified in one flow or multiple flows and / or blocks. Figure 1 The device that implements the functions specified in one flow or multiple flows and / or blocks.

[0173] These computer program instructions can also be stored in a computer readable memory that can guide the computer or other programmable data processing device to work in a specific way, so that the instructions stored in the computer readable memory produce a product including instruction devices, which implement the functions specified in the flowcharts and / or block diagrams. Figure 1 The device that implements the functions specified in one flow or multiple flows and / or blocks. Figure 1the function specified in the one or more blocks.

[0174] These computer program instructions can also be loaded into computer or other programmable data processing devices, so that a series of operation steps are performed on the computer or other programmable data processing devices to generate computer-implemented processes, so that the instructions executed on the computer or other programmable data processing devices provide processes for implementing the flow Figure 1 the flow or flows and / or blocks Figure 1 the function specified in the one or more blocks.

[0175] It should be noted that the above detailed description of the specific implementation can enable those skilled in the art to have a more comprehensive understanding of the present application, but in no way limits the present application. Therefore, although the present application has been described in detail in the present specification and examples, those skilled in the art should understand that the present application can still be modified or replaced by equivalents; and all technical solutions and improvements which do not deviate from the spirit and scope of the present application are covered in the protection scope of the present application. Any reference signs in the claims should not be regarded as limiting the claims involved. Simple changes or equivalent replacements of the technical solutions which can be obviously obtained by those skilled in the art within the technical scope disclosed by the present application all belong to the protection scope of the present application.

Claims

1. A three-phase foam flow injection simulation method for enhanced oil recovery in high water content reservoirs, characterized in that: include: A computational model for an immiscible three-phase flow is constructed within the lattice Boltzmann framework, and a continuous surface force model and a recoloring algorithm are introduced into the model. 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, enabling the model to capture changes in the interfaces of droplets and bubbles in porous media. The recoloring algorithm is used to separate the phases of fluid along the interface normal while maintaining mass conservation for each phase, thereby generating a recognizable interface. When executing the collision step in 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 modify the immiscible three-phase flow calculation model. The inter-foam repulsion force is calculated by the inter-foam repulsion force model to curb the non-physical aggregation of foams. After the collision step is completed, the distribution function in moment space is first restored by inverse transformation to obtain the total density distribution function. Then, the density distribution function is assigned to the three different fluids using a recoloring algorithm to ensure that different fluids are separated along the interface normal and the mass conservation of each phase is maintained. 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. The required three-phase fluid physical properties are input and the auxiliary parameters necessary for simulation execution are provided to run the immiscible three-phase flow calculation model. Multiple sets of data characterizing the three-phase flow state during the foam injection flooding process are obtained through time-marching simulation.

2. The three-phase foam flow injection simulation method for enhanced oil recovery in high water-cut 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 in 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-cut 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 move away from each other, the specific calculation is as follows: The repulsive force generated by point Q on point P is: Where: K is the 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 at point P rep 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 foams into account in the resultant force F, the repulsive force between foams 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-cut reservoirs according to claim 1, characterized in that: The re-coloring operation of the re-coloring algorithm is performed on the density distribution function in the velocity space. Therefore, 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 recolored: Where: x represents spatial position, t is time, 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-cut reservoirs according to claim 1, characterized in that: The wetting boundary condition is applied at the fluid-solid boundary using the contact angle geometric formula, specifically including: At any local solid wall surface, the computational domain is divided into a fluid domain and a 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 follows: for a solid grid point, if any grid point among the 8 surrounding grid points is a fluid grid point, it is considered 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 and water phases 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 and l represent any two phases of the three-phase fluid o, w, and 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 at point S S and the weighted contact angle θ k After that, we can make two lines 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 of the characteristic line and the grid line in the solid domain is an invalid intersection; 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 point S1 or S2 on the horizontal or vertical grid line are all fluid grid points; if the adjacent grid point is not a fluid grid point, the characteristic line needs to be extended until the above conditions are met; Assume 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. Through flow field calculation, the fluid volume fractions at points A, B, C, and D are all known. Therefore, 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: Indicates the position coordinates of the intersection S1, Indicates the position coordinates of the intersection S2, x A Indicates 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; Indicates 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 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 judged. or Assign values ​​to the solid boundary grid 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 the oil, water, and gas three-phase fluids that match 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 construction module is used to construct an immiscible three-phase flow computational model within the lattice Boltzmann framework and introduce a continuous surface force model and a recoloring algorithm into the immiscible three-phase flow computational model. 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, enabling the model to capture changes in the interfaces of droplets and bubbles in porous media. The recoloring algorithm is used to separate the phases of fluid along the interface normal while maintaining the conservation of mass of each phase, thereby generating a recognizable 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. 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 recoloring algorithm to ensure that the different fluids are separated along the interface normal and maintain the conservation of mass of each phase fluid. The contact angle geometric formula is used to apply wetting boundary conditions 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. Through time-marching simulation, multiple sets of data characterizing the three-phase flow state during the foam injection flooding process are obtained.

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