A method and system for numerical simulation of radial flow of fluid in porous media based on an acidizing model implemented in a DBF framework

Abstract

The application relates to a numerical simulation method and system for radial flow of fluid in porous media based on an acidification model of a DBF framework, which comprises the following steps: step 1, constructing an acidification model for describing fluid flow, solute reaction transport and rock property change; step 2, carrying out grid division on a simulation time and a circular ring type solving region; then defining different variables of the acidification model at different positions of the grid cells; step 3, carrying out numerical discretization on the acidification model constructed in step 1 to form a large linear equation group; finally, solving in combination with set parameters to realize numerical simulation of acid etching wormholes. The application provides an interlaced grid finite difference method for numerical simulation of radial flow of fluid in porous media based on the acidification model of the DBF framework, numerical discretization and efficient solving are carried out on the acidification model, relatively accurate numerical simulation of the acid etching wormhole process is realized, and then an important guiding role is played on oil and gas exploitation.

Description

TECHNICAL FIELD

[0001] The application relates to a method and system for numerical simulation of radial flow of fluid in porous media based on an acidification model of a DBF framework, and belongs to the technical field of numerical simulation of acid-etched wormholes in carbonate matrix acidification. BACKGROUND

[0002] Acidification treatment is a widely used stimulation technique in the oil industry, aiming to increase the effective connectivity between the wellbore and the carbonate reservoir, and then efficiently extract the oil in the reservoir to the surface. Matrix acidification is considered an effective and efficient acidification technique for carbonate formations, which can increase the permeability and porosity of carbonate rocks and improve the acidification efficiency and oil and gas production.

[0003] Acid-etched wormholes refer to the injection of acid into the wellbore, and after the acid flows into the carbonate reservoir, it reacts with the minerals in the formation, changes the rock structure, increases the permeability and porosity of the carbonate rock, and then generates a series of irregular holes like earthworms. As high-porosity flow channels, acid-etched wormholes can change the flow characteristics of fluid in porous media, improve acidification efficiency, and increase oil and gas production. In addition, the growth trajectory and trend of wormholes are closely related to the effect of carbonate matrix acidification, therefore, in order to achieve good acidification effect, numerical simulation of the process of acid-etched wormholes is crucial, which can guide oil and gas production to some extent.

[0004] Most of the numerical simulations of acid etching wormholes are realized by solving acidizing models based on the Darcy framework in the Cartesian coordinate system. Kou and Sun et al. [Comput. Methods Appl. Mech. Engrg., 2016, 298: 279-302] proposed a globally conservative hybrid finite element method for solving acidizing models based on the Darcy-Forchheimer framework in the Cartesian coordinate system. Li and Rui [J. Sci. Comput., 2018, 74:1115-1145] proposed a block-centered finite difference method for solving compressible acidizing models based on the Darcy framework in the Cartesian coordinate system and [J. Fluid Mech., 2019, 872: 438-471] proposed a globally conservative finite difference method for simulating Darcy-Brinkman-Forchheimer acidizing models in the Cartesian coordinate system. Yang et al. [J. Comput. Phys., 2023, 473] proposed a nonlinear complementarity simulator with CPR preconditioner for simulating acid etching wormholes based on the Darcy-Brinkman framework in the Cartesian coordinate system. In addition, Guo et al. [J. Sci. Comput., 2021, 89(1)] proposed a high-order bound-preserving finite difference method for solving incompressible acidizing models based on the Darcy framework in the Cartesian coordinate system. Xia et al. [Chinese Science Bulletin, 2024, 1-30] gave the algorithmic study of compressible acidizing models based on the Darcy framework in the two-dimensional polar coordinate system.

[0005] Acidizing models based on the Darcy framework can well describe the acid etching wormhole process when the porosity does not change significantly; however, once the porosity evolves significantly with time, the flow velocity of the fluid in the high porosity region will increase, and the acidizing model under the Darcy framework is not sufficient to characterize the above acid etching wormhole process. In order to better realize the numerical simulation of acid etching wormholes, it is necessary to establish acidizing models based on the Darcy-Brinkman-Forchheimer framework. At the same time, the actual matrix acidizing is to inject acid into a circular wellbore, and the flow of acid in this process is radial from the wellbore to the carbonate reservoir, so it is further necessary to establish a Darcy-Brinkman-Forchheimer acidizing model in the polar coordinate system to more accurately simulate the acid etching wormhole process in actual industrial applications.

[0006] The staggered grid is characterized by storing different variables at different positions of the grid, that is, for the established acidification model about pressure, velocity, concentration, and porosity, the pressure, concentration, and porosity are defined at the center of the grid cell, and the velocity is defined at the midpoint of the four edges of the grid cell. The acidification model in polar coordinates is numerically discretized by using the staggered grid finite difference method, and the original physical properties of the model, such as mass conservation and momentum conservation, can be well maintained, and the numerical solution process is simple and efficient.

[0007] The acid-etch wormhole process can effectively describe the acidification effect of the matrix, and for the radial flow of the fluid in the carbonate reservoir, a simple and efficient staggered grid finite difference method is used to numerically discretize and solve the acidification model close to the actual acidification model, which can effectively improve the numerical simulation precision and has certain guiding significance for oil and gas exploitation. SUMMARY

[0008] The first aspect of the application provides a numerical simulation method for radial flow of fluid in porous media based on the DBF framework of the acidification model; comprising:

[0009] Step 1: constructing an acidification model for describing fluid flow, solute reaction transport and rock property change;

[0010] Step 2: meshing the simulation time and the annular solving region; and defining different variables of the acidification model at different positions of the grid cell;

[0011] Step 3: numerically discretizing the acidification model constructed in step 1 to form a large linear equation group; finally, combining with the set parameters to solve, and realizing the acid-etch wormhole numerical simulation.

[0012] According to the application, the acidification model for describing fluid flow, solute reaction transport and rock property change is constructed; comprising:

[0013] A set of two-dimensional polar coordinate-based acidification models based on the DBF framework is established, as shown below:

[0014] ; (1)

[0015] ; (2)

[0016]

[0017] ; (3)

[0018] ; (4)

