Simulation of Acid–rock Heterogeneous Flow Reaction Based on the Lattice Boltzmann Method

Matrix acidizing is an essential strategy to maintain or increase the productivity or injectivity of hydrocarbon wells. However, for sandstone reservoirs, the heterogeneous flow reaction mechanism of acid–rock in porous media is very complex because of their complex mineral and chemical compositions. It is often difficult to match real formation conditions by experimental simulation. Also, traditional numerical simulation methods have the disadvantages of complex boundary processing and low computational efficiency. In this study, the lattice Boltzmann method (LBM) was used to establish the heterogeneous flow reaction model of acid–rock from a new perspective, which was solved by MATLAB to obtain the distribution of temperature, concentration of various substances, porosity, and permeability. The simulation results indicate that with increases in injection time and injection speed, the temperature and mass transfer distance of the acid will also increase. Changing the injection time had a more obvious influence on the transfer of temperature and mass than did changing the injection speed. The increasing rates of porosity and permeability in the middle of the flow channel were the highest. The fast-reaction mineral content, hydrofluoric acid injection concentration, and acid injection time had a great influence on the acidizing effect, whereas the slow-reaction mineral content, acid injection temperature, and injection speed had little influence on the acidizing effect. The results suggest that to improve the acidizing effect, priority should be given to improve the HF concentration and acid dose. It will be important for further guiding the optimization of acidizing process design parameters.


Introduction
Acidizing is an important procedure for stable and increased production of oil and gas wells and for facilitating water injection in wells. The purpose of acidizing is to dissolve some of the formation minerals and plugs by injecting suitable acid and additives into a reservoir, thus increasing porosity and permeability, improving reservoir seepage conditions, and increasing oil and gas well productivity [1]. As more low-porosity and low-permeability porous sandstone reservoirs are discovered, more acidification is required to increase production. Porous sandstone reservoirs consist of sand and cements. Because the sand contains feldspar, quartz, mica, and other components, and the cements contain clay, carbonate, metal oxide, and other components, porous media have many mineral components and complex structures. During acidizing, most mineral components react with the acid in liquid-solid heterogeneous chemical reactions. The mechanism of the heterogeneous chemical reactions of these acid-rocks is very complex, and some of the products are precipitates, which cause secondary damage to the formation if not properly treated [2].
Due to the complexity of the heterogeneous flow reaction of acid-rock in porous media, blind chemical stimulation greatly reduces the success rate of sandstone acidizing and seriously hinders the stimulation effect [3]. Therefore, to simulate a flow reaction before an acidizing operation, it is necessary to adopt certain experimental and numerical simulation methods. However, matching real formation conditions using experimental methods is often difficult, and traditional numerical simulation methods have some shortcomings, such as complex boundary processing and low computational efficiency. The lattice Boltzmann method (LBM) is a mesoscopic method of fluid dynamics derived from kinetic principles, whose microscopic nature is not limited by continuity assumptions. The LBM also has the mesoscopic characteristic of being recovered to the macroscopic model, and has many advantages such as a simple form and high parallel efficiency. Many scholars have researched mechanisms and models of acid-rock heterogeneous flow reactions, as well as flow reactions based on the LBM.
Labrid studied the heterogeneous flow reaction between a mud acid system and aluminosilicate at room temperature [4]. They discovered that the products generated by the heterogeneous flow reaction will react with hydrofluoric acid (HF) again, which will form a series of complex ions of aluminum fluoride and silicon fluoride. Gdanski traversed the entire heterogeneous flow reaction process between a mud acid system and sandstone minerals, including the first reaction between HF and aluminosilicate, the second reaction between H 2 SiF 6 and aluminosilicate, and the third reaction between HF and fluorosilicate precipitate [5]. Most of the mechanisms researched have studied the reaction between only a single mineral and an acid solution, and did not consider the influence of multicomponent minerals on the heterogeneous flow reaction of acid-rock. Rodoplu et al. established a model of the mechanism of wormhole formation in the process of acid-rock heterogeneous flow reactions based on laboratory experiments [6]. Li et al. established a mathematical model of acid-rock heterogeneous flow reaction based on the generalized distribution model, which can be applied to in situ matrix acidizing operations [7]. Most acidrock heterogeneous flow reaction models were macroscopic, but no matter how optimized, the macroscopic models were always limited by the assumption of continuous media and the difficulties of tracking phase interfaces. Many scholars have also studied the flow reaction in porous media by the LBM. Kang [10]. By combining the two, they studied the process of H + dissolving CaCO 3 under the influence of different temperatures. Tian et al. replaced the construction of complex porous media by adding a partial rebound term to the particle velocity evolution equation [11]. They also studied the heterogeneous flow reaction of CaCO 3 and H + in porous media containing both cracks and matrices.
To date, although some scholars have studied flow reactions with an LBM, few have studied multireaction processes between sandstone and acids. Our study focused on the development trend of "explaining macroscopic problems with mesoscopic and microcosmic models" and used the LBM as a special method to establish the corresponding heterogeneous flow reaction model of acid-rock. We solved the model using MATLAB and analyzed the distribution of various physical quantities after acidizing. We believe this study will be important for further guiding the optimization of acidizing process design parameters.

