A Multicomponent Thermal Fluid Numerical Simulation Method considering Formation Damage

Multicomponent thermal fluid huff and puff is an innovative heavy oil development technology for heavy oil reservoirs, which has been widely used in offshore oilfields in China and has proved to be a promising method for enhancing oil recovery. Components of multicomponent thermal fluids contain many components, including carbon dioxide, nitrogen, and steam. Under high temperature and high pressure conditions, the complex physical and chemical reactions between multicomponent thermal fluids and reservoir rocks occur, which damage the pore structure and permeability of core. In this paper, the authors set up a reservoir damage experimental device, tested the formation permeability before and after the injection of multiple-component thermal fluids, and obtained the formation damage model. The multicomponent thermal fluid formation damage model is embedded in the component control equation, the finite difference method is used to discretize the control equation, and a new multielement thermal fluid numerical simulator is established. The physical simulation experiment of multicomponent thermal fluid huff and puff is carried out by using the actual sand-packed model. By comparing the experimental results with the numerical simulation results, it is proved that the new numerical simulation model considering formation damage proposed in this paper is accurate and reliable.


Introduction
Multicomponent thermal fluid-enhanced oil recovery technology is based on the combustion injection principle of aerospace engine. The mixture of air and fuel is pressurized and injected into the multicomponent thermal fluid launcher. Then, start the igniter. The fuel combustion produces combustion products (N 2 , CO 2 ) at 2000°C. At the same time, water is injected into the multicomponent thermal fluid launcher to form nitrogen, carbon dioxide, and steam. The mixture of nitrogen, carbon dioxide, and steam is ejected at high speed through the nozzle and enters the reservoir. The schematic diagram of the multicomponent thermal fluid device is shown in Figure 1. Multicomponent thermal fluid technology does not require a smoke extraction process. It can fully utilize heat energy. Burning the same fuel, multicomponent thermal fluid technology can produce 16 percent more steam than conventional steam generators. Thermal efficiency of multicomponent thermal fluid technology can reach more than 95%. The main oil-increasing mechanism of multicomponent thermal fluid includes the dissolution viscosity reduction mechanism, heating viscosity reduction mechanism, supplementing formation energy mechanism, and increasing steam sweep efficiency mechanism. Due to the above characteristics of multicomponent thermal fluid, multicomponent thermal fluid technology has been widely used in the development of heavy oil reservoirs, and it has been widely used in Shengli oilfield and Bohai oilfield and achieved good oil-increasing results.
As a new emerging technology, multicomponent thermal fluid technology has achieved good application results, but the theoretical research on the seepage mechanism of multicomponent thermal fluid is not perfect.
Multicomponent thermal fluids contain many complex components, and temperature changes dramatically during seepage process, which brings great difficulties to numerical simulation. At present, many scholars have established numerical simulation methods for the particularity of multicomponent thermal fluids. Hou et al. [1] establish a two-zone thermal recovery seepage model for horizontal wells, which can calculate the productivity of horizontal wells with multicomponent thermal fluid huff and puff. Comparing with the results of CMG software, it is found that the calculation error is less than 5%. Dong et al. [2] take into account the phase and mass changes of fluid flowing in horizontal wellbore, couple fluid flow in reservoir with multistage fracturing wellbore flow in horizontal wellbore, and establish a coupled calculation model which can couple the flow in wellbore and reservoir. Sun et al. [3] establish heat conduction equation of multicomponent thermal fluids for concentric dualtubing wells and horizontal wells [4]. The new equation is solved by finite difference method. Type curves of SMTF flow in CDTW are obtained by finite difference method and iteration technique. Sun et al. studied the heat and mass transfer characteristics of superheated steam and noncondensable gas (SNG) flow in a horizontal wellbore [5]. In addition, Sun et al. conduct a series of works to study flow and heat transfer characteristics of superheated steam [6][7][8].
Based on the actual injection data of multithermal fluids in Nanpu Oilfield, Bohai Sea, China, Yang and Sun [9] established the equivalent numerical simulation model of multithermal fluids by historical fitting method on the basis of the conventional numerical simulation model of thermal recovery. Feng et al. [10] used the traditional thermal recovery numerical simulation model to analyze the sensitivity of the influencing factors of multicomponent thermal fluids to enhance oil recovery. It was considered that the oil saturation and oil viscosity were the important factors affecting the development effect of multicomponent thermal fluid technology. Zheng [11] used the traditional numerical simulation method of thermal recovery to analyze the law of gas channeling between wells in the process of multicomponent thermal fluid huff and puff, and proposed that multiwell multicomponent thermal fluid huff and puff technology and high-temperature foam plugging technology can effectively reduce the probability of gas channeling. Li et al. [12] used the traditional thermal numerical simulation method to analyze the influence of multicomponent thermal fluid injection parameters on enhanced oil recovery. He considered that the most sensitive parameters affecting multicomponent thermal fluid technology were injection-production ratio and injection temperature.
At present, the numerical simulation methods of multielement thermal fluids can be divided into two categories. One is an analytical method based on the theory of percolation mechanics, which is mainly used to calculate the temperature and pressure distribution in wellbore, but the prerequisite of the model is homogeneous model, which cannot be used in complex heterogeneous reservoirs. The second is the traditional thermal numerical simulation model. Although this model can be used to study the injection parameters of multicomponent thermal fluids in heterogeneous reservoirs, the traditional thermal recovery model cannot consider the damage of multicomponent thermal fluids to reservoir permeability [13][14][15].
In the process of injecting multicomponent thermal fluids to develop heavy oil, multicomponent thermal fluids react with reservoir rocks in a complex physical and chemical way, which will cause damage to reservoir permeability. At present, scholars have not carried out experiments on the influence of multicomponent thermal fluids on reservoir permeability. However, some scholars have carried out experiments on reservoir damage caused by steam. During steam huff and puff or steam drive, the pH value of steam reaches 8-9. After repeated huff and puff or long-time steam drive, great changes have taken place in pore structure and physical properties of reservoirs, which has a serious impact on the later development of oilfields [16]. There are two viewpoints on the influence of steam on reservoir permeability. Some scholars believe that steam reduces reservoir permeability. Through physical simulation experiments, it is proved that steam can swell clay and block pore throat passage [17][18][19]. Steam dissolves quartz particles and silicates in reservoirs, and mineral dissolution increases the content of SiO 2 , Al 2 O 3 , and K 2 O in the solution. These metal ions will produce zeolite, calcite, and calcite, and plug the pore throat [20]. In addition to reacting with reservoir rocks, steam also interacts with crude oil to produce asphalt precipitation, which ultimately reduces reservoir permeability [21][22][23]. Many scholars have established formation damage mathematical models according to the mechanism of steam damage to reservoirs. The main mathematical models describing 2 Geofluids steam damage to reservoirs are Wojtanowicz et al. model [24,25], Khilar and Fogler model [26], Gruesbeck and Collins model [27], Ohen and Civan model [28], Sharma and Yortsos model [29], and Rege and Fogler model [30]. These models can describe the deposition of incompressible solids on the pore surface and their migration and blockage. However, the variables in the above model are not easily obtained by experimental methods, such as the influence coefficient of rock particle migration on porosity, the influence coefficient of clay expansion on porosity, and the influence coefficient of solid particle migration and deposition on effective porosity. Moreover, the relationship between reservoir permeability and injection fluid temperature is not directly established in the above model nor does it directly describe the relationship between reservoir permeability and cumulative injection fluid PV number. Therefore, the above model cannot be directly used to describe the change of reservoir permeability after steam injection.
In order to establish a mathematical model that can calculate the permeability variation due to injecting multicomponent thermal fluid, a new type of test equipment has been developed. We conduct a series of reservoir damage experiments using multicomponent thermal fluid and analyze the law of multicomponent thermal fluid damage to reservoir permeability at different temperatures and different injection volumes. According to the component characteristics of multicomponent thermal fluids, a compositional numerical simulation model was established, and the formation damage model was embedded in the compositional numerical simulation model. The governing equation was discretized by finite difference method, and a new multicomponent thermal fluid numerical simulator was established. In order to verify the accuracy of the multicomponent thermal fluid numerical simulator established in this paper, the physical simulation experiment of multicomponent thermal fluid huff and puff is carried out by using the actual sand-packed model. Through comparison, it is found that the experimental results of physical simulation are in good agreement with the numerical simulation results.