[0019] Wherein, equations (1) and (2) represent the fluid flow process based on the DBF framework, i.e., the acidification model used to describe the fluid flow; equation (1) is a vector equation representing the momentum conservation equation, and the right-hand side of equation (1) is: The term represents the Darcy term, used to characterize the Darcy flow phenomenon in porous media; the second term on the left side of formula (1): The Brinkman term describes the transition flow between boundaries; the last term on the left side of formula (1): The term represents the Forchheimer term, also known as the inertial term, which is used to describe the significant inertial effect of the fluid when the flow rate is high; Formula (2) is the mass conservation equation; Formula (3) is the acid concentration reaction transport equation, which is used to describe the acidification model for solute reaction transport; Formula (4) reflects the change of rock porosity over time, which is used to describe the acidification model for changes in rock properties.

[0020] in, , It is the flow radius, in units of ; It is a circular region; It is time. For the final moment, the unit is ; It is the fluid velocity vector. These are the fluid velocities. At the polar radius Direction and polar angle Components of direction, in units of ; , Indicates polar radius The boundary of direction; It is the pressure of the fluid, and the unit is . ; It is the concentration of acid, in units of ; It refers to the porosity of rocks, a dimensionless quantity; velocity. ,pressure ,concentration Porosity All four variables are unknowns; It is the density of the fluid, in units of... ; It is the viscosity of the fluid, measured in units of... ; It is a pseudo compression factor, with units of , This causes a slight change in fluid density during the solute dissolution process. This represents a positive number to ensure that the coefficient matrix is ​​invertible; is the local mass transport coefficient, with units of ; is the injection concentration, with units of ; is the acid's dissolution ability, with units of ; is the rock's density, with units of ; is the Forchheimer coefficient, dimensionless; is the rock's permeability, with units of ; where is the injection rate, is the production rate, with units of ; the positive definite matrix is the acid's diffusion coefficient in the porous medium, with units of , and are the components of in the and directions; for simplicity, we assume is a diagonal matrix, where is the molecular diffusion rate, is the identity matrix, denotes a diagonal matrix with as diagonal entries; are given functions, while is a function of the rock's porosity;

[0021] is the acid's concentration at the fluid (acid) and solid (rock) interface, with units of , according to first-order kinetic reactions, the acid's concentration at the fluid-solid interface is related to the acid's concentration in the fluid by the following relation:

[0022] ; (5)

[0023] where is the surface reaction rate constant, with units of ;

[0024] The pore-scale model that characterizes the rock's property variations is as follows:

[0025] ; (6)

[0026] ; (7)

[0027] wherein formula (6) is established by Carman-Kozeny, and is used to depict the relationship between the porosity and the permeability of the rock, represents the permeability, and are the initial porosity and the initial permeability of the rock respectively; the porosity and the permeability are calculated wherein, is the interface area (specific surface area) of the unit volume medium for the reaction, and the unit is , is the initial interface area;

[0028] The boundary and initial conditions are as follows:

[0029] (8)

[0030] wherein the two conditions in the first row: are boundary conditions, indicating that the velocity and the gradient of the concentration are both 0 at the inlet boundary and the outlet boundary; the following four are all initial conditions, , , , represent the initial distribution functions of the velocity, the pressure, the concentration and the porosity in the region respectively;

[0031] The acidification model formula (1)-(4) based on the DBF framework in the two-dimensional polar coordinates is combined with the pore-scale model formula (6)-(7) and the initial boundary condition formula (8) to construct the acidification model of the application; the flow properties of the non-Darcy seepage of the porous medium in the acidification process of the carbonate rock matrix are comprehensively considered, and compared with the prior art, the acidification model established in the application can more accurately depict the acid corrosion wormhole process;

[0032] In order to simplify the equation and the writing of the following discrete format, an auxiliary variable is defined, and combined with formula (2), (5)-(7), the acid concentration reaction transmission formula (3) is simplified and converted into the following formula (9):

[0033] (9).

[0034] According to the application, the simulation time and the annular solving region are meshed; then different variables of the acidification model are defined at different positions of the grid cells; including:

[0035] The total simulation time (the final moment) is evenly divided into parts, and the time step is and the first moment ; in order to better simulate the situation in reality, for a given simulation area , namely the annular solution area, non-uniform grid division is carried out, the area length in the direction is divided into parts, the area length in the direction is divided into parts, forming grid units, and the grid points are: and ; the midpoint of the grid unit edge in the direction and the direction, and the division step , , , are respectively:

[0036]

[0037] The staggered grid is characterized in that different variables can be stored in different positions of the grid unit; for the nonlinear strong coupling acidification model established in step 1 about pressure , velocity , concentration , porosity , the numerical solution of pressure , concentration , porosity is defined at the center of the grid unit; the numerical solution of velocity is defined at the midpoint of the four edges of the grid unit, wherein the component of velocity about the direction, namely the velocity in the direction, is defined at the midpoint of the tangential edge of the grid unit, and the component of velocity about the direction is defined at the midpoint of the radial edge of the grid unit.

[0038] According to the preferred embodiment of the present application, the acidification model constructed in step 1 is numerically discretized to form a large linear equation group; finally, the solution is obtained by combining the set parameters, and the acid etching wormhole numerical simulation is realized; including:

[0039] In order to simplify the writing of the subsequent discrete format, first, the variable simplification rule is defined, the discrete function with independent variable and function value at the appropriate discrete node is represented as , represents the moment, and the function is at the grid point where are the grid points are the coordinates in directions and directions, take , take ; note that the superscript indicates the time instant and the subscript indicates the spatial location ; the superscript can be omitted if no ambiguity arises ; if is a vector function, then we also need to use the superscript to distinguish which component in which direction is considered; the notation for the other variables is similar

[0040] The symbolic definition of the difference quotient in place of the derivative is given by , , , , for a function about , , the derivative at a grid point, it is defined as follows:

[0041]

[0042] where , denotes the approximation of the derivative at the midpoint of the edge of the grid cell by the difference quotient of the function values at the center of the grid cell; , denotes the approximation of the derivative at the center of the grid cell by the difference quotient of the function values at the midpoint of the edge of the grid cell;

[0043] The definitions of the interpolation operator and the square root mean operator, including:

[0044] For a point , assume ; note that the direction is a periodic boundary, the bilinear interpolation operator is defined as follows:

[0045] When , there are two extrapolation formulas as follows:

[0046]

[0047] where,​ ;

[0048] For a discrete function , define a piecewise constant function on , satisfying: ;

[0049] For a vector function , whose components are pairs of discrete functions , define the interpolation operator as:

[0050] ;

[0051] where denote the components of the interpolation operator in the direction and the direction, respectively, denotes the function value of at point , and denotes the function value of at point ;

[0052] Let be the norm function of vector , and have ;

[0053] Define the square root mean operator and as follows:

[0054]

[0055] Specifically, it is expressed as:

[0056]

[0057] Use to represent the corresponding numerical solution of , that is, P is the numerical solution of pressure p, W is the numerical solution of fluid velocity u, is the numerical solution of acid concentration , Q is the numerical solution of auxiliary variable q, and T is the numerical solution of rock porosity ; where the superscript indicates the time layer; the velocity and auxiliary variable are vectors, so the superscript is used to distinguish the direction.

[0058] According to the present application, the acidification model constructed in step 1 is preferably numerically discretized by using backward Euler in time and staggered grid finite difference method in space, to obtain the staggered grid finite difference format in polar coordinates, i.e. equations (10)-(16), including:

[0059] A, known and , represent the value of the concentration , the velocity component about the direction , and the porosity at the grid point at the time , obtained according to the equation of the change of porosity with time (4), as follows:

[0060] ; (10)

[0061] wherein, , ; represents the value of the porosity of the rock at the grid point at the initial time;

[0062] B, known , wherein represents the value of the velocity component about the direction at the grid point at the time , represents the value of the velocity component about the direction at the grid point at the time , represents the value of the pressure of the fluid at the grid point at the time ; According to the momentum conservation equation (1) and the mass conservation equation (2) describing the fluid flow, we obtain , as follows:

[0063]

[0064] (11)

[0065] (12)

[0066] ; (13) ​​​

[0067] where, represents an interpolation operator, represents the permeability of the rock, represents the Forchheimer coefficient;

[0068] C, known , according to the acid concentration reaction transport equation formula (9) and auxiliary variable definition, get , as follows:

[0069] (14)

[0070] ; (15)

[0071] ; (16)

[0072] where, , , are numerical approximations of respectively;

[0073] D, initial boundary conditions:

[0074] (17)

[0075] In combination with the above numerical format, using the initial boundary value conditions given in step D, namely formula (17), from start time loop, solving the explicit equation, namely formula (10) and large linear equations, namely formula (11)-(13) and (14)-(16), repeat the three steps of A, B, C, until that is, the numerical solution of the four unknown variables, namely the porosity , pressure , velocity , concentration at the final time, through the staggered grid finite difference method to solve the acidification model, get the numerical solution, can simulate the distribution of rock porosity at any discrete time point, and then draw the evolution of rock porosity with time.

[0076] A computer device comprising a memory and a processor, the memory stores a computer program, the processor executes the computer program to realize the steps of the acidification model based on the DBF framework to realize the numerical simulation method of radial flow of porous medium fluid.

[0077] A computer readable storage medium, having stored thereon a computer program, the computer program being executed by a processor to implement steps of a numerical simulation method for radial flow of fluid in porous media based on an acidification model of a DBF framework.

[0078] The second aspect of the present application provides a system for numerical simulation of radial flow of fluid in porous media based on an acidification model of a DBF framework, comprising:

[0079] An acidification model establishing module configured to construct an acidification model for describing fluid flow, solute reaction transport and rock property change;

[0080] A grid division module configured to divide a simulation time and a ring-shaped solving region into grids, and then define different variables of the acidification model at different positions of the grid cells;

[0081] A numerical solving module configured to numerically discretize the acidification model constructed by the acidification model establishing module to form a large linear equation system, and finally solve the equation system in combination with set parameters to realize numerical simulation of acid etching wormholes.

[0082] The present application has the following beneficial effects:

[0083] 1. Most of the existing acidification models are mathematical models based on the Darcy framework in the Cartesian coordinate system. However, considering that the porosity of the rock changes over time during the acidification process of the carbonate rock matrix and the flow direction of the acid liquid flowing into the carbonate rock matrix through the wellbore is radial, the present application establishes an acidification model based on the Darcy-Brinkman-Forchheimer framework in the polar coordinate system, which can more accurately depict the flow phenomenon of non-Darcy seepage in porous media during the acidification process of the carbonate rock matrix, and thus realize more accurate numerical simulation of acid etching wormholes, thereby playing an important guiding role in oil and gas exploitation.

[0084] 2. The present application proposes to use the staggered grid finite difference method to discretely solve the acidification model in the polar coordinate system, and one of its major features is that different variables can be stored at different positions of the grid. This method not only can well maintain the original physical properties of the model, such as mass conservation and momentum conservation, but also can maintain the high precision of the four unknown quantities of porosity, pressure, velocity and concentration under non-uniform grids, and the numerical solving process is simple and efficient. BRIEF DESCRIPTION OF DRAWINGS

[0085] Figure 1 A framework diagram of the numerical simulation method for radial flow of fluid in porous media based on the acidification model of the DBF framework of the present application;

[0086] Figure 2 A framework diagram of the numerical simulation method for radial flow of fluid in porous media based on the acidification model of the DBF framework of the present application; The non-uniform grid at the moment;

[0087] Figure 3 For the present application, when the initial porosity and permeability of the matrix obeys a uniform distribution and the acid is continuously injected into the matrix through a circular wellbore, the distribution of the rock porosity at the moment is as follows: The distribution diagram at the moment;

[0088] Figure 4 For the present application, when the initial porosity and permeability of the matrix is only large at individual points and the acid is continuously injected into the matrix through a circular wellbore, the distribution of the rock porosity at the moment is as follows: The distribution diagram at the moment;

[0089] Figure 5 For the present application, when the initial porosity and permeability of the matrix obeys a uniform distribution and the acid is continuously injected into the matrix through a circular wellbore, the distribution of the rock porosity at the moment is as follows: The distribution diagram at the moment;

[0090] Figure 6 For the present application, when the initial porosity and permeability of the matrix is only large at individual points and the acid is continuously injected into the matrix through a circular wellbore, the distribution of the rock porosity at the moment is as follows: The distribution diagram at the moment. DETAILED DESCRIPTION

[0091] The present application is further described below by way of examples and in conjunction with the accompanying drawings, but is not limited thereto.

[0092] Example 1

[0093] A numerical simulation method for radial fluid flow in porous media based on an acidification model of a DBF framework, as shown in the figure, comprises: Figure 1 Step 1: constructing an acidification model for describing fluid flow, solute reaction transport and rock property change;

