Discrete fracture model-based carbonate rock two-phase acidification numerical simulation method

Abstract

The invention discloses a carbonate rock two-phase acidification numerical simulation method based on a discrete fracture model. The method comprises the following steps: S1, constructing a carbonate rock two-phase acidification numerical simulation physical model based on an embedded discrete fracture model; s2, dividing a matrix grid and a fracture grid, giving physical parameters of a matrix, a natural fracture and an acid solution, and loading the physical parameters into a computational domain grid; s3, constructing a two-phase seepage model, an acid liquid chemical model and a saturation model in a reservoir matrix and a crack based on the EDFM; s4, solving a saturation model by adopting a Newton-Raphson iterative algorithm, solving an acid liquid chemical model by adopting an operator splitting method, and updating physical property parameters such as porosity and permeability according to a calculation result; s5, judging whether the wormhole of the current time step breaks through or not: if not, repeating the steps S3-S5, and judging whether the wormhole of the next time step breaks through or not; if yes, a carbonate rock two-phase acidification numerical simulation result based on the embedded discrete fracture model is obtained.

Description

Technical Field

[0001] The invention relates to the technical field of acidizing, and in particular to a numerical simulation method for two-phase acidizing of carbonate rocks based on a discrete fracture model. Background Art

[0002] Acidizing is a common technique for increasing porosity and permeability in carbonate reservoirs. Improving fluidity through acidizing depends mainly on the morphology of acid-etched wormholes. Well-developed natural fractures play a key role in guiding the direction of fluid flow in the reservoir, thereby affecting the morphology of wormholes in carbonate reservoirs. In general, reservoir flow involves multiphase dynamics, and the acidizing characteristics of two-phase flow and single-phase flow are also different. Therefore, it is crucial to study the effect of complex fracture networks on wormhole morphology during two-phase acidizing of carbonates. There are a lot of studies on the effect of natural fractures on wormholes during acidizing, and experimental methods can directly observe the effect of fractures on the morphology of wormholes formed after acidizing. However, the scale of the research target is limited, and experiments that are too small (core) or too large (reservoir) are difficult to conduct. Compared with experimental methods, numerical simulation methods provide an economical and flexible method for studying acidizing processes under complex conditions. At present, scholars often use three models to simulate fractures: (1) single fracture model; (2) pseudo-fracture model; (3) discrete fracture network model. These models face difficulties in preprocessing or low computational efficiency when calculating multiphase acidizing problems in rocks with complex fracture networks, a problem that becomes more serious as the number of physical fields involved in the acidizing process increases.

[0003] In summary, in the two-phase acidizing process of carbonate rocks with complex fracture networks, cases of overall coupling of multiple physical fields are rare. Therefore, the present invention provides a numerical simulation method for two-phase acidizing of carbonate rocks based on a discrete fracture model, which can accurately simulate the two-phase acidizing process of carbonate rocks and provide a new method for two-phase acidizing of complex fracture reservoirs, which has important application significance for improving the simulation of two-phase acidizing of complex fracture reservoirs with multi-physical field coupling. Summary of the invention

[0004] In view of the above problems, the present invention aims to provide a numerical simulation method for two-phase acidizing of carbonate rocks based on a discrete fracture model.

[0005] The technical solution of the present invention is as follows:

[0006] The numerical simulation method of carbonate rock two-phase acidizing based on discrete fracture model includes the following steps:

[0007] S1: Construction of a physical model for numerical simulation of two-phase acidizing of carbonate rocks based on the embedded discrete fracture model (EDFM);

[0008] S2: Divide the matrix grid and fracture grid, give the physical parameters of the matrix, natural fractures and acid, and load the physical parameters into the computational domain grid;

[0009] S3: Based on EDFM, two-phase flow model, acid chemical model and saturation model are constructed in reservoir matrix and fractures;

[0010] S4: Use the Newton-Raphson iterative algorithm to solve the saturation model, the operator splitting method to solve the acid chemical model, and update the porosity, permeability and other physical parameters based on the calculation results;

[0011] S5: Determine whether the wormhole in the current time step has been broken through:

[0012] If there is no breakthrough, steps S3-S5 are repeated to determine whether the wormhole in the next time step has been broken through, thereby determining the effect of the fracture on the morphology of the wormhole formed after acidification;

[0013] If a breakthrough is achieved, the numerical simulation results of carbonate two-phase acidizing based on the embedded discrete fracture model are obtained.

[0014] Preferably, in step S2, the physical parameters include the reservoir matrix porosity φ m 、Fracture porosity φ f , matrix permeability K m , crack permeability K f and acid solubility α c .

[0015] Preferably, in step S3, the two-phase seepage model is constructed based on EDFM, and the crossflow coefficient needs to be solved, including the crossflow coefficient between the matrix and the fracture and between the fractures. It is a key parameter affecting the fluid crossflow velocity between the media and is related to the intersection relationship between the fracture and the matrix grid. The calculation formula is as follows:

[0016]