Multicomponent Thermal Fluid Formation
Damage Experiment 2.1. Experimental Device. The experimental device is a hightemperature and high-pressure core displacement process. The main equipment includes an ISCO 260D high-pressure constant velocity pump, three intermediate container bottles, a steam generator, a constant temperature box, five pressure gauges, a high-temperature and high-pressure core holder, a DBR high-pressure back pressure valve, a condenser, and a production liquid metering system. The experimental device is shown in Figure 2.

Experimental Fluid and Experimental
Core. In the experiment, the volume of gas varies with the change of pressure. In the experimental scheme, the injection volume of CO 2 and N 2 is the volume under injection pressure. The steam injection rate is equivalent to the actual flow rate of multicomponent thermal fluid huff and puff in the field. The injection volume of steam and gas is based on the same heat quantity. Before the physical simulation, it is necessary to calculate the similarity criterion and scale the experimental conditions of physical simulation with the actual reservoir conditions in equal proportion. The principle of equal heat quantity is to ensure that the heat quantity of the multielement thermal fluid injected into the core is equal to the heat quantity injected into the reservoir by the actual injection well. The specific heat quantity of water varies obviously with temperature and pressure, which need to be determined according to injection conditions. The enthalpy values of steam (or hot water) at 140°C, 200°C, and 300°C are about 1020 kJ/kg, 158 kJ/kg, and 260 kJ/kg at the actual injection pressure. According to these enthalpy values, the injection amount and injection rate of steam, CO 2 , and N 2 are determined. The injection conditions of multicomponent thermal fluids are determined according to the relationship between the temperature of the multicomponent fluid generator and the composition of the fluid.
The scheme design of multicomponent thermal fluid reservoir damage experiment is shown in Table 1. The experimental core is the actual core from a Xinjiang oil field as shown in Figure 3. After the core is washed and dried, the basic parameters of the core are measured. The basic parameters of the core are shown in Table 1. The reservoir temperature is 80°C and the reservoir pressure is 25 MPa.