[0094] Step 2: meshing the simulation time and the annular solving region; and then defining different variables of the acidification model at different positions of the grid cells;

[0095] Step 3: numerically discretizing the acidification model constructed in Step 1 to form a large linear equation system; and finally solving in combination with the set parameters to realize numerical simulation of acid erosion wormholes.

[0096] Example 2

[0097] The numerical simulation method for radial fluid flow in porous media based on the acidification model of the DBF framework according to Example 1 is different in that:

[0098]

[0099] ​​An acidification model for describing fluid flow, solute reaction transport and rock property change is established; including:

[0100] A set of acidification models in two-dimensional polar coordinates based on the DBF framework is established as follows:

[0101] ; (1)

[0102] ; (2)

[0103]

[0104] ; (3)

[0105] ; (4)

[0106] wherein, formula (1) and (2) represent the fluid flow process based on the DBF framework, i.e. the acidification model for describing fluid flow; formula (1) is a vector equation, representing the momentum conservation equation, the right end of formula (1): represents the Darcy term, used to describe the Darcy percolation phenomenon of porous media; the second term on the left end of formula (1): represents the Brinkman term, used to describe the transition flow between boundaries; the last term on the left end of formula (1): represents the Forchheimer term, also known as the inertial term, used to describe the significant inertial effect of fluid when the flow rate is high; formula (2) is the mass conservation equation; formula (3) is the acid concentration reaction transport equation, i.e. the acidification model for describing solute reaction transport; formula (4) reflects the change of rock porosity with time evolution, i.e. the acidification model for describing rock property change;

[0107] wherein, , is the flow radius, with the unit of ; is a circular ring area; is time, is the final time, with the unit of ; is the fluid velocity vector, are the components of the fluid velocity in the polar radius direction and the polar angle direction, with the unit of ; , represent the boundary in the polar radius direction; is the fluid pressure, with the unit of ; is the concentration of the acid, in units of ; is the porosity of the rock, dimensionless; velocity , pressure , concentration , porosity All four of these variables are unknown; is the density of the fluid, in units of ; is the viscosity of the fluid, in units of ; is a pseudo-compressibility, in units of , causes a small change in the density of the fluid during the dissolution of the solute, is a small positive number to ensure that the coefficient matrix is invertible; is the local mass transport coefficient, in units of ; is the injection concentration, in units of ; is the dissolution ability of the acid, in units of ; is the density of the rock, in units of ; is the Forchheimer coefficient, dimensionless; is the permeability of the rock, in units of ; where is the injection rate, is the production rate, in units of ; positive definite matrix is the diffusion coefficient of the acid in the porous medium, in units of , and are components of the velocity in the and directions; for simplicity, assume that is a diagonal matrix, where is the molecular diffusion rate, is the identity matrix, denotes a diagonal matrix with on the diagonal; are given functions, while is a function of the porosity of the rock;

[0108] is the concentration of the acid at the fluid (acid) and solid (rock) interface, in units of , according to first-order kinetic reaction, the concentration of the acid at the fluid-solid interface Concentration of acid in the fluid There is a relationship as follows:

[0109] (5)

[0110] wherein, is a surface reaction rate constant, and the unit is ;

[0111] The pore-scale model depicting the change of rock properties is as follows:

[0112] (6)

[0113] (7)

[0114] wherein, formula (6) is established by Carman-Kozeny, and is used to depict the relationship between rock porosity and permeability, denotes permeability, and are initial porosity and initial permeability of the rock respectively; the porosity and the permeability are calculated to obtain wherein, is the interface area (specific surface area) of a unit volume of medium for reaction, and the unit is , is the initial interface area;

[0115] The boundary and initial conditions are as follows:

[0116] (8)

[0117] wherein, the two conditions in the first row: are boundary conditions, indicating that the velocity and the concentration gradient are both 0 at the inlet boundary and the outlet boundary; the following four are all initial conditions, , , , respectively denote the initial distribution functions of the velocity, the pressure, the concentration and the porosity in the region ;

[0118] The acidification model formula (1)-(4) based on the DBF framework in the two-dimensional polar coordinates and the pore-scale model formula (6)-(7) and the initial boundary condition formula (8) are combined together to construct the acidification model of the application; the flow properties of the non-Darcy seepage of the porous medium in the acidification process of the carbonate rock matrix are comprehensively considered, and compared with the prior art, the acidification model established in the application can more accurately depict the acid erosion wormhole process;

[0119] For simplifying the equations and the writing of the following discrete schemes, define an auxiliary variable where, and combine equations (2), (5)-(7), the concentration reaction transport equation (3) of acid is simplified and transformed as follows:

[0120] (9).

[0121] Divide the simulation time and the annular solving region into grids; define different variables of the acidification model at different positions of the grid cells; including:

[0122] Divide the total simulation time (the final time) into parts, then the time step and the time ; in order to better simulate the actual situation, for a given simulation region , i.e. the annular solving region, perform non-uniform grid division, divide the length of the region in the direction into parts, divide the length of the region in the direction into parts, form grid cells; the grid point coordinates are: and ; the midpoint , of the grid cell edges in the direction and the direction and the division step , , are respectively:

[0123]

[0124] The staggered grid is characterized by storing different variables at different positions of the grid cells; for the nonlinear strong coupling acidification model of pressure , velocity , concentration , and porosity established in step 1, i.e. equations (1)-(4), define the numerical solution of pressure , concentration , and porosity at the center of the grid cell; define the numerical solution of velocity at the midpoint of the four edges of the grid cell, wherein the component of velocity in the direction is The velocity in the direction is defined at the midpoint of the tangential edge of the grid cell, and the velocity component in the direction is defined at the midpoint of the radial edge of the grid cell. The velocity in the direction is defined at the midpoint of the tangential edge of the grid cell, and the velocity component in the direction is defined at the midpoint of the radial edge of the grid cell.

[0125] The acidification model constructed in step 1 is numerically discretized to form a large linear equation system, and finally combined with the set parameters to solve, realizing the numerical simulation of acid etching wormhole; including:

[0126] To simplify the writing of the subsequent discrete format, first define the variable simplification rules, including:

[0127] The independent variable is and the discrete function with function value at the appropriate discrete node is represented as , represents at time , the value of the function at the grid point ; wherein are the coordinates of the grid point in the direction and the direction, , , ; let , wherein the superscript indicates the time , and the subscript indicates the spatial position , and the superscript can usually be omitted without ambiguity ; if is a vector function, then on this basis, the superscript is also needed to distinguish the component in which direction; the simplification of the remaining variables is similar;

[0128] The symbolic definition of replacing the derivative with the difference quotient is given, denoted as , , , , is the derivative of the function with respect to , , at a certain grid point, then its definition is as follows:

[0129]

[0130] wherein , indicates that the derivative at the midpoint of the grid cell edge is approximated by the difference quotient of the function value at the center of the grid cell; , denotes the difference quotient of the function values at the midpoints of the edges of the grid cell to approximate the derivative at the center of the grid cell;

[0131] The definitions of the interpolation operator and the square root mean operator include:

[0132] For a point , assume , it is worth noting that The direction is the periodic boundary, , the value of is defined as follows bilinear interpolation operator :

[0133] When , there are two point extrapolation formulas as follows:

[0134]

[0135] Where, ;

[0136] For a discrete function , define a piecewise constant function on , which satisfies: ;

[0137] Then, for a vector function , whose components are a pair of discrete functions , define the interpolation operator as follows:

[0138] ;

[0139] Where, denotes the components of the interpolation operator in the direction and the direction, respectively, denotes the function value of at the point , and denotes the function value of at the point ;

[0140] Let be the norm function of the vector , and ;

[0141] Define the square root mean operator and as follows:

[0142]

[0143] is expressed as:

[0144]

[0145] is expressed as: The corresponding numerical solution, i.e. P is the numerical solution of pressure p, W is the numerical solution of fluid velocity u, is the numerical solution of the concentration of acid Q is the numerical solution of auxiliary variable q, T is the numerical solution of rock porosity ; wherein the superscript represents the time layer; the velocity and auxiliary variable are vectors, so the superscript is used to distinguish the direction;

[0146] The acidification model constructed in step 1 is numerically discretized by using backward Euler in time and staggered grid finite difference method in space to obtain the staggered grid finite difference format in polar coordinates, i.e. formula (10)-(16), including:

[0147] A, known and , represent the value of the concentration , porosity at the grid point , according to the equation of porosity changing with time formula (4), to obtain , as follows:

[0148] ; (10)

[0149] wherein, , ; represents the initial time, the value of rock porosity at the grid point ;

[0150] B, known , wherein represents the value of the component of velocity about direction at the grid point , represents the value of the component of velocity about direction at the grid point , represents the value of the component of velocity about direction at the grid point , represents the value of the component of velocity about ​​​moment , the pressure of the fluid at the grid point ;

[0151] According to the momentum conservation equation formula (1) and the mass conservation equation formula (2) describing the fluid flow, we get , as follows:

[0152] (11)

[0153] (12)

[0154] ; (13)

[0155] wherein, represents an interpolation operator, represents the permeability of the rock, represents the Forchheimer coefficient;

[0156] C, known , according to the acid concentration reaction transmission equation formula (9) and the definition of auxiliary variable , we get , as follows:

[0157] (14)

[0158] ; (15)

[0159] ; (16)

[0160] wherein, , , are numerical approximations to , respectively;

[0161] D, initial boundary conditions:

[0162] (17)

[0163] Combined with the above numerical format, using the initial boundary value conditions given in step D, i.e. formula (17), starting from time loop, solving the explicit equation, i.e. formula (10), and the large linear equation system, i.e. formula (11)-(13) and (14)-(16), repeating steps A, B, C, until , i.e. the final time four unknown variables, i.e. porosity , pressure , velocity , concentration The numerical solution of the acidification model is obtained by the staggered grid finite difference method, and the distribution of the rock porosity at any discrete time point can be simulated, and then the evolution of the rock porosity with time can be drawn.

[0164] The numerical simulation of the acidification and acid-etching wormhole in the carbonate rock matrix is realized through steps 1, 2 and 3; based on the flow characteristics of the non-Darcy seepage of the porous medium, a two-dimensional polar coordinate system Darcy-Brinkman-Forchheimer nonlinear strong coupling acidification model is established, and the staggered grid finite difference method is used for numerical discretization and efficient solution, so that the numerical simulation of the acid-etching wormhole process is more practical and more accurate.

[0165] In order to verify the accuracy of the acidification model established in the application, the feasibility of the numerical solution algorithm and the numerical simulation effect, four numerical experiments are listed, and the computer software is used to write programs for numerical calculation and simulation. First, Table 1 shows some physical parameters required for numerical simulation and their values. Then in the four numerical examples, the initial porosity, permeability distribution and different acid injection range with different characteristics are given, and the numerical simulation effect verifies the effectiveness of the technical scheme of the application.

[0166] Table 1

[0167]

[0168] In the numerical simulation of the acid-etching wormhole reservoir, the simulation time interval is , and a total of about 46 days are simulated, the time step is , the simulation area radius is , the wellbore radius is , that is , and a non-uniform grid of units is used, as shown in Figure 2 , wherein each small quadrilateral represents a grid unit, and the vertices of all grid units are called grid points. The initial pressure is set to , the initial concentration is 0. The initial permeability is uniformly distributed in the range , the initial porosity is uniformly distributed in the range , and the acid injection and production rates are as follows:

[0169]

[0170] When the initial porosity and initial permeability of the carbonate rock matrix obey a uniform distribution, the acid is continuously injected into the carbonate rock matrix through a circular wellbore, the acid reacts with the rock to change the structure of the rock, and the porosity and permeability of the rock are increased as the matrix acidification process advances. Since the high-porosity region has less resistance to fluid, the acid preferentially flows into this region, forming randomly distributed high-conductivity channels (branching out of the dominant wormhole), and the acid-etched wormhole numerical simulation results are as shown in Figure 3 .

[0171] Example 3

[0172] In this example, the simulation time interval is , and a total of about 24 days are simulated, with a time step of . The simulation area radius is , the wellbore radius is , i.e. , and a non-uniform grid of elements is used, as shown in FIG. 2. The initial pressure is set to , and the initial concentration is 0. The initial porosity and permeability and the injection and production rates are as follows:

[0173]

[0174]