[0017] Where: CI is the cross-flow coefficient between matrix and crack, dimensionless; Tf i,j is the overall cross-flow coefficient between fracture unit i and fracture unit j, dimensionless; L i,k is the length of crack element k within matrix element i, m; <d> i,k is the average distance between matrix unit i and crack unit k, m; α i , α j are the crossflow coefficient of fracture unit i and the crossflow coefficient of fracture unit j, dimensionless; FA i represents the aperture of fracture unit i, m; λ t,i is the total fluidity of fracture unit i, dimensionless; dx f is the numerical discretization spacing of fracture element i, m.

[0018] Preferably, in step S3, the two-phase seepage model constructed in the reservoir matrix and in the fracture based on EDFM includes two core equations: a multiphase pressure equation and a multiphase velocity control equation;

[0019] The multiphase pressure equation includes a multiphase matrix pressure control equation and a multiphase fracture pressure control equation;

[0020] The multiphase matrix pressure governing equation:

[0021]

[0022] The multiphase fracture pressure control equation:

[0023]

[0024] Where: K m , K f are matrix permeability scalar and fracture permeability scalar, m 2 ; t is the total mobility, 1 / (Pa·s); P m , P f are the pressure of matrix and crack, Pa; Q m , Q f are the source phases of matrix and fracture, kg / (sm 2 );ψ mf , fm , ff are the cross-flow rates between matrix and crack, crack and matrix, and crack and crack, respectively, m 3 / s; Represents the divergence operator. The subscript m of each physical quantity represents that the physical quantity belongs to the matrix, and f represents that it belongs to the crack.

[0025] Among them, the cross-flow equation is calculated by the following formula:

[0026]

[0027] The cross-flow equation between fracture and matrix is:

[0028]

[0029] The cross-flow equation between cracks is:

[0030] ψ ff =Tf i,j (P i -P j ) (8)

[0031] Where: K represents permeability, m 2 ; μ represents viscosity, Pa·s; P i , P j Respectively represent the pressure of fracture unit i and the pressure of fracture unit j, Pa.

[0032] Further, according to Darcy's law, the multiphase speed control equation can be derived as follows:

[0033]

[0034] In the formula: U, U w , U o Respectively represent the total seepage velocity, water phase seepage velocity and oil phase seepage velocity, m / s; λ w , o Represent the relative fluidity of water phase and oil phase respectively, dimensionless.

[0035] Preferably, in step S3, the chemical model of the matrix and the acid solution in the fracture is constructed, including a two-phase acidification chemical reaction equation and a chemical reaction concentration equation;

[0036] The two-phase acidification chemical reaction equation is:

[0037]

[0038] where r s Represents the reaction rate on the matrix interface. The reaction between HCl and CaCO3 is approximately an irreversible first-order reaction, so the above formula can also be defined as:

[0039] r s =k s C s (13)

[0040] Where: C s represents the reaction concentration at the interface, kmol / m 3 ;k s Indicates the reaction rate, m / s. k s The calculation is as follows:

[0041]

[0042] Where: C f is defined as the pore center concentration, which is similar to C s Equal, kmol / m 3 ; R represents the universal gas constant, 8.314 J / (mol·K); T represents temperature, K.

[0043] Further, the chemical reaction concentration equation includes a chemical reaction concentration equation within the matrix and a chemical reaction concentration equation within the crack;

[0044] The chemical reaction concentration equation in the matrix is:

[0045]

[0046] The chemical reaction concentration equation in the crack is:

[0047]

[0048] Where: φ m ,φ f are the porosity of matrix and fracture, dimensionless; t is time, s; S m,w , S f,w are the saturation of water phase in the matrix and fracture, dimensionless; C f,m , C f,f Respectively represent the concentration of acid in the matrix and the center of the crack, kmol / m 3 ; U m,w , U f,w are the seepage velocities of water phase in the matrix and fractures, m / s; are the effective diffusion coefficients of the acid in the matrix and cracks, respectively, m 2 / s;a v,m 、a v,f are the pore specific surface areas of the matrix and cracks, m 2 / m 3 ; R(C s,m )、R(C s,f ) are the acid consumption rates of the matrix and fractures, respectively, which represent the amount of material flowing to a unit rock area per unit time, kmol / (s·m 2 ).

[0049] The acid consumption rate calculation formula is:

[0050]

[0051] Where: R(C s ) is the chemical reaction term, the unit is kmol / (sm 2 ), k c represents the local mass transfer coefficient of the matrix unit, which depends on the local pore structure, reaction rate and local velocity of the fluid, m / s; D e is the effective diffusion coefficient of the acid, m 2 / s;r p is the pore radius, m; Sh is the Sherwood number, dimensionless; D m is the molecular diffusion coefficient, m 2 / s; Sh ∞ is the asymptotic Sherwood number of the pores, which depends on the porous structure and is dimensionless; d h is the hydraulic diameter of the tube-filled pore, m; x is the length of the tube-filled pore, m; Sc is the Schmidt number, dimensionless; m is the ratio of the pore diameter to the pore length, dimensionless; Re p is the pore Reynolds number, dimensionless.