Establishing the Evolution Equation
In the LBM, an evolution equation can be used to describe the evolution process of a fluid particle velocity distribution function with discrete velocity on a fixed lattice. The evolution equation, also known as a lattice Boltzmann equation, originates from the Boltzmann equation. The Boltzmann equation is expressed as [12] ( ) The Boltzmann equation is discretized to obtain the evolution equation, which is expressed as The lattice Boltzmann equation consists of a collision step and a flow step. Through these two steps, the physical quantity of this time step can be transferred to the next time step.
Collision step: Streaming step: Ω i is a collision operator that indicates the effect of collisions between microscopic fluid particles in the velocity distribution function. Therefore, whether the collision operator is selected properly plays a decisive role for the LBM in accurately describing the physical laws of fluid flow.
Lattice ζ can be regarded as a closed system satisfying the conservation of mass, momentum, and internal energy. Through these three conservation equations, macroscopic physical quantities can be obtained:

Transformation from Lattice Space to Physical Space
In the LBM, the computational flow field is in a virtual lattice space, and all physical quantities are solved using dimensionless lattice units, which differ from an actual flow field in a physical space. To make LBM analysis meaningful, it is necessary to transform the lattice space into a physical space by finding the relations between lattice space parameters and physical space parameters. Currently, a more mature transformation method is dimensionless analysis.
In fluid mechanics theory, the computational flow field and the actual flow field need to meet certain similarity criteria, and the values of dimensionless characteristic numbers such as the Reynolds number (Re) remain unchanged in both flow fields: The Reynolds number [13] can be calculated according to the size l p , time t p , and kinematic viscosity ν p of the actual flow field. The dimensionless size and dimensionless time of the actual flow field are both 1. If the number of lattices in the computational flow field is set to N, then the dimensionless lattice step length △l p of the actual flow field is The length conversion parameter l r and the time conversion parameter t r are introduced to represent the length ratio and the time ratio respectively of the actual flow field to the computational flow field, which can be expressed as p r l l x δ ∆ = (8) and p r t t t δ ∆ = (9) In the LBM, square lattices are generally selected for calculation. The computational flow field lattice step length δx and time step length δt of a square lattice are both 1, thus the length conversion parameter l r can be calculated. In addition, the relation between the actual flow field velocity and the computational flow field velocity is In the dimensionless analysis, the actual flow field velocity u p is 1, from which the time conversion parameter t r can be calculated. According to the definition of a Reynolds number, the kinematic viscosity of the flow field can be calculated: The rest of the physical quantities can also be converted based on the dimensions, to achieve the conversion from lattice space to physical space.