Experimental Steps
(1) Saturate the core with formation water. The formation water is used to test the water permeability of the core (2) The injection rate of multicomponent hot fluids is set to 1 mL/min. Adjust the temperature of the steam generator to be the same as the design temperature. Adjust the temperature of thermostat to be the same as that of reservoir. When the temperature of steam generator and thermostat meets the design requirements, multielement thermal fluids are injected. Record the pressure values during injection of multicomponent hot fluids  Figure 4, the initial permeability test results of the core and the permeability test results after injection of multicomponent thermal fluids are plotted. It can be seen from Figure 4 that the permeability values of four groups of cores decrease after injection of 24 PV multicomponent thermal fluids at different temperatures. It shows that the permeability of cores can be changed by multicomponent thermal fluids at different temperatures, which can damage   4 Geofluids the permeability of reservoirs. The higher the temperature of multicomponent thermal fluid is, the more obvious the damage of core permeability is. When the temperature of multicomponent thermal fluid is 573 K, the core permeability after injection of multicomponent thermal fluid is only 45% of the initial permeability. When the temperature of the multicomponent thermal fluid is 473 K, the core permeability after injection of the multicomponent thermal fluid is 58% of the initial permeability. When the temperature of multicomponent thermal fluid is 413 K, the core permeability after injection of multicomponent thermal fluid is 72% of the initial permeability. The damage to reservoir caused by multielement thermal fluids is not only related to the temperature of multielement thermal fluids but also to the injection amount of multielement thermal fluids. During the experiment, the core permeability with different volumes of injected multicomponent thermal fluids is obtained. The test results are shown in Figure 5. The permeability value of core decreases with the volume of injected multicomponent thermal fluid increasing. When the volume of multielement thermal fluid injection is less than 6 PV, core permeability decreases rapidly with the increase of multielement thermal fluid injection; when the volume of multielement thermal fluid injection is larger than 6 PV, core permeability changes little with the volume of multicomponent thermal fluid injection, but there is a small fluctuation up and down. It shows that the damage of multicomponent thermal fluid to core permeability is not instantaneous, and the damage of multicomponent thermal fluid to core needs a certain action time. At the same time, it shows that the multicomponent thermal fluid destroys the original cements, makes the particles migrate in the core, and makes the particles migrate to different positions, which will have a great or small impact on the core permeability. Therefore, after the injection of 6 PV multicomponent thermal fluid, the core permeability will fluctuate slightly.
In order to quantify the damage of multicomponent thermal fluids to core permeability, according to the experimental results, a multicomponent nonlinear regression algorithm is used to obtain the formation damage model of multicomponent thermal fluid reservoir permeability, as shown in Equation (1). The values of the parameters in Equation (1) are shown in Table 2.
In order to quantify the damage of multicomponent thermal fluids to core permeability, a formation damage model of multicomponent thermal fluids is obtained by using multivariate nonlinear regression algorithm based on experimental results. The formation damage model is shown in Equation (1). The applicable conditions of Equation (1) are that the displacing fluid is multicomponent thermal fluid, and the porous media is sandstone with permeability of sandstone that ranges 500 md to 1500 md. The values of undetermined coefficients in Equation (1) are shown in Table 2. In order to verify the fitting accuracy of Equation (1), the calculation results of Equation (1) and the experimental results are drawn in Figure 5. From Figure 5, it can be seen that Equation (1) is basically consistent with the experimental results.