[0175] When the initial porosity and permeability of the carbonate rock matrix are large at only a few points, the acid is continuously injected into the carbonate rock matrix through a circular wellbore, and as the acidification time evolves, the acid extends radially outward faster at the points where the initial porosity and permeability are large, thereby forming a narrow high-conductivity channel. This trend is consistent with the expected results, and the acid-etched wormhole numerical simulation results are as shown in Figure 4 .

[0176] Example 4

[0177] In this example, the simulation time interval is , and a total of about 58 days are simulated, with a time step of . The simulation area radius is , the wellbore radius is , i.e. , and a non-uniform grid of elements is used, as shown in FIG. 2. The initial pressure is set to , and the initial concentration is 0. The initial permeability obeys a uniform distribution in the range , and the initial porosity obeys a uniform distribution in the range uniform distribution, and the injection and production rates of the acid are as follows:

[0178]

[0179] When the initial porosity and permeability of the carbonate rock matrix obey a uniform distribution, if only through a circular wellbore, the acid is continuously injected into the carbonate rock matrix, with the evolution of the acidizing time, the acid continuously diffuses and preferentially flows into the high porosity area, and due to the fact that the acid extends radially outward faster in the high porosity area, a randomly distributed high conductivity channel (a branched out dominant wormhole) is gradually formed, and the acid etching wormhole numerical simulation result is as shown in Figure 5

[0180] Example 5

[0181] In this example, the simulation time interval is , and a total of about 92 days is simulated, with a time step of . The simulation area radius is , the wellbore radius is , i.e. , and a non-uniform grid of units is used, as shown in Fig. 2. The initial pressure is set to , and the initial concentration is 0. The initial porosity, permeability, and injection and production rates are as follows:

[0182]

[0183]

[0184] When the initial porosity and permeability of the carbonate rock matrix are large in only two continuous areas along the radial direction, similar to the existence of two small cracks in the carbonate rock matrix, if only through a small range of wellbores near the two areas, acid is continuously injected into the matrix with the above characteristics, due to the fact that the porosity and permeability at the cracks are large, equivalent to a pipeline, a certain conductivity effect is achieved, then the acid preferentially flows to this area, with the evolution of the acidizing time, a large amount of acid flows into the two high porosity channels and forms a branched out branch, and the final acid etching wormhole numerical simulation result is as shown in Figure 6

[0185] The numerical simulation results of the above examples effectively demonstrate the accuracy and feasibility of the acidizing model and numerical solving method of the present application. Therefore, the technical scheme of the present application not only can more accurately numerically simulate the acid etching wormhole process of the carbonate rock matrix acidizing, but also can play an important guiding role in oil and gas exploitation.

[0186] Example 6

[0187] ​​A computer device comprises a memory and a processor, the memory stores a computer program, and the processor implements the steps of the numerical simulation method for radial flow of fluid in porous media based on the DBF framework acidification model according to any one of embodiments 1-5 when executing the computer program.

[0188] Embodiment 7

[0189] A computer readable storage medium, having stored thereon a computer program, the computer program, when executed by a processor, implements the steps of the numerical simulation method for radial flow of fluid in porous media based on the DBF framework acidification model according to any one of embodiments 1-5.

[0190] Embodiment 8

[0191] A system for numerical simulation of radial flow of fluid in porous media based on a DBF framework acidification model, comprising:

[0192] An acidification model establishing module configured to construct an acidification model for describing fluid flow, solute reaction transport and rock property change;

[0193] A grid division module configured to divide a simulation time and a ring-type solving region into grids, and then define different variables of the acidification model at different positions of the grid cells;

[0194] A numerical solving module configured to numerically discretize the acidification model constructed by the acidification model establishing module to form a large linear equation system, and finally solve the equation system combined with set parameters to realize numerical simulation of acid etching wormholes.

Claims

1. A method for numerical simulation of radial flow in porous media based on DBF framework for acidizing model implementation, characterized in that, Comprising: Step 1: constructing an acidification model for describing fluid flow, solute reaction transport and rock property variation; Comprising: A set of two-dimensional polar coordinate-based acidification models under the DBF framework are established as follows: wherein, formula (1) and (2) represent the fluid flow process based on the DBF framework, i.e., the acidification model for describing fluid flow; formula (1) is a vector equation, representing the momentum conservation equation, and the right end of formula (1): represents the Darcy term, used to depict the Darcy seepage phenomenon of the porous medium; the second term on the left end of formula (1): represents the Brinkman term, used to depict the transition flow between boundaries; the last term on the left end of formula (1): represents the Forchheimer term, used to describe the significant inertial effect of the fluid when the flow rate is high; formula (2) is the mass conservation equation; formula (3) is the acid concentration reaction transport equation, i.e., the acidification model for describing solute reaction transport; formula (4) reflects the change of rock porosity with time evolution, i.e., the acidification model for describing the change of rock properties. where (r, θ) ∈ Ω = {(r, θ) | r in ≤ r ≤ r out , 0 ≤ θ ≤ 2π}, t ∈ J = [0, T], r is the flow radius with the unit of m; θ is the flow angle; Ω ∈ R 2 is a circular ring region; t is the time, T is the final time with the unit of s; u is the fluid velocity vector, u r , u θ are the components of the fluid velocity u in the radial direction and the angular direction, respectively, with the unit of m / s; r in , r out represent the boundary in the radial direction; p is the fluid pressure with the unit of Pa; c f is the concentration of acid with the unit of mol / m 3 ; φ is the porosity of rock, dimensionless; the four variables of velocity u, pressure p, concentration c f , and porosity φ are all unknown; ρ is the fluid density with the unit of kg / m 3 ; μ is the fluid viscosity with the unit of Pa·s; ε is a pseudo-compression coefficient with the unit of 1 / Pa, ε represents a positive number; k c is the local mass transport coefficient with the unit of m / s; c I is the injection concentration with the unit of mol / m 3 ; α is the dissolution ability of acid with the unit of kg / mol; ρ s is the density of rock with the unit of kg / m 3 ; is the Forchheimer coefficient, dimensionless; K(φ) is the permeability of rock with the unit of m 2 ; f = f I + f P , where f I is the injection rate and f P is the production rate, with the unit of m / s; the positive definite matrix D is the diffusion coefficient of acid in porous media with the unit of m 2 / s, D r and D θ are the components of D in the r and θ directions; it is assumed that D = d mol I = diag(D ll ) is a diagonal matrix, where d mol is the molecular diffusion rate, I is the unit matrix, and diag(D ll ), (l = 1, 2) represents the diagonal matrix with D ll as the diagonal elements; c s is the concentration of acid at the fluid-solid interface, in mol / m 3 , according to the first-order kinetic reaction, the concentration of acid at the fluid-solid interface c s and the concentration of acid in the fluid c f have the following relationship: where k s is the surface reaction rate constant with units of m / s; The pore-scale model for depicting rock property variation is as follows: wherein formula (6) is used to depict the relationship between rock porosity and permeability, K represents permeability, and φ0and K0are initial porosity and initial permeability of the rock, respectively; porosity and permeability are calculated to obtain a v wherein a v is the interface area per unit volume of the medium for reaction, with a unit of 1 / m, and a0is the initial interface area; The boundary and initial conditions are as follows: where the two conditions in the first row: u = 0, are boundary conditions, indicating that the gradients of velocity and concentration are both 0 at the inlet and outlet boundaries; the latter four are initial conditions, u0(r, θ), p0(r, θ), c f0 (r, θ), φ0(r, θ) represent the initial distribution functions of velocity, pressure, concentration, and porosity in the region Ω, respectively. The two-dimensional polar coordinate-based acidification model formula (1)-(4) under the DBF framework is combined with the pore-scale model formula (6)-(7) and the initial boundary condition formula (8) to form an acidification model; Define a helper variable where, In combination with equations (2), (5)-(7), the concentration reaction transport equation (3) for acid simplifies to the following: Step 2: grid partitioning for the simulation time and the annular solving region; then defining different variables of the acidification model at different positions of the grid cells; Step 3: numerical discretization of the acidification model constructed in Step 1 to form a large linear equation system; finally, combining with the set parameters to solve, realizing the numerical simulation of acid-etched wormholes.