[0052] Preferably, in step S3, constructing the saturation model in the reservoir matrix and in the fracture mainly includes a two-phase saturation equation in the matrix and a two-phase saturation equation in the fracture;

[0053] The two-phase saturation equation in the matrix is:

[0054]

[0055] The two-phase saturation equation in the fracture is:

[0056]

[0057] Where: S o,m , S w,m are the oil phase saturation and water phase saturation of the matrix, dimensionless; f o 、f w are the coefficients of oil phase and water phase seepage velocity, dimensionless; U m , U f is the flow velocity of the matrix and the flow velocity of the fracture, m / s; Q m,o , Q m,w Source phase of oil phase and water phase in the matrix, kg / (sm 2 );S o,f , S w,f are the oil phase saturation and water phase saturation of the fracture, respectively, dimensionless; Q f,o , Q f,w The source phases of oil phase and water phase in the fracture, m 3 / s.

[0058] Preferably, in step S4, the Newton-Raphson iterative algorithm is used to solve the saturation model, the operator splitting method is used to solve the acid chemical model, and the physical parameters such as porosity and permeability are updated according to the calculation results as shown below;

[0059] The Newton-Raphson iterative algorithm is:

[0060]

[0061] H(S)=∑(fU) f ·S f (27)

[0062]

[0063] J(S t )ΔS t+Δt =-F(S t ),S t+Δt =S t +ΔS t+Δt (29)

[0064]

[0065] Where: F(S), H(S), J(S) are all matrices related to saturation, dimensionless. φ is porosity, dimensionless; V p is the flow rate, m 3 / s; Q is the flow rate, m 3 / s;S f is the mesh surface normal vector, dimensionless. p It is the flow rate during oil and gas production. From the model perspective, it is generally the flow rate of production wells and injection wells. In this article, Q = 0.

[0066] The operator splitting method is:

[0067]

[0068] Where: When the operator splitting method processes the source phase, it is split into two parts, S u represents a part of it; |d| represents the distance between the centers of two adjacent grids, m. C f It is defined as the pore center concentration, p represents the physical quantity at the grid center, and f represents the interpolation of the physical quantity at the grid center onto the surface.

[0069] The method for updating the porosity and permeability calculation according to the calculation results is:

[0070]

[0071] Where: a v is the pore specific surface area, m 2 / m 3 ; α c is the acid solubility, kg / kmol; ρ s is the rock density, kg / m 3 ; K init , K new are the initial and post-dissolution permeabilities, m 2 ; φ init ,φ new are the initial porosity and the porosity after dissolution, dimensionless; β is an empirical parameter, dimensionless; r p,init 、r p,new are the initial and post-dissolution pore radii, m; a v,init 、a v,new are the initial and post-dissolution pore specific surface areas, m 2 / m 3 ; |U| represents the modulus of seepage velocity, m / s; Represent the effective diffusion coefficients in the x and y directions, respectively, m 2 / s; α os , L , T Both are 0.25, which are dimensionless constants related to the pore structure.

[0072] Furthermore, the breakthrough condition of whether the wormhole is broken is defined as the acid injection pressure being reduced to 1% of the initial acid injection pressure.

[0073] Furthermore, the physical parameters include the reservoir matrix heterogeneous porosity φ m =[0.05, 0.04], fracture porosity φ f =0.99, matrix heterogeneity permeability K m =[1e -14 , 1e -16 ]m 2 , crack permeability K f =1e -10 m 2 and acid solubility α c =50kg / kmol.

[0074] The beneficial effects of the present invention are:

[0075] The carbonate two-phase acidizing model constructed in the present invention based on the embedded discrete fracture model method takes into account the characteristics of two-phase seepage, strong reservoir heterogeneity, and complex fracture development, and realizes accurate simulation of the carbonate two-phase acidizing process. It provides a new method for considering the two-phase acidizing of complex fractured reservoirs, which can further realistically simulate the process of carbonate two-phase acidizing, and has important application significance for improving the two-phase acidizing simulation of complex fractured reservoirs coupled with multiple physical fields. BRIEF DESCRIPTION OF THE DRAWINGS

[0076] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, 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 paying creative labor.

[0077] Figure 1 is a flow chart of the overall method of an embodiment of the present invention;

[0078] Figure 2 A physical model diagram used in calculations in an embodiment of the present invention;

[0079] Figure 3 This is a porosity distribution cloud diagram of two-phase acidizing without natural fractures in an embodiment of the present invention;

[0080] Figure 4 This is a porosity distribution cloud diagram of two-phase acidizing of natural fractures developed in an embodiment of the present invention;

[0081] Figure 5 This is a porosity distribution cloud diagram of single-phase acidizing without natural fractures in an embodiment of the present invention;