Multicomponent Thermal Fluid Compositional Numerical Simulation Model considering Formation Damage
The compositional numerical simulation method of multicomponent thermal fluid must consider the temperature change, phase state change, the interaction between multicomponent According to the principle of mass conservation, the equation of mass conservation of water components is as follows:  (1)).
The mass conservation equation of gas components is as follows: The mass conservation equation of light oil components is as follows: The mass conservation equation of heavy oil components is as follows: Considering the damage of multicomponent thermal fluids to formation permeability, absolute permeability is no longer constant in the compositional numerical simulation model. In this paper, an explicit time-varying permeability method is used to calculate the absolute permeability. According to the absolute permeability, PV number, and temperature of the T − 1 time step, the absolute permeability of the T time step is calculated. The specific calculation method is shown in the following equation: K t ð Þ = beta1 × PV 5 + beta2 × PV 4 + beta3 × PV 3 + beta4 À × PV 2 + beta5 × PV + beta6 × T 3 + beta7 × T 2 In the process of numerical simulation of multivariate thermal fluids, the change of formation temperature must be taken into account, and the energy conservation equation should be established. The energy conservation equation is shown in the following equation: According to the law of conservation of energy, the equation of conservation of energy for multicomponent thermal fluids can be obtained, as shown in the following equation: Constraint equations are classified into four main categories. They are saturation constraint equation, mole fraction constraint equation, phase equilibrium equation, and capillary force equation. Among them, the saturation constraint equation is as follows: The molar fraction constraint equations are as follows: The phase equilibrium equations are as follows: The capillary force equations are as follows: When the production rate is constant, the internal boundary condition is expressed as follows: When the bottom hole pressure is constant, the internal boundary condition is expressed as follows: Initial conditions are as follows: 3.2. Differential Discrete Method for Compositional Numerical Simulation Model of Multicomponent Thermal Fluids. In order to discretize the governing equation, the numerical simulation model is divided into N x × N y × N z uniform grids. Let the space step of the grid (i, j, k) be Δx, Δy, Δz and the time step be Δt. The superscript n denotes t = n ⋅ Δt moment. The thermal conductivity of the fluid is taken as the harmonic average of two adjacent grids, as shown in the following equations: The velocity of fluid in the grid is expressed by the following equations: In the formula, the ratio of absolute permeability to viscosity is equal to the harmonic average of two adjacent grids.
Relative permeability takes the grid value in the direction of flow. The calculation method of relative permeability is shown in the following equations: In summary, the governing equations can be discretized by difference method. The discrete form of the mass conservation equation for water components is as follows: The discrete form of mass conservation equation for dissolved gas components is as follows: 8 Geofluids The difference discrete form of mass conservation equation for light oil components is as follows: The discrete form of mass conservation equation for heavy oil components is as follows: The discrete form of the energy conservation equation is as follows:

Solution Method of Multicomponent Thermal Fluid
Numerical Simulation Model. Each grid node contains four discrete governing equations. Each discrete governing equation has strong nonlinearity. In order to improve computational stability and speed up convergence, the Newton-Raphson iteration method is used to linearize the nonlinear equation. The discrete equation of four components can be written as follows: In Equation (34), vector F represents the discrete equation of four components, and vector X represents the variable to be desired. By using the Newton-Raphson  Figure 6: The flowchart of the multicomponent thermal fluid high-temperature and high-pressure experimental device.
Equation (35) can be expressed as follows: In Equation (36), J is the Jacobian determinant and can be calculated by the following equation: Let F m ðX + δXÞ = 0 and ignore oðδX 2 Þ, then you can get the following: Determinant J is a large sparse matrix, which is solved by conjugate gradient method. After the solution δX is obtained, the convergence is judged. If it is not convergent, the value of the next new iteration step is obtained by taking δX as the increment of the initial value of the previous iteration.

Verification of Multicomponent Thermal Liquid Compositional Numerical Simulation Model
In order to simulate and study the huff and puff and displacement effect of multicomponent thermal fluids under high pressure, we specially constructed a multicomponent thermal fluid high-pressure experimental device with four temperature measuring points and two pressure measuring points.
The flowchart of the multicomponent thermal fluid hightemperature and high-pressure experimental device is shown in Figure 6.
11 Geofluids intermediate vessel, oil, gas, and water metering device, metering pump, flowmeter, pressure meter, and data acquisition and control system. In order to prevent sand production, it is necessary to install screens at the upper and lower ends of the sand-packed model.
Experimental steps are shown as follows: (1) Clean the actual oil sand   (4) Saturate the model with crude oil until the produced liquid does not contain water, and the liquid production speed is stable, and place the model for 12 hours (5) Set back pressure to 20 MPa. Injection of 350 mL steam + N 2 /CO 2 from the production end recorded changes in injection rate and pressure (6) After injection of 350 mL multicomponent thermal fluid, close the well for more than 2 hours according to the design time (7) Constant injection pressure at injection end (injection oil) and stable back pressure at production end, oil production rate, and water production rate are recorded until the oil production rate is 28 mL/(MPa min) (equivalent to the speed of cold production) (8) Return pressure to empty, record oil production rate, water production rate and gas production rate, and pressure According to the geometrical size (diameter 2.5 cm and length 40 cm) of the physical simulation sand-packed model, a three-dimensional numerical simulation model is established. The basic parameters of the numerical simulation model are shown in Table 3. Mesh of numerical simulation model is shown in Figure 7. The relative permeability curve used in the numerical simulation model is shown in Figures 8 and 9.
In this paper, we set the injection parameters and production parameters of numerical simulation of multicomponent thermal fluid huff and puff according to the process of physical simulation experiment. The oil production index, gas production index, and water production index are calculated by using the numerical simulation model considering reservoir damage and the numerical simulation model without considering reservoir damage, respectively. The results of physical simulation experiment and numerical simulation are shown in Figure 10.
As can be seen from Figure 10, the oil production index, water production index, and gas production index all decrease with the increase of cumulative oil production. Firstly, the traditional numerical simulation method (without considering the damage of multicomponent thermal fluid to formation permeability) is used to calculate the multicomponent thermal fluid huff and puff process. The calculation results show that the oil production index, water production index, and gas production index calculated by the traditional numerical simulation method are significantly higher than the experimental results. In particular, the oil production index calculated by traditional numerical simulation method is twice the oil production index obtained by laboratory test. Keep all the basic parameters of the numerical simulation model unchanged. A new numerical simulation method (considering the damage of multicomponent thermal fluids to formation permeability) presented in this paper is used to simulate the multicomponent thermal fluid huff and puff process. The calculation results show that the results of the new numerical simulation model (considering the damage of multicomponent thermal fluids to formation permeability) are closer to the experimental results. There are three main reasons for the decline of oil production index, water production index, and gas production index. The first reason is the change of saturation. According to relative permeability curves (as shown in Figures 8 and 9), the change of water saturation, oil saturation, and gas saturation will lead to the change of water phase effective permeability, oil phase effective permeability, and gas phase effective permeability. The second reason is the formation pressure drop. Because the production process of multicomponent thermal fluid huff and puff is a depletion development process, the decline of formation pressure will lead to the decline of oil production, water production, and gas production. The third reason is that the pore structure of the reservoir is damaged by multicomponent thermal fluids, and the permeability of the reservoir decreases. Traditional numerical simulation method can simulate the influence of saturation and pressure on production index, but it cannot simulate the influence of reservoir permeability damage on production index. Traditional numerical simulation methods can simulate the effects of saturation and pressure changes on production indicators, but cannot simulate the effects of formation damage on production indicators. Based on the traditional numerical simulation method, the formula of multicomponent thermal fluid formation damage is introduced into the compositional numerical simulation model. The results of compositional numerical simulation model considering reservoir damage are very close to the experimental results, which shows that the proposed compositional numerical simulation method considering formation damage is accurate and reliable.