2. The method of claim 1, wherein the acidizing model is implemented based on a DBF framework for numerical simulation of radial fluid flow in porous media. The simulation time J = [0, T] and the annular solution region Ω = {(r, θ) | r in ≤ r ≤ r out , 0 ≤ θ ≤ 2π} are meshed. Then defining different variables of the acidification model at different positions of the grid cells; comprising: Divide the total simulation time T into N equal parts, then the time step is... And at the nth time t n =nΔt, n≤N; For the given simulation region Ω, i.e., the annular solution region, perform non-uniform meshing, dividing the region length in the r direction into N... r The region length in the θ direction is divided into N parts. θ Parts, forming N r ×N θ There are 1 grid cell, and the grid points are: And there are Let r be the midpoint of the grid cell edge in the r and θ directions. i θ j and the step size of the subdivision They are respectively: The nonlinear strongly coupled acidification model, i.e., equations (1)-(4), established for step 1 regarding pressure p, velocity u, concentration c f , and porosity φ, defines the numerical solution of pressure p, concentration c f , and porosity φ at the center of the grid cell; and defines the numerical solution of velocity u at the midpoints of the four edges of the grid cell, wherein the component of velocity with respect to r direction, i.e., velocity in r direction, is defined at the midpoints of the tangential edges of the grid cell, and the component of velocity with respect to θ direction is defined at the midpoints of the radial edges of the grid cell.

3. The method of claim 2, wherein the acidizing model is implemented based on a DBF framework for numerical simulation of radial fluid flow in porous media. Numerical discretization of the acidification model constructed in Step 1 to form a large linear equation system; finally, combining with the set parameters to solve, realizing the numerical simulation of acid-etched wormholes; comprising: A variable definition simplification rule represents a discrete function with arguments r, θ, t and function values at appropriate discrete nodes as g(r, θ, t), g(r l ,θ m ,t n ) represents the value of the function g at grid point (r l ,θ m ) at time t n ; where r l , θ m are the coordinates of the grid point (r l ,θ m ) in the r direction and the θ direction respectively, and l takes i, m takes j, (i, j are integers); note where the superscript represents the n-th time t n , and the subscript represents the spatial position (r l ,θ m ); The symbolic definition of the derivative is given by the difference quotient The derivative of a function g with respect to r, θ, t at a grid point is defined as follows: where d r , d θ denotes the approximation of the derivative at the midpoint of the edge of the grid cell using the difference quotient of the function values at the center of the grid cell; denotes the approximation of the derivative at the center of the grid cell using the difference quotient of the function values at the midpoint of the edge of the grid cell; The definitions of the interpolation operator and the square root mean operator, including: For a point (r, θ), assume r ∈ [r i ,r i+1 ],θ∈[θ j ,θ j+1 ], i = 1, 2, ..., N r -1,j=1,2,…,N θ The θ direction is the periodic boundary. Use p i,j The value is defined as follows for the bilinear interpolation operator Π h p: When Or Then, there are two extrapolation formulas as follows: where p 1,j = p(r1, θ j ), p 2,j = p(r2, θ j ), For the discrete function q i,j = q(r i , θ j ), define a piecewise constant function on Ω satisfying: Π h q(x, y) = q i,j (x, y) for (x, y) e Ω i,j ; For a vector function V = (V r ,V θ ), whose components V r , V θ are a pair of discrete functions and the interpolation operator I1 is defined as follows: wherein, wherein denote the components of the interpolation operator I1in the r and θ direction, respectively, denote the function value of V r at the point denote the function value of V denote the function value of V θ at the point denote the function value of V Let N(U, V) be the norm function of the vector (U, V) and have Definition of the square root mean operator S r and S θ as shown in the following equation: Specifically represented as: use to indicate The corresponding numerical solutions are: P is the numerical solution for pressure p, W is the numerical solution for fluid velocity u, and Ψ is the concentration c of the acid. f Numerical solution, Auxiliary variable The numerical solution is T, which is the numerical solution of rock porosity φ; where the superscript n represents the time layer; the velocity and auxiliary variables are vectors, so the superscripts r and θ are used to distinguish the directions.