[0082] Figure 6 This is a porosity distribution cloud diagram of single-phase acidizing of natural fractures developed in an embodiment of the present invention. DETAILED DESCRIPTION

[0083] The present invention is further described below in conjunction with the accompanying drawings and embodiments. It should be noted that, in the absence of conflict, the embodiments in this application and the technical features in the embodiments can be combined with each other. It should be pointed out that, unless otherwise specified, all technical and scientific terms used in this application have the same meanings as those generally understood by ordinary technicians in the technical field to which this application belongs. The words "including" or "comprising" and the like used in the disclosure of the present invention mean that the elements or objects appearing before the word cover the elements or objects listed after the word and their equivalents, without excluding other elements or objects.

[0084] like Figure 1 As shown, the present invention proposes a numerical simulation method for two-phase acidizing of carbonate rock based on a discrete fracture model, comprising the following steps:

[0085] S1: A physical model for numerical simulation of two-phase acidizing of carbonate rocks is constructed based on the embedded discrete fracture model (EDFM).

[0086] S2: Divide the matrix grid and fracture grid, give the physical parameters of the matrix, natural fractures and acid, and load the physical parameters into the computational domain grid;

[0087] In a specific embodiment, the physical parameters include the reservoir matrix heterogeneous porosity φ m =[0.05, 0.04], fracture porosity φ f =0.99, matrix heterogeneity permeability K m =[1e -14 , 1e -16 ]m 2 , crack permeability K f =1e - 10 m 2 and acid solubility α c =50kg / kmol.

[0088] S3: Based on EDFM, a two-phase seepage model, acid fluid chemistry model, and saturation model are constructed in the reservoir matrix and fractures.

[0089] Furthermore, the method further comprises the steps of:

[0090] S4: The Newton-Raphson iterative algorithm is used to solve the saturation model, the operator splitting method is used to solve the acid chemical model, and the physical properties such as porosity and permeability are updated according to the calculation results.

[0091] Furthermore, the method further comprises the steps of:

[0092] S5: Determine whether the wormhole in the current time step has been broken through:

[0093] If not, repeat steps S3-S5 to determine whether the wormhole is broken through in the next time step;

[0094] If a breakthrough is achieved, the numerical simulation results of two-phase acidizing of carbonate rocks based on the embedded discrete fracture model are obtained to determine the influence of fractures on the morphology of wormholes formed after acidizing.

[0095] The breakthrough condition of whether the wormhole is broken is defined as the acid injection pressure being reduced to 1% of the initial acid injection pressure.

[0096] In a specific embodiment, constructing a two-phase seepage model based on EDFM requires solving the crossflow coefficient, including the crossflow coefficient between the matrix and the fracture and between the fractures. It is a key parameter affecting the fluid crossflow velocity between the media and is related to the intersection relationship between the fracture and the matrix grid. The calculation formula is as follows:

[0097]

[0098] Where: CI is the cross-flow coefficient between matrix and crack, dimensionless; Tf i,j is the overall cross-flow coefficient between fracture unit i and fracture unit j, dimensionless; L i,k is the length of crack element k within matrix element i, m; <d> i,k is the average distance between matrix unit i and crack unit k, m; α i , α j are the crossflow coefficients of fracture unit i and fracture unit j, dimensionless, α j It is also calculated using formula (2), where i is replaced by j; FA i represents the aperture of fracture unit i, m; λ t,i is the total fluidity of fracture unit i, dimensionless; dx f is the numerical discretization spacing of fracture element i, m.

[0099] Table 1 shows the sources of the parameters in formulas (1) to (3)

[0100] constant Measured Calculated CI √ <![CDATA[Tf i,j ]]> √ <![CDATA[L i,k ]]> √ <![CDATA[ <d> i,k ]]> < / d> √ <![CDATA[α i ,a j ]]> √ <![CDATA[FA i ]]> √ <![CDATA[λ t,i ]]> √ <![CDATA[dx f ]]> √

[0101] In a specific embodiment, in step S3, the two-phase seepage model constructed in the reservoir matrix and in the fracture based on EDFM includes two core equations: a multiphase pressure equation and a multiphase velocity control equation;

[0102] The multiphase pressure equation includes a multiphase matrix pressure control equation and a multiphase fracture pressure control equation;

[0103] The multiphase matrix pressure governing equation:

[0104]

[0105] The multiphase fracture pressure control equation:

[0106]

[0107] Where: K m , K f are matrix permeability scalar and fracture permeability scalar, m 2 ; t is the total mobility, 1 / (Pa·s); P m , P f are the pressure of matrix and crack, Pa; Q m , Q f are the source phases of matrix and fracture, kg / (sm 2 );ψ mf , fm , ff are the cross-flow rates between matrix and crack, crack and matrix, and crack and crack, respectively, m 3 / s; Represents the divergence operator.

[0108] Table 2 shows the sources of various parameters in formulas (4) and (5)

[0109]

[0110] Among them, the cross-flow equation is calculated by the following formula:

[0111]

[0112] The cross-flow equation between fracture and matrix is:

[0113]

[0114] Fracture-fracture crossflow equation:

[0115] ψ ff =Tf i,j (P i -P j ) (8)

[0116] Where: K represents permeability, m 2 ; μ represents viscosity, Pa·s; P i , P j Respectively represent the pressure of fracture unit i and the pressure of fracture unit j, Pa.

[0117] Table 3 shows the sources of various parameters in formulas (6) to (8)

[0118] constant Measured Calculated k √ μ √ <![CDATA[P i 、P j ]]> √

[0119] According to Darcy's law, the multiphase velocity control equation can be derived as:

[0120]

[0121] In the formula: U, U w , U o Represent the total seepage velocity, water phase seepage velocity and oil phase seepage velocity, m / s; is the gradient operator, which means the gradient of the calculated pressure P; w , o They represent the relative fluidity of the water phase and the oil phase respectively. They are dimensionless. In the subscripts of each physical quantity, w represents the water phase and o represents the oil phase.

[0122] Table 4 shows the sources of various parameters in formulas (9) to (11)

[0123]

[0124] In a specific embodiment, in step S3, the chemical model of the matrix and the acid fluid in the fracture is constructed, including a two-phase acidification chemical reaction equation and a chemical reaction concentration equation;

[0125] The two-phase acidification chemical reaction equation is:

[0126]

[0127] where r s Represents the reaction rate on the matrix interface. The reaction between HCl and CaCO3 is approximately an irreversible first-order reaction, so the above formula can also be defined as:

[0128] r s =k s C s (13)

[0129] Where: C s represents the reaction concentration at the interface, kmol / m 3 ;k s Indicates the reaction rate, m / s. k s The calculation is as follows:

[0130]

[0131] Where: C f is defined as the pore center concentration, which is similar to C s Equal, kmol / m 3 ; R represents the universal gas constant, 8.314 J / (mol·K); T represents temperature, K.

[0132] Table 5 shows the sources of various parameters in formulas (12) to (14)

[0133] Variables in the numerical simulation process Measured Calculated <![CDATA[C s ]]> √ <![CDATA[k s ]]> √ T √

[0134] The chemical reaction concentration equation includes a chemical reaction concentration equation within the matrix and a chemical reaction concentration equation within the crack;

[0135] The chemical reaction concentration equation in the matrix is:

[0136]

[0137] The chemical reaction concentration equation in the crack is:

[0138]

[0139] Where: φ m ,φ f are the porosity of matrix and fracture, dimensionless; t is time, s; S m,w , S f,w are the saturation of water phase in the matrix and fracture, dimensionless; C f,m , C f,f Respectively represent the concentration of acid in the matrix and the center of the crack, kmol / m 3 ; U m,w , U f,w are the seepage velocities of water phase in the matrix and fractures, m / s; are the effective diffusion coefficients of the acid in the matrix and cracks, respectively, m 2 / s;a v,m 、a v,f are the pore specific surface areas of the matrix and cracks, m 2 / m 3 ; R(C s,m )、R(C s,f ) are the acid consumption rates of the matrix and fractures, respectively, which represent the amount of material flowing to a unit rock area per unit time, kmol / (s·m 2 ).

[0140] The acid consumption rate calculation formula is:

[0141]

[0142]

[0143] Where: R(C s ) is the chemical reaction term, the unit is kmol / (sm 2 ), calculate R(C s,m )、R(C s,f ) also uses formula (17); k c represents the local mass transfer coefficient of the matrix unit, which depends on the local pore structure, reaction rate and local velocity of the fluid, m / s; D e is the effective diffusion coefficient of the acid, m 2 / s;r p is the pore radius, m; Sh is the Sherwood number, dimensionless; D m is the molecular diffusion coefficient, m 2 / s; Sh ∞ is the asymptotic Sherwood number of the pores, which depends on the porous structure and is dimensionless; d h is the hydraulic diameter of the tube-filled pore, m; x is the length of the tube-filled pore, m; Sc is the Schmidt number, dimensionless; m is the ratio of the pore diameter to the pore length, dimensionless; Re p is the pore Reynolds number, dimensionless.

[0144] Table 6 shows the sources of various parameters in formulas (17) to (19)

[0145] constant Measured Calculated <![CDATA[k c ]]> √ <![CDATA[D e ]]> √ <![CDATA[r p ]]> √ <![CDATA[D m ]]> √ <![CDATA[d h ]]> √ Sh √ x √ m √ Sc √ <![CDATA[Re p ]]> √

[0146] In a specific embodiment, in step S3, constructing the saturation model in the reservoir matrix and in the fracture mainly includes a two-phase saturation equation in the matrix and a two-phase saturation equation in the fracture;

[0147] The two-phase saturation equation in the matrix is:

[0148]

[0149] The two-phase saturation equation in the fracture is:

[0150]

[0151]