Mathematical Model Based on the LBM
In the mathematical model based on the LBM, the Bhatnagar-Gross-Krook (BGK) model is most widely applied because it uses a simple linearized collision operator and therefore can accurately and efficiently reflect the motion law of fluid particles. Linearized collision operators are represented as The collision matrix and the equilibrium distribution function are two basic elements of the collision operator.
The single relaxation time model, also called the lattice Bhatnagar-Gross-Krook (LBGK) model, which uses a single relaxation time to determine the collision matrix, is currently the simplest and most widely used BGK model. In this study, the LBGK model was selected as the mathematical model. Its evolution equation is: The LBGK model can also be classified according to different equilibrium distribution functions. Among them, the n-dimensional and b-velocities (DnQb) model is the most representative; its equilibrium distribution function is: where c s is the lattice sound velocity, which can be calculated according to the lattice constant c: Compared with other models of the same dimension, the two-dimensional and nine-velocity (D2Q9) model has the most discrete velocity, so it also has better stability and accuracy. In this study, the D2Q9 model was used as the basic model for simulation. Figure 1 shows the D2Q9 model. In Eqs. 17 and 18, c i is the discrete velocity and ω is the weight coefficient. The discrete velocity and weight coefficient of the DnQb model can be calculated by the Gauss-Hermite integral.
The discrete velocity of the D2Q9 model is: The corresponding weight coefficient is The D2Q9 model follows the conservation of mass and momentum, and the macroscopic density and velocity in the model can be calculated by the following formulas: The Mach number (Ma) was introduced to indicate the ratio of the velocity at a point in the flow field to the local sound velocity at that point. The expression can be written as: When the Mach number is less than 0.3, the flow can be called a low Mach number flow. Obviously, the flow of acid in porous media is a low Mach number flow. In this case, through multiscale analysis, the LBGK model can be reduced to the macroscopic Navier-Stokes equation [14] by using the Chapman-Enskog expansion method in the kinetic theory: When the density of the fluid is close to a constant value, it can be regarded as a constant. Equation (23) is the standard incompressible Navier-Stokes equation. The kinematic viscosity ν and macroscopic pressure p in the model can be given by

Establishment of Heterogeneous Flow Reaction Model for Acid-Rock
The LBM is applied to the calculation of heat and mass transfer based on the LBM's mathematical model. In this study, a double-distribution-function (DDF) model was used to simulate the heat transfer process, and a convectiondispersion equation (CDE) model was used to simulate the mass transfer process. The heat transfer coefficients and the chemical reaction rates of various physical quantities were calculated according to the characteristics of an acid-rock heterogeneous flow reaction. We established an acid-rock heterogeneous flow reaction model combined with initial conditions and boundary conditions in porous media based on the LBM. Through simulation, the distribution of temperature, concentration of various substances, porosity, and permeability by acidizing were obtained and analyzed.

Double Distribution Function Model Simulates the Heat Transfer Process
There is a certain temperature difference between the acid solution and the formation matrix, and heat transfer will occur during the injection of the acid solution. The model for calculating heat transfer based on the LBM is called the thermal lattice Boltzmann model, which can be divided into a multispeed (MS) model and a DDF model.
The stability of an MS model is poor-generally applicable only to a situation where the Prandtl number is constant-and the MS model's application is also limited [15].
In the DDF model, both the temperature and seepage fields are calculated at the same velocity, which has several advantages, including a simple structure and a difficult divergence of calculation results.

DDF Model
The DDF model was proposed by Shan [16], who believed that temperature was a passive variable that changed with the movement of a fluid. Therefore, temperature can be regarded as an additional component in a fluid system. Also, a fluid system consisting of temperature and velocity is a two-component system. In this two-component system, each component is controlled by a similar equation and is operated at the same velocity.
A single relaxation time D2Q9 model was used to calculate the heat transfer, and the temperature evolution equation was: The temperature relaxation time τ T can be expressed as where K is the temperature transfer coefficient. The equilibrium distribution function is c s is the lattice sound velocity, which can be calculated according to the lattice constant c: The discrete velocity of the D2Q9 model is The corresponding weight coefficient is The model abides by the conservation of mass and momentum, and the macroscopic density and macroscopic velocity in the model can be calculated by The macroscopic temperature can be expressed as Through multiscale analysis, the DDF model can be returned to the correct heat transfer equation by using the Chapman-Enskog expansion method in the kinetic theory:

Temperature Field Model
Before using the LBM to establish a temperature field model, the following assumptions must be made about the real acid-rock heterogeneous flow reaction: The heat exchange between formations in the direction of well depth is ignored.
The specific volume heat, thermal conductivity, and density of the acid fluid and the formation rock do not change during the process of the acid fluid flow.
The heat transfer process is carried out uniformly. The acid solution flowing in the formation will cause a heat transfer with the formation rocks. However, the higher the temperature of the solution, the higher the reaction rate of the acid-rock heterogeneous flow, and the temperature will have a certain influence on the reaction. Arrhenius established the relation between the acid-rock reaction rate constant, the reaction activation energy, and the temperature. The reaction rate constant at a reservoir temperature can be obtained through the reaction rate constant measured in laboratory experiments: The heat transfer coefficient between acids and rocks can be calculated by The initial conditions of the temperature field are The boundary conditions of the temperature field are

Convection Dispersion Equation Model Simulates the Mass Transfer Process
During an acid-rock heterogeneous flow reaction, the concentration distribution of acid and minerals in porous media is not uniform, and mass transfer will occur. The model for calculating mass transfer based on the LBM is relatively simple, and currently only the CDE model is used to deal with the convection-diffusion equation.
The CDE model is based on the DDF model. In the same way as the temperature field can be seen as being driven by the movement of fluid, a change of the concentration field can be seen as being driven the same way. The calculation of the concentration and seepage fields in the CDE model uses the same velocity.

CDE Model
The CDE model can be subdivided into a single relaxation time model, a multiple relaxation time model, and an n-dimensional nonlinear model with a complex form [17][18][19].
In this study, the single relaxation D2Q9 model was still used for the mass transfer calculation, and the concentration evolution equation in the model is [20] The equilibrium distribution function is The calculation methods of seepage parameters such as discrete velocity, weight coefficient, macroscopic velocity, and macroscopic density are consistent with the standard generalized lattice Boltzmann equation (GLBE) model. Macroscopic concentration can be expressed as Zeng believed that, to distinguish whether the chemical reaction occurs at different locations in the same space, the evolution equation where the chemical reaction occurs in the same space must be added with the chemical reaction source term Q C , whereas the evolution equation where the chemical reaction did not occur did not need to be added with Q C [17].
The variable n s is the distribution density of solid particles in porous media, which can reflect the amount of solid phase and thus indirectly indicate the degree of chemical reaction. The introduction of n s can indicate not only the existence and strength of the chemical reaction to unify the evolution equation, but also the influence of the solid phase dissolution on the evolution equation.
For the mass transfer model of a chemical reaction, the chemical reaction concentration source terms Q C and n s are introduced into the evolution equation, and the evolution equation is rewritten as The chemical reaction concentration source term Q C can be expressed as By using the Chapman-Enskog expansion method in kinetic theory and through multiscale analysis, the model can be restored to obtain the correct convection-diffusion equation, namely The concentration relaxation time τ C can be expressed as

Concentration Field Model
Before using the LBM to establish a concentration field model, the following assumptions should be made about the real heterogeneous flow reaction of acid rocks: The influences of molecular diffusion and gravity are not considered.
The carbonate rocks in the reservoir have been treated with a prefluid before acidizing with a mud acid system, and the HCl in the mud acid system is not directly involved in the reaction during acidizing. The reservoir rocks can be divided into fast-reaction minerals and slow-reaction minerals according to the rate of heterogeneous flow reaction between the reservoir rocks and the mud acid system.
The complex acid-rock heterogeneous flow reaction is simplified by referring to the "two acids and three minerals" model. Feldspar and clay minerals that react with a mud acid system react faster and belong to fast-reaction minerals, whereas quartz reacts slowly and belongs to slow-reaction minerals.
By considering fast-reaction minerals and slow-reaction minerals as a whole, the reaction equation of them with a mud acid system can be unified as follows, respectively: The H 2 SiF 6 produced by the reaction between two minerals and the mud acid system will react with the fast-reaction minerals to form a Si(OH) 4 silica gel precipitate, which will cause secondary damage to the minerals. The chemical reaction equation can be written as 4  For the chemical reaction equations mentioned above, Da-Motta et al. [21] gave the values of the stoichiometric coefficient δ σ , as shown in Table 1. According to the molar concentration equilibrium principle of acid-rock heterogeneous flow reactions, the chemical reaction rates of HF, H 2 SiF 6 , fast-reaction minerals, slow-reaction minerals, and silica gel precipitate concentrations can be derived, and then the chemical reaction source terms of corresponding substances can be obtained.
For HF, the chemical reaction rate is ( ) For H 2 SiF 6 , the chemical reaction rate is For fast-reaction minerals, the chemical reaction rate is ( ) For slow-reaction minerals, the chemical reaction rate is ( ) For silica gel precipitate, the chemical reaction rate is The initial conditions of the concentration field are The boundary conditions of the concentration field are ( ) 10 0, 0, 0, ,

