It is crucial for vuggy carbonate formation to enhance near wellbore permeability to connect formation, release damage and improve productivity. Carbonate minerals can be dissolved by hydrochloric acid and the conductivity path will be formed between formation and wellbore  . The wormhole is formed when the acid dissolving reacts with the carbonate mineral and its shape plays an important role in acidizing. Many researchers have revealed the relationship between acid type, injection velocity, acid volume and wormhole development rule. The acid/rock reaction mechanism indicates that the reaction is determined by convection, diffusion and mass transfer      . T. A. Kremleva  considered that the controlling acid injection volume and reaction time are the most significant influencing factors. Mohan K. R   (2005) researched matrix acidizing in carbonate formation by combining the simulation and laboratory experiment results. A two-dimensional numerical model was established based on Darcy flow and microscopic flow in pore. And the effects of porosity heterogeneity, injection rate and DamKahler number on wormhole propagation patterns were analyzed. The acid-etched wormhole can be formed only with appropriate acid injection volume, which will connect formation and improve permeability. Vemuri Balakotaiah  established a radial numerical model for wormhole based on fractal geometry. All the models mentioned above are for acidizing simulation in single medium. In fact, in the formation of vuggy carbonate, acid flow in the vug system is described by Navier-Stokes equation, while that in matrix is described by Darcy law. Omer Izgec  deduced pressure field and acidizing models for vuggy carbonate formation according to the above theory. Actual 3-D porosity model was also established by experiment and CT scanning technology. Flow pattern in fracture-vug media can be accurately described by the equivalent continuum model (Li Yajun) combined with discrete medium model     . In this paper, a numerical model for vuggy carbonate formation is deduced based on equivalent permeability, which avoids distractions of pressure field and seepage variation between fracture-vug and matrix. The simulation result provides an important guidance to optimize the operation parameters of matrix acidizing in carbonate formation.
2. Model Establishment and Deduction
2.1. Velocity Field Equation is
Since the acid flow pattern conforms to Darcy laws, the velocity can be described by equivalent permeability  :
Then, the pressure field equation is:
Assuming porosity will not change per time step, then (3) is expressed as a steady-state equation:
2.2. The Convection Diffusion Equation
Based on the microelement analysis (Figure 1), Darcy flow and the microscopic flow in pore, the molarity-balance equation of acid/rock reaction is deduced:
When (5) is divided by, then：
Assuming Δx and Δy approach zero, then (6) is expressed as:
Figure 1. Unit control volume model.
Diffusion coefficient De, which is composed of mass transfer and convection diffusion, can be described by Peclet coefficient Pe  :
Then the convection diffusion equation can be simplified as:
2.3. The Porosity Equation
The porosity variation can be calculated by acid dissolving capacity number α:
By simplifying (10), the porosity equation is expressed as:
(2), (4), (9), (11) constitute the numerical model of matrix acidizing in carbonate formation:
The acid concentration in rock is zero before being acidized. The pressure boundary and concentration boundary are closed when acid is injected into the entry end of rock with the constant volume. And the pressure boundary of the outlet end is a constant and the acid concentration is zero.
Initial conditions are:
Boundary conditions are:
Modan et al.  put forward the following relations:
K―equivalent permeability, can describe heterogeneity and anisotropy;
vx―Darcy velocity (x - z boundary), m/min;
vY―Darcy velocity (y - z boundary), m/min;
av―Pore specific surface, m2/m3;
α―acid dissolving capacity number, mol/mol;
ρs―Rock density, kg/m3;
F―formation resistivity coefficient;
Qinj―acid injection volume, m3/min;
Δz―core thickness, m;
rp―throat pore diameter, m;
Sh∞―approximate Sherwood number, Sh∞ = 2;
Rep―Reynolds number, Rep ≈ 0.2;
Sc―Schmidt number, Sc = v/Dm;
kc―mass transfer coefficient, m/min;
Cf―acid mass concentration, kg/m3;
Cs―acid mass concentration in vuggy surface, kg/m3.
3. Difference Solution
The pressure equation can be simplified from (4) by the centered difference as the nine-diagonal equation:
The convection diffusion equation can be simplified from (9) by the centered difference as the five-diagonal equation:
The step of the central point is calculated by the arithmetic method:
The porosity equation can be simplified from (11) by the forward difference as a linear equation:
1) Calculate equivalent permeability tensors of each grid at initial time;
2) Integrate the equivalent permeability tensors into Equation (16), and solve the five-diagonal matrix to achieve pressure distribution in the core;
3) According to Equation (17) and the pressure distribution, the acid concentration in the core is received.
4) The reservoir porosity is solved by combining the reaction equation Equation (17) and the porosity equation Equation (18).
5) Return to step. 1 to calculate the equivalent permeability tensor of each grid next time, recycling the above steps until the end of the injection.
4. Comparison and Analysis
Omer Izgec et al.  cut a core into 150 ~ 200 slices and scanned the rock slices by CT. A three-dimensional model of porosity distribution and wormhole propagation was established before and after the acidizing.
In this paper, the acid consumption for wormhole breakthrough and differential pressure varying with time is calculated by the proposed model. The following simulation results are based on the parameters in Table 1. The results are similar to the records by experiment. BC21and BC10, as the experimental subjects of Omer Izgec, indicate that this model is accurate and provides a significant guidance to the carbonate formation or vuggy carbonate formation acidizing. The simulation model program image is illustrated as Figure 2.
Table 1. Core Physical Parameters and Test Parameters.
Figure 2. The simulation model program image.
By the CT scanning and computer simulation, the porosity distribution of BC21 and BC10 before and after acidizing is shown in Figure 3(A) and Figure 4(A). Major wormhole and branch wormholes can be observed in the figures. The distribution of macropores before acidizing is marked by red arrows in Figure 3(A-1), while the propagating trajectory of the major wormhole is marked by red arrows in Figure 3(A-2). Since the two trajectories are similar, the flow channel is depended on the rock, minerals and the macropores distribution.
Figure 3. Porosity distribution of core BC21 before and after acidizing.
Figure 4. Porosity distribution of core BC10 before and after acidizing.
Porosity distribution of BC21 and BC10 before and after acidizing by simulation is shown in Figure 3(B) and Figure 4(B). The major wormhole and branch wormholes diagrammed in Figure 3(B-2) and Figure 4(B-2) are the same as the result of Omer Izgec, which indicates that the simulated result shows a good agreement with those of corresponding experiments. The differential pressure drop of the entry end is depended on the wormhole propagation, therefore the differential pressure drop tends to be gentle (Figure 5 and Figure 6) after forming the major wormhole channel.
The differential pressure drops of experimental and simulation result are similar (Figure 5 and Figure 6). The differential pressure drop tends to be gentle when acid volume is 0.028 m3 (BC21) and 0.048 m3 (BC10), which is different from Mohan’s viewpoint   that wormhole breaks through the core when the differential pressure of the entry end is the 0.01 time of the initial differential pressure. It is concluded in this paper that the cross-point of the tangent to the curve is PV (PV is defined as the acid consumption when the major wormhole channel is formed). When the acid volume reaches PV, the major wormhole is formed (Figure 7(A)). When the acid volume is 1.5 times of PV (1.5 times of PV divided by pore volume of core is defined as acid breakthrough volume), the wormhole can break through the core completely (Figure 7(B)), and the shape will not change even continuing injecting acid. And its end will be etched (Figure 7(C)).
Figure 5. Differential pressure drop of core BC21.
Figure 6. Differential pressure drop of core BC10.
Figure 7. Wormhole shape for different injecting volume.
Furthermore, the acid type, concentration and injection velocity influence the wormhole development. More intensive studies will be developed in the follow-up work.
1) In this paper, a novel numerical model for carbonate formation and vuggy carbonate formation acidizing is deduced based on equivalent permeability, which avoids distractions of pressure field and seepage variation between fracture-vuggy and matrix. The model is verified to be correct by the computer simulation.
2) According to the experimental and simulation results, the acid flow channel is depended on the rock, minerals and the macropores distribution. If continuing to inject acid after the wormhole breaking through the core completely, the formed wormhole will not change significantly but the end will be dissolved.
3) However, the main reason for the difference between experimental results and simulations is that, the simulation results calculated through our models are based on 2-D vision. However, the experimental results come from 3-D core flooding experiments.
4) PV is defined as the acid consumption when the major wormhole channel is formed while 1.5 times of PV divided by pore volume of core is defined as the acid breakthrough volume.
5) Vug size and its distribution are the main factors for the formed wormhole and for the wormhole to break through the rock. The simulation result can provide an important guidance to optimize the operation parameters of matrix acidizing in the vuggy carbonate formation.