[0152] Where: S o,m , S w,m are the oil phase saturation and water phase saturation of the matrix, dimensionless; f o 、f w are the coefficients of oil phase and water phase seepage velocity, dimensionless; U m , U f is the flow velocity of the matrix and the flow velocity of the fracture, m / s; Q m,o , Q m,w Source phase of oil phase and water phase in the matrix, kg / (sm 2 );S o,f , S w,f are the oil phase saturation and water phase saturation of the fracture, respectively, dimensionless; Q f,o , Q f,w Source phases of oil phase and water phase in fractures, kg / (sm 2 ).

[0153] S4: The Newton-Raphson iterative algorithm is used to solve the saturation model, the operator splitting method is used to solve the acid chemical model, and the physical properties such as porosity and permeability are updated according to the calculation results.

[0154] In a specific embodiment, in step S4, the Newton-Raphson iterative algorithm is used to solve the saturation model, the operator splitting method is used to solve the acid chemical model, and the physical parameters such as porosity and permeability are updated according to the calculation results as shown below:

[0155] The Newton-Raphson iterative algorithm is:

[0156]

[0157] H(S)=∑(fU) f ·S f (27)

[0158]

[0159] J(S t )ΔS t+Δt =-F(S t ),S t+Δt =S t +ΔS t+Δt (29)

[0160]

[0161] Where: F(S), H(S), J(S) are all matrices related to saturation, dimensionless. φ is porosity, dimensionless; V p is the flow rate, m 3 / s; Q is the flow rate, m 3 / s; f represents the interpolation from the grid center to the grid surface, not a physical quantity; S f is the mesh surface normal vector, dimensionless. p It is the flow rate during oil and gas production. From the model perspective, it is generally the flow rate of production wells and injection wells. In this article, Q = 0.

[0162] The operator splitting method is:

[0163]

[0164] Where: When the operator splitting method processes the source phase, it is split into two parts, S u represents a part of it; |d| represents the distance between the centers of two adjacent grids, m. The above formula is used to solve the acid concentration, which affects R(C s ), and then affect the porosity update. C f It is defined as the pore center concentration, p represents the physical quantity at the grid center, and f represents the interpolation of the physical quantity at the grid center onto the surface.

[0165] The method for updating the porosity and permeability calculation according to the calculation results is:

[0166]

[0167] Where: a v is the pore specific surface area, m 2 / m 3 ; α c is the acid solubility, kg / kmol; ρ s is the rock density, kg / m 3 ; K init , K new are the initial and post-dissolution permeabilities, m 2 ; φ init ,φ new are the initial porosity and the porosity after dissolution, dimensionless; β is an empirical parameter, dimensionless; r p,init 、r p,new are the initial and post-dissolution pore radii, m; a v,init 、a v,new are the initial and post-dissolution pore specific surface areas, m 2 / m 3 ; |U| represents the modulus of seepage velocity, m / s; Represent the effective diffusion coefficients in the x and y directions, respectively, m 2 / s; α os , L , T Both are 0.25, which are dimensionless constants related to the pore structure.

[0168] S5: Determine whether the wormhole in the current time step has been broken through:

[0169] If not, repeat steps S3-S5 to determine whether the wormhole is broken through in the next time step;

[0170] If a breakthrough is achieved, the numerical simulation results of two-phase acidizing of carbonate rocks based on the embedded discrete fracture model are obtained to determine the influence of fractures on the morphology of wormholes formed after acidizing.

[0171] In a specific embodiment, the breakthrough condition for determining whether the wormhole is broken is defined as the acid injection pressure being reduced to 1% of the initial acid injection pressure.

[0172] In this embodiment, if Figure 2 The physical model used in this embodiment is shown in Figure 1. The two-phase acidizing process problem is studied as an example. The problem area is a core scale with a size of 0.3m×0.1m (length×height). The rock contains natural fractures and is saturated with liquid (water (acid) or oil). The acid is injected from the left side of the rock to corrode the channels inside the rock to increase its porosity and permeability. The matrix area of ​​the physical model is meshed using a structured uniform grid. Complex fractures are randomly displayed and meshed using a structured uniform grid ( Figure 2 It is a schematic diagram, not the actual number of grids divided in the embodiment), and the natural fractures are located inside the matrix.

[0173] like Figure 3-Figure 6 As shown in the figure, it is a porosity distribution cloud map of carbonate rock acidification process under different conditions. It can be seen that when there is no natural fracture, the wormhole distribution characteristics are close to a straight line; when there is a natural fracture, the acid fluid preferentially enters the natural fracture and forms a dissolution zone at the entrance of the natural fracture; when there is also a natural fracture, the wormhole width of the single-phase flow is significantly larger than the wormhole width of the two-phase flow, which indicates that the oil phase in the reservoir will have a more obvious hindering effect on the acid fluid. It can be seen that this method can accurately simulate the oil-water two-phase flow acidification process of fractured carbonate.