Porosity and Permeability Distribution Model
In the process of an acid-rock heterogeneous flow reaction, the concentration of rock minerals decreases, and the porosity of the formation increases after the acid solution reacts with the rock minerals.
After the concentration of the three rock minerals-feldspar, clay, and quartz-are calculated, the formation porosity can be worked out according to the principle of mineral volume balance: When minerals are dissolved in an acid solution, the permeability of the minerals will increase. Labrid proposed that the permeability with the same trend has the relation [4]

Model Solution and Analysis
In our simulation, the initial porosity was 20%, the initial permeability was 100×10 −3 µm 2 , the injection rate was 0.025m/s, and the dynamic viscosity was 1mPa·s. The values of the physical quantities used in the simulation were taken from Xue [22]. The specific data are shown in Table 2.

Temperature Distribution
We used the DDF model to simulate the temperature transfer of acid in porous media using the related physical quantities in Table 2. The temperature distribution is shown in Figure 2. The simulation results show that the effect of the temperature transfer between fluid and rock is minor and affects only the rock temperature at the entrance of the porous media. The distance of the temperature transfer is short.
In an acidizing operation, acid injection displacement and acid dosage are the two most important issues to be considered. Injection displacement can be reflected in the injection velocity, whereas acid dosage can be reflected in the injection time and velocity. Therefore, we analyzed the temperature distribution after the reaction by using different injection times and velocities. The resulting temperature distribution curves are shown in Figure 3.  The results indicate that the temperature transfer distance increases with increases in injection times and injection velocities. With the same acid dosage, the effect of a changing injection time on temperature transfer is more obvious than that of a changing injection velocity.

Concentration Distribution
Similarly, we simulated the mass transfer of an acid-rock heterogeneous flow reaction in porous media according to the relevant physical quantities in Table 2 using the CDE model. The concentration distribution of acids and minerals after the reaction was analyzed under different injection times and velocities. The corresponding concentration distribution curves are shown in Tables 3 and 4.

Concentration distribution curves at different injection velocities
Rapid-reactio n minerals Slow-reactio n minerals

Silica gel precipitate
The concentration distribution of HF shows a downward trend. With increases in injection times and velocities, the mass transfer distance increased gradually. The concentration distribution of H 2 SiF 6 increased first and then decreased.
Because most H 2 SiF 6 cannot react with formation minerals near an HF reaction front, the concentration of H 2 SiF 6 was the highest here, and the highest point gradually increased and moved to the right with increases in injection times and velocities. The concentration distribution of the fast-reaction minerals decreased first and then increased. Because those minerals continued to react with H 2 SiF 6 , the lowest point was near the front of the HF reaction and gradually decreased and moved to the right with increases in injection times and velocities. The concentration distribution of the slow-reaction minerals showed an upward trend. With increases in injection times and velocities, the mass transfer distance increased gradually. Silica gel precipitate and H 2 SiF 6 are reactants and generators of each other, and their concentration distribution was similar. They both increased first and then decreased. To some extent, the formation of the silica gel precipitate caused the pores that had been improved to be reblocked.
On the whole, HF concentration was the highest, and the concentration of other objects was low. With increases in injection times and velocities, the mass transfer distance of all objects increased. As it was with temperature transfer, the effect of changing injection times on mass transfer was more obvious than that of changing injection velocities.  To further verify the accuracy of the proposed model, we did a comparison analysis with the latest nonhomogeneous sandstone acidizing model [22]. The results are shown in Figure 4. From the concentration distribution curves of temperature and HF, it can be seen that the model in this study matches well with the nonhomogeneous sandstone acidizing model and has good accuracy.