Conclusion
(1) A multielement thermal fluid reservoir damage experimental device was established to test core permeability before and after multielement thermal fluid injection at different temperatures and different injection volumes. It was found that with the increase of multielement thermal fluid temperature, the damage degree of core permeability increased (2) The damage of multicomponent thermal fluid to core permeability needs a certain action time. When the injection volume of multicomponent thermal fluid is less than 6PV, the core permeability decreases rapidly with the increase of the injection volume of multicomponent thermal fluid. When the injection volume of multicomponent thermal fluid is larger than 6PV, the core permeability changes little with the increase of the injection volume of multicomponent thermal fluid Quality of outflow unit (kg) M 3 : Quality of source/sink (kg) o, w, g: Oil, water, gas ϕ: Porosity (dimensionless) KðtÞ: Absolute permeability at t-step (m 2 ) t: Time (s) K rg , K ro , K rw : Relative permeability of gas, oil, and water phases (dimensionless) ρ g , ρ o , ρ w : Molar density of gas, oil, and water phases (mol/m 3 ) μ g , μ o , μ w : Viscosity of gas, oil, and water phases (Pa·s) P g , P o , P w : Pressure of gas phase, oil phase, and water phase (Pa) q j g , q j o , q j w : Volume flow rate of gas, oil, and water phases in well no. J (m 3 /s) S g , S o , S w : Saturation of gas, oil, and water phases (dimensionless) ξ: Gravity acceleration (m/s 2 ) Z: Reservoir depth (m) Temperature (°C) λ: Thermal conductivity of fluid (W/(m·°C)) U g , U o , U w : Internal energy of gas, oil, and water phases (J/mol) ρ R : Molar density of rocks (mol/m 3 ) U R : Internal energy of rock (J/mol) K 1 : Phase equilibrium constants of dissolved gas components (dimensionless) K 2 : Phase equilibrium constants of light oil components (dimensionless) P sat : Saturated vapor pressure of steam (MPa) P cow : Oil-water phase capillary force P cog : Oil-gas phase capillary force.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that they have no conflicts of interest.