[0174] The above description is only a preferred embodiment of the present invention and does not constitute any form of limitation to the present invention. Although the present invention has been disclosed as a preferred embodiment as above, it is not intended to limit the present invention. Any technician familiar with the profession can make some changes or modifications to equivalent embodiments of equivalent changes using the technical contents disclosed above without departing from the scope of the technical solution of the present invention. However, any simple modification, equivalent change and modification made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solution of the present invention still falls within the scope of the technical solution of the present invention.< / d> < / d>

Claims

1. A numerical simulation method for two-phase acidizing of carbonate rocks based on a discrete fracture model, characterized in that: The following steps are involved: S1: Construction of a physical model for numerical simulation of two-phase acidizing of carbonate rocks based on the embedded discrete fracture model (EDFM); S2: Divide the matrix grid and fracture grid, give the physical parameters of the matrix, natural fractures and acid, and load the physical parameters into the computational domain grid; S3: Based on EDFM, two-phase flow model, acid chemical model and saturation model are constructed in reservoir matrix and fractures; S4: Use the Newton-Raphson iterative algorithm to solve the saturation model, the operator splitting method to solve the acid chemical model, and update the porosity, permeability and other physical parameters based on the calculation results; S5: Determine whether the wormhole in the current time step has been broken through: If not, repeat steps S3-S5 to determine whether the wormhole is broken through in the next time step; If a breakthrough is achieved, the numerical simulation results of two-phase acidizing of carbonate rocks based on the embedded discrete fracture model are obtained to determine the influence of fractures on the morphology of wormholes formed after acidizing.

2. The method for numerical simulation of carbonate rock two-phase acidizing based on discrete fracture model according to claim 1, characterized in that: In step S2, the physical parameters include the reservoir matrix porosity φ m 、Fracture porosity φ f , matrix permeability K m , crack permeability K f and acid solubility α c .

3. The method for numerical simulation of carbonate rock two-phase acidizing based on discrete fracture model according to claim 1, characterized in that: In step S3, the construction of a two-phase seepage model based on EDFM requires solving the crossflow coefficient, including the crossflow coefficient between the matrix and the fracture and between the fractures. It is a key parameter affecting the fluid crossflow velocity between the media and is related to the intersection relationship between the fracture and the matrix grid. The calculation formula is as follows: Where: CI is the cross-flow coefficient between matrix and crack, dimensionless; Tf i,j is the overall cross-flow coefficient between fracture unit i and fracture unit j, dimensionless; L i,k is the length of crack element kj within matrix element i, m; <d> i,k is the average distance between matrix unit i and crack unit k, m; α i , α j are the crossflow coefficient of fracture unit i and the crossflow coefficient of fracture unit j, dimensionless;< / d> FA i represents the aperture of fracture unit i, m; λ t,i is the total fluidity of fracture unit i, dimensionless; dx f is the numerical discretization spacing of fracture element i, m.

4. The method for numerical simulation of two-phase acidification of carbonate rock based on discrete fracture model according to claim 3, characterized in that: In step S3, the two-phase seepage model constructed in the reservoir matrix and in the fracture based on EDFM includes two core equations: a multiphase pressure equation and a multiphase velocity control equation; The multiphase pressure equation includes a multiphase matrix pressure control equation and a multiphase fracture pressure control equation; The multiphase matrix pressure governing equation: The multiphase fracture pressure control equation: Where: K m , K f are matrix permeability scalar and fracture permeability scalar, m 2 ; λ t is the total liquidity, 1 / (Pa·s); represents the divergence operator; P m , P f are the pressure of matrix and crack, Pa; Q m , Q f are the source phases of matrix and fracture, kg / (sm 2 ).