Porosity and Permeability Distribution
After the reaction, the porosity and permeability at various locations of the porous media increased by different amounts. Controllable factors and uncontrollable formation factors can affect the increase of porosity and permeability after acidizing and also influence the acidizing effect.

Controllable Factors
By changing injection times and velocities, the increasing distribution curves of porosity and permeability respectively were obtained by simulation. The results are shown in Figure  5 and Figure 6. The fast-reaction minerals were consumed by H 2 SiF 6 near an HF reaction front, where the increase of porosity and permeability was the highest. With increases in injection times and velocities, the porosity and permeability increased gradually, and the influence of changing injection times on porosity and permeability was more obvious.
By changing the acid injection temperature and the HF injection concentration, we obtained the distribution curves of the porosity increases by simulation. The results are shown in Figure 7.
The results show that with increases in acid injection temperatures and HF injection concentrations, the porosity increased gradually. The effect on the porosity of changing the HF injection concentration is more obvious.
However, if the HF injection concentration is too high, the HF will excessively dissolve the reservoir minerals and easily cause pore collapse, blocking the formation. Therefore, the concentration must be analyzed in combination with the actual situation.

Uncontrollable Factors
By changing the formation heterogeneity (that is, changing the porosity deviation coefficient), the distribution curves of the porosity and permeability obtained by simulation increase. The results are shown in Figure 8. The results show that variations in the deviation coefficient had little effect on porosity and permeability. On the whole, porosity after acidizing was higher when the deviation coefficient was small, whereas permeability after acidizing was higher when the deviation coefficient was large.
If the content of fast-reaction minerals and slow-reaction minerals in the formation are changed, the distribution curves of porosity obtained by simulation increase. The results are shown in Figure 9.
The results indicate that the porosity increased gradually with an increase in reactive minerals in the formation. The effect of the fast-reaction mineral content on porosity was more obvious.
By comparing and analyzing the above factors, it can be found that the content of fast-reaction minerals, HF injection concentration, and acid injection time had a great influence on the acidizing effect, whereas the content of slow-reaction minerals, acid injection temperature, and injection velocity had little influence on the acidizing effect.
Excluding uncontrollable formation factors, the results suggest that to improve the acidizing effect, priority should be given to improve the HF concentration and acid dose.

Conclusions
Because the reaction mechanism of an acid-rock heterogeneous flow in porous media is very complex, it is often difficult to match real formation conditions by experimental methods, and traditional numerical simulation methods have the disadvantages of complex boundary treatments and inefficient calculations. We established a model of acid-rock heterogeneous flow reaction based on the LBM. The LBM has a micro nature that is not limited by the continuity assumption, and the mesoscopic characteristics that can be recovered to the macro model. The proposed model has the advantages of a simple form and high parallel efficiency.
The DDF model was used to simulate the heat transfer process based on the mathematical model of the LBM, and the CDE model was used to simulate the mass transfer process. The simulation results show that the effect of temperature transfer between fluid and rock was slight, and the distance of the temperature transfer was short. With an increase in injection time and injection velocity, the temperature transfer distance also increased. Compared with the injection velocity, the effect of the injection time on the temperature transfer was more obvious.
The concentration of HF was the highest, and the concentration of other substances was low. With an increase in injection time and velocity, the mass transfer distance of all substances increased. Similar to the temperature transfer, the effect on the mass transfer of changing the injection time was more obvious than the effect on the injection velocity.
Because the fast-reaction minerals in the middles of the channels were consumed by H 2 SiF 6 , the increasing rate of porosity and permeability in the middles of the channels was the highest. Fast-reaction mineral content, HF injection concentration, and acid injection time had a great influence on the acidizing effect, whereas slow-reaction mineral content, acid injection temperature, and injection velocity had little influence on the acidizing effect. We suggest that priority be given to improve HF concentration and acid dose to improve the acidizing effect.