4. The method of claim 3, wherein the acidizing model is implemented based on a DBF framework for numerical simulation of radial fluid flow in porous media. Numerical discretization of the acidification model constructed in Step 1 using the backward Euler in time and the staggered grid finite difference method in space to obtain the staggered grid finite difference format in polar coordinates, comprising: A, known and denotes the (n-1)th time t n-1 , the value of the concentration Ψ, the porosity T at the grid point (r i , θ j ) is obtained from the equation (4) for the time variation of the porosity as follows: wherein φ 0,i,j denotes the value of the rock porosity φ at the grid point (r i ,θ j ) at the initial time instant B. Known in This represents the (n-1)th time t. n-1 The velocity component W in the r direction r At grid points The value at that location, This represents the (n-1)th time t. n-1 The component of velocity W in the θ direction θ At grid points The value at that location, This represents the (n-1)th time t. n-1 The fluid pressure P at grid point (r) i ,θ j The value at () According to the momentum conservation equation formula (1) and the mass conservation equation formula (2) describing fluid flow, the following is obtained As follows: wherein Π h represents an interpolation operator, K represents the permeability of the rock, F(Π h T n ) represents the Forchheimer coefficient; C, known According to the reaction transport equation formula (9) and the definition of the auxiliary variable , we get As follows: wherein are numerical approximations of respectively. D. Initial boundary conditions: Combining the above numerical format, using the initial boundary condition given in Step D, i.e. formula (17), starting the time loop from n=1, n=1, 2, …, N, sequentially solving the explicit equation, i.e. formula (10), and the large linear equation system, i.e. formula (11)-(13) and (14)-(16), repeating the three steps of A, B, and C until n=N, i.e. obtaining the numerical solutions of the four unknown variables, i.e. porosity T, pressure P, velocity W, and concentration Ψ at the final time.

5. A computer device comprising a memory and a processor, the memory storing a computer program, characterized in that, The processor executes the computer program to realize the steps of the DBF framework-based acidification model of any one of claims 1-4 to realize the steps of the numerical simulation method of radial fluid flow in porous media.

6. A computer-readable storage medium having stored thereon a computer program, characterized in that, The computer program is executed by the processor to realize the steps of the DBF framework-based acidification model of any one of claims 1-4 to realize the steps of the numerical simulation method of radial fluid flow in porous media.

7. A system for numerical simulation of radial flow in porous media based on DBF framework acidizing model implementation, characterized in that, Comprising: The acidification model establishment module is configured to construct an acidification model for describing fluid flow, solute reaction transport and rock property variation; Comprising: A set of two-dimensional polar coordinate-based acidification models under the DBF framework are established as follows: wherein, formula (1) and (2) represent the fluid flow process based on the DBF framework, i.e., the acidification model for describing fluid flow; formula (1) is a vector equation, representing the momentum conservation equation, and the right end of formula (1): represents the Darcy term, used to depict the Darcy seepage phenomenon of the porous medium; the second term on the left end of formula (1): represents the Brinkman term, used to depict the transition flow between boundaries; the last term on the left end of formula (1): represents the Forchheimer term, used to describe the significant inertial effect of the fluid when the flow rate is high; formula (2) is the mass conservation equation; formula (3) is the acid concentration reaction transport equation, i.e., the acidification model for describing solute reaction transport; formula (4) reflects the change of rock porosity with time evolution, i.e., the acidification model for describing the change of rock properties. where (r, θ) ∈ Ω = {(r, θ) | r in ≤ r ≤ r out , 0 ≤ θ ≤ 2π}, t ∈ J = [0, T], r is the flow radius with the unit of m; θ is the flow angle; Ω ∈ R 2 is an annular region; t is time, T is the final time with the unit of s; u is the fluid velocity vector, u r , u θ are the components of the fluid velocity u in the radial direction r and the angular direction θ, respectively, with the unit of m / s; r in , r out represent the boundaries in the radial direction r; p is the fluid pressure with the unit of Pa; c f is the concentration of acid with the unit of mol / m 3 ; φ is the porosity of rock, dimensionless; the four variables of velocity u, pressure p, concentration c f , and porosity φ are all unknown; ρ is the fluid density with the unit of kg / m 3 ; μ is the fluid viscosity with the unit of Pa·s; ε is a pseudo-compression coefficient with the unit of 1 / Pa, ε represents a positive number, k c is the local mass transport coefficient with the unit of m / s; c I is the injection concentration with the unit of mol / m 3 ; α is the dissolution ability of acid with the unit of kg / mol; ρ s is the density of rock with the unit of kg / m 3 ; is the Forchheimer coefficient, dimensionless; K(φ) is the permeability of rock with the unit of m 2 ; f = f I + f P , where f I is the injection rate, f P is the production rate, with the unit of m / s; the positive definite matrix D is the diffusion coefficient of acid in porous media with the unit of m 2 / s, D r and D θ are the components of D in the r and θ directions; it is assumed that D = d mol I = diag(D ll ) is a diagonal matrix, where d mol is the molecular diffusion rate, I is the unit matrix, diag(D ll ), (l = 1, 2) represents the diagonal matrix with D ll as the diagonal elements; c s is the concentration of acid at the fluid-solid interface, in mol / m 3 , according to the first order kinetic reaction, the concentration of acid at the fluid-solid interface c s is related to the concentration of acid in the fluid c f by the following relationship: where k s is the surface reaction rate constant with units of m / s; The pore-scale model for depicting rock property variation is as follows: wherein formula (6) is used to depict the relationship between rock porosity and permeability, K represents permeability, and φ0and K0are initial porosity and initial permeability of the rock, respectively; porosity and permeability are calculated to obtain a v wherein a v is the interface area per unit volume of the medium for reaction, with a unit of 1 / m, and a0is the initial interface area; The boundary and initial conditions are as follows: where the two conditions in the first row: u = 0, are boundary conditions, indicating that the gradients of velocity and concentration are both 0 at the inlet and outlet boundaries; the latter four are initial conditions, u0(r, θ), p0(r, θ), c f0 (r, θ), φ0(r, θ) represent the initial distribution functions of velocity, pressure, concentration, and porosity in the region Ω, respectively. The two-dimensional polar coordinate-based acidification model formula (1)-(4) under the DBF framework is combined with the pore-scale model formula (6)-(7) and the initial boundary condition formula (8) to form an acidification model; Define a helper variable where, In combination with equations (2), (5)-(7), the concentration reaction transport equation (3) for acid simplifies to the following: The grid partitioning module is configured to perform grid partitioning for the simulation time and the annular solving region; then defining different variables of the acidification model at different positions of the grid cells; The numerical solution module is configured to: perform numerical discretization on the acidification model built by the acidification model building module to form a large linear equation set; and finally solve in combination with set parameters to realize acid etching burrow numerical simulation.

Images