5. The method for numerical simulation of carbonate rock two-phase acidizing based on discrete fracture model according to claim 4, characterized in that: The ψ mf , mf , ff are the cross-flow rates between matrix and crack, crack and matrix, and crack and crack, respectively, m 3 / s; Among them, the cross-flow equation is calculated by the following formula: The cross-flow equation between fracture and matrix is: The cross-flow equation between cracks is: ψ ff (Tf i,j (P i -P j ) (8) Where: K represents permeability, m 2 ; μ represents viscosity, Pa·s; P i , P j Respectively represent the pressure of fracture unit i and the pressure of fracture unit j, Pa.

6. The method for numerical simulation of carbonate rock two-phase acidizing based on discrete fracture model according to claim 5, characterized in that: According to Darcy's law, the multiphase velocity control equation can be derived as: In the formula: U, U w , U o Respectively represent the total seepage velocity, water phase seepage velocity and oil phase seepage velocity, m / s; λ w , o Represent the relative fluidity of water phase and oil phase respectively, dimensionless.

7. The method for numerical simulation of carbonate rock two-phase acidizing based on discrete fracture model according to claim 6, characterized in that: In step S3, the chemical model of the matrix and the acid fluid in the fracture is constructed, including a two-phase acidification chemical reaction equation and a chemical reaction concentration equation; The two-phase acidification chemical reaction equation is: where r s Represents the reaction rate on the matrix interface; the reaction between HCl and CaCO3 is approximately an irreversible first-order reaction, so the above formula can also be defined as: r s =k s C s (13) Where: C s represents the reaction concentration at the interface, kmol / m 3 ;k s Indicates the reaction rate, m / s; k s The calculation is as follows: Where: C f is defined as the pore center concentration, which is similar to C s Equal, kmol / m 3 ; R represents the universal gas constant, 8.314 J / (mol·K); t represents temperature, K.

8. The method for numerical simulation of carbonate rock two-phase acidizing based on discrete fracture model according to claim 7, characterized in that: The chemical reaction concentration equation includes a chemical reaction concentration equation within the matrix and a chemical reaction concentration equation within the crack; The chemical reaction concentration equation in the matrix is: The chemical reaction concentration equation in the crack is: Where: φ m ,φ f are the porosity of matrix and fracture, respectively, dimensionless; t is time, s; S m,w , S f,w are the saturation of water phase in the matrix and fracture, dimensionless; C f,m , C f,f Respectively represent the concentration of acid in the matrix and the center of the crack, kmol / m 3 ; U m,w , U f,w are the seepage velocities of water phase in the matrix and fractures, m / s; are the effective diffusion coefficients of the acid in the matrix and cracks, respectively, m 2 / s;a v,m 、a v,f are the pore specific surface areas of the matrix and cracks, m 2 / m 3 ; R(C s,m )、R(C s,f ) are the acid consumption rates of the matrix and fractures, respectively, which represent the amount of material flowing to a unit rock area per unit time, kmol / (s·m 2 ); The acid consumption rate calculation formula is: Where: k c represents the local mass transfer coefficient of the matrix unit, which depends on the local pore structure, reaction rate and local velocity of the fluid, m / s; D e is the effective diffusion coefficient of the acid, m 2 / s;r p is the pore radius, m; Sh is the Sherwood number, dimensionless; D m is the molecular diffusion coefficient, m 2 / s; Sh ∞ is the asymptotic Sherwood number of the pores, which depends on the porous structure and is dimensionless; d h is the hydraulic diameter of the tube-filled pore, m; x is the length of the tube-filled pore, m; Sc is the Schmidt number, dimensionless; m is the ratio of the pore diameter to the pore length, dimensionless; Re p is the pore Reynolds number, dimensionless.

9. The method for numerical simulation of carbonate rock two-phase acidizing based on discrete fracture model according to claim 8, characterized in that: In step S3, constructing the saturation model in the reservoir matrix and in the fracture mainly includes the two-phase saturation equation in the matrix and the two-phase saturation equation in the fracture; The two-phase saturation equation in the matrix is: The two-phase saturation equation in the fracture is: Where: S o,m , S w,m They are the oil phase saturation of the matrix and the water phase saturation of the matrix, dimensionless; f o 、f w are the coefficients of oil phase and water phase seepage velocity, dimensionless; U m , U f is the flow velocity of the matrix and the flow velocity of the fracture, m / s; Q m,o , Q m,w Source phase of oil phase and water phase in the matrix, kg / (sm 2 );S o,f , S w,f are the oil phase saturation and water phase saturation of the fracture, respectively, dimensionless; Q f,o , Q f,w Source phases of oil phase and water phase in fractures, kg / (sm 2 ).

10. The method for numerical simulation of carbonate rock two-phase acidizing based on discrete fracture model according to claim 9, characterized in that: In step S4, the Newton-Raphson iterative algorithm is used to solve the saturation model, the operator splitting method is used to solve the acid chemical model, and the physical parameters such as porosity and permeability are updated according to the calculation results as shown below: The Newton-Raphson iterative algorithm is: H(S)=∑(fU) f ·S f (27) J(S t )ΔS t+Δt =-F(S t ),S t+Δt =S t +ΔS t+Δt (29) Where: F(S), H(S), J(S) are all matrices related to saturation, dimensionless; φ is porosity, dimensionless; V p is the flow rate, m 3 / s; Q is the flow rate, m 3 / s;S f is the mesh surface normal vector, dimensionless; The operator splitting method is: Where: When the operator splitting method processes the source phase, it is split into two parts, S u represents a part of it; |d| represents the distance between the centers of two adjacent grids, m; The porosity and permeability calculation methods are updated according to the calculation results: Where: a v is the pore surface area, m 2 / m 3 ; α c is the acid solubility, kg / kmol; ρ s is the rock density, kg / m 3 ; K init , K new are the initial and post-dissolution permeabilities, m 2 ; φ init ,φ new are the initial porosity and the porosity after dissolution, dimensionless; β is an empirical parameter, dimensionless; r p,init 、r p,new are the initial and post-dissolution pore radii, m; a v,init 、a v,new are the initial and post-dissolution pore specific surface areas, m 2 / m 3 ; |U| represents the modulus of the seepage velocity, m / s; Represent the effective diffusion coefficients in the x and y directions, respectively, and m 2 / s; α os , L , T Both are 0.25, which are dimensionless constants related to the pore structure.

Images