# A Multicomponent Thermal Fluid Numerical Simulation Method considering Formation Damage

Authors: Xinan Yu; Xiaoping Li; Shuoliang Wang; Yi Luo
Journal: Geofluids
 (2021)

---

## Abstract

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.


---

## Body

## 1. 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 (N2, CO2) 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.Figure 1
Schematic diagram of multicomponent thermal fluid generator.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 dual-tubing 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–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–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–19]. Steam dissolves quartz particles and silicates in reservoirs, and mineral dissolution increases the content of SiO2, Al2O3, and K2O 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–23]. Many scholars have established formation damage mathematical models according to the mechanism of steam damage to reservoirs. The main mathematical models describing 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.


## 2. Multicomponent Thermal Fluid Formation Damage Experiment

### 2.1. Experimental Device

The experimental device is a high-temperature 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 Figure2.Figure 2
Multicomponent thermal fluid formation damage experiment.


### 2.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 CO2 and N2 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, CO2, and N2 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 Table1. 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.Table 1
Multicomponent thermal fluid formation damage experimental scheme.
ParameterFluidsNo. 1 core groupNo. 2 core groupNo. 3 core groupNo. 4 core groupPermeability (10-3 μm2)125612471290934926957703716745523538540Porosity (%)35.333.934.730.229.531.529.530.629.628.527.325.6Steam temperature (K)413473543423483553433493563453523573Velocity (mL/min)Steam0.70.450.970.640.430.910.600.410.790.500.40.4CO20.0320.0320.0320.0320.0320.0320.0320.0320.0320.0320.0320.032N20.150.150.150.150.150.150.150.150.150.150.150.15Figure 3
Pictures of experimental core.


### 2.3. 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(3)
When the volume of multicomponent thermal fluid reaches the PV number required by the scheme design (2 PV, 4 PV, 6 PV, 8 PV, 10 PV, 12 PV, 14 PV, 16 PV, 18 PV, 20 PV, 22 PV, and 24 PV), the formation water is injected and the permeability of the core is tested for the second time(4)
According to the experimental scheme, the temperature of multicomponent thermal fluid is changed and the above experimental steps are repeated. The temperatures of multicomponent thermal fluids are 413 K, 473 K, 543 K, 423 K, 483 K, 553 K, 433 K, 493 K, 563 K, 453 K, 523 K, and 573 K, respectively


### 2.4. Experimental Results

In Figure4, 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 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.Figure 4
Experimental results of formation damage experiments of multicomponent thermal fluids at different temperatures.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 Figure5. 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.Figure 5
Experimental results of formation damage experiments of multicomponent thermal fluids at different injection volume of multicomponent thermal fluids (the scatter point is the experimental results and the solid line is the results calculated by Equation (1)).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.Table 2
Parameter table of formation damage model.
ParameterValueParameterValueParameterValuebeta1−7.199799E−07beta25.652760E−05beta3−1.709992E−03beta42.478296E−02beta5−1.711991E−01beta6−7.771773E−08beta75.550085E−05beta8−1.430857E−02beta92.29120371In 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.
(1)KKi=beta1×PV5+beta2×PV4+beta3×PV3+beta4×PV2+beta5×PV+beta6×T3+beta7×T2+beta8×T+beta9.


## 2.1. Experimental Device

The experimental device is a high-temperature 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 Figure2.Figure 2
Multicomponent thermal fluid formation damage experiment.


## 2.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 CO2 and N2 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, CO2, and N2 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 Table1. 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.Table 1
Multicomponent thermal fluid formation damage experimental scheme.
ParameterFluidsNo. 1 core groupNo. 2 core groupNo. 3 core groupNo. 4 core groupPermeability (10-3 μm2)125612471290934926957703716745523538540Porosity (%)35.333.934.730.229.531.529.530.629.628.527.325.6Steam temperature (K)413473543423483553433493563453523573Velocity (mL/min)Steam0.70.450.970.640.430.910.600.410.790.500.40.4CO20.0320.0320.0320.0320.0320.0320.0320.0320.0320.0320.0320.032N20.150.150.150.150.150.150.150.150.150.150.150.15Figure 3
Pictures of experimental core.


## 2.3. 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(3)
When the volume of multicomponent thermal fluid reaches the PV number required by the scheme design (2 PV, 4 PV, 6 PV, 8 PV, 10 PV, 12 PV, 14 PV, 16 PV, 18 PV, 20 PV, 22 PV, and 24 PV), the formation water is injected and the permeability of the core is tested for the second time(4)
According to the experimental scheme, the temperature of multicomponent thermal fluid is changed and the above experimental steps are repeated. The temperatures of multicomponent thermal fluids are 413 K, 473 K, 543 K, 423 K, 483 K, 553 K, 433 K, 493 K, 563 K, 453 K, 523 K, and 573 K, respectively


## 2.4. Experimental Results

In Figure4, 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 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.Figure 4
Experimental results of formation damage experiments of multicomponent thermal fluids at different temperatures.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 Figure5. 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.Figure 5
Experimental results of formation damage experiments of multicomponent thermal fluids at different injection volume of multicomponent thermal fluids (the scatter point is the experimental results and the solid line is the results calculated by Equation (1)).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.Table 2
Parameter table of formation damage model.
ParameterValueParameterValueParameterValuebeta1−7.199799E−07beta25.652760E−05beta3−1.709992E−03beta42.478296E−02beta5−1.711991E−01beta6−7.771773E−08beta75.550085E−05beta8−1.430857E−02beta92.29120371In 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.
(1)KKi=beta1×PV5+beta2×PV4+beta3×PV3+beta4×PV2+beta5×PV+beta6×T3+beta7×T2+beta8×T+beta9.


## 3. 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 thermal fluid and reservoir rock, and the interaction between multicomponent thermal fluid and reservoir fluid. Considering the characteristics of multicomponent thermal fluids, the following basic assumptions are made for the model:(1)
There are oil, water, and gas three-phase fluids in the fluid(2)
Fluid flow conforms to Darcy’s law(3)
There is no chemical reaction in the process of fluid flow(4)
The fluid consists of three phases (oil phase, water phase, and gas phase) and four components (light oil component, heavy oil component, water component, and gas component). The phase change of components obeys the principle of phase equilibrium(5)
Fluid seepage is an isothermal seepage process(6)
The way of heat transfer is convection and heat transfer


### 3.1. Basic Governing Equations for Compositional Numerical Simulation of Multicomponent Thermal Fluids

The compositional numerical simulation model consists of four kinds of equations: mass conservation equation, energy conservation equation, constraint equation, and boundary condition equation. Among them, the mass conservation equation is based on the principle of mass conservation. For any component, the principle of mass conservation can be expressed as follows:(2)ΔM=M1−M2+M3.According to the principle of mass conservation, the equation of mass conservation of water components is as follows:(3)∂∂tϕρwSw+ρgy1Sg+∇⋅ρwKt⋅krwμw∇Pw−ρwξ∇Z+ρgy1Kt⋅krgμg∇Pg−ρgξ∇Z=∑j=1Nρwqwj+ρgy1qgj.The mass conservation equation of gas components is as follows:(4)∂∂tϕρox1So+ρoy2Sg+∇⋅ρox1Kt⋅kroμo∇Po−ρoξ∇Z+ρgy2Kt⋅krgμg∇Pg−ρgξ∇Z=∑j=1Nρox1qoj+ρgy2qgj.The mass conservation equation of light oil components is as follows:(5)∂∂tϕρox2So+ρoy3Sg+∇⋅ρox2Kt⋅kroμo∇Po−ρoξ∇Z+ρgy3Kt⋅krgμg∇Pg−ρgξ∇Z=∑j=1Nρox2qoj+ρgy3qgj.The mass conservation equation of heavy oil components is as follows:(6)∂∂tϕρox3So+∇⋅ρox3K⋅kroμo∇po−ρoξ∇Z=∑j=1Nρox2qoj.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 theT−1 time step, the absolute permeability of the T time step is calculated. The specific calculation method is shown in the following equation:
(7)Kt=beta1×PV5+beta2×PV4+beta3×PV3+beta4×PV2+beta5×PV+beta6×T3+beta7×T2+beta8×T+beta9×Kt−1.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:(8)ΔQ=Q1+Q2+Q3−Q4.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:(9)∂∂tϕρwSwUw+ρoSoUo+ρgSgUg+1−ϕρRUR+∇⋅ρwUwVw+ρoUoVo+ρgUgVg=∇⋅λ∇T+∑j=1NρwUwqwj+ρoUoqoj+ρgUgqgj.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:(10)Sw+So+Sg=1.The molar fraction constraint equations are as follows:(11)x1+x2+x3=1,(12)y1+y2+y3=1.The phase equilibrium equations are as follows:(13)K1p,T=y2x1,(14)K2p,T=y3x2,(15)y1=PsatPg.The capillary force equations are as follows:(16)Pw=Po−PcowSw,(17)Pg=Po−PcogSg.When the production rate is constant, the internal boundary condition is expressed as follows:(18)qI=const.When the bottom hole pressure is constant, the internal boundary condition is expressed as follows:(19)pwf=const.Initial conditions are as follows:(20)P=Pit=0,T=Tit=0.


### 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 intoNx×Ny×Nz 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:(21)λi±1/2,j,kn+1=2λi,j,kn+1λi±1,j,kn+1λi,j,kn+1±λi±1,j,kn+1,λi,j±1/2,kn+1=2λi,j,kn+1λi,j±1,kn+1λi,j,kn+1±λi,j±1,kn+1,λi,j,k±1/2n+1=2λi,j,kn+1λi,j,k±1n+1λi,j,kn+1±λi,j,k±1n+1.The velocity of fluid in the grid is expressed by the following equations:(22)Vi±1/2,j,kn+1=±Kxμi±1/2,j,kn+1kr,i±1/2,j,kn+1Pi,j,kn+1−Pi±1,j,kn+1Δx,(23)Vi,j±1/2,kn+1=±Kxμi,j±1/2,kn+1kr,i,j±1/2,kn+1Pi,j,kn+1−Pi,j±1,kn+1Δy,(24)Vi,j,k±1/2n+1=±Kxμi,j,k±1/2n+1kr,i,j,k±1/2n+1Pi,j,kn+1−Pi,j,k±1n+1Δz.In the formula, the ratio of absolute permeability to viscosity is equal to the harmonic average of two adjacent grids.(25)Kxμi±1/2,j,kn+1=2Kx/μi±1,j,kn+1Kx/μi,j,kn+1Kx/μi±1,j,kn+1+Kx/μi,j,kn+1,Kyμi,j±1/2,kn+1=2Ky/μi,j±1,kn+1Ky/μi,j,kn+1Ky/μi,j±1,kn+1+Ky/μi,j,kn+1,Kzμi,j,k±1/2n+1=2Kz/μi,j,k±1n+1Kz/μi,j,kn+1Kz/μi,j,k±1n+1+Kz/μi,j,kn+1.Relative permeability takes the grid value in the direction of flow. The calculation method of relative permeability is shown in the following equations:(26)kr,i+1/2,j,kn+1=kr,i,j,kn+1Vi+1/2,j,kn+1≥0,kr,i+1,j,kn+1Vi+1/2,j,kn+1<0,kr,i−1/2,j,kn+1=kr,i−1,j,kn+1Vi−1/2,j,kn+1≥0,kr,i,j,kn+1Vi−1/2,j,kn+1<0,(27)kr,i,j+1/2,kn+1=kr,i,j,kn+1Vi,j+1/2,kn+1≥0,kr,i,j+1,kn+1Vi,j+1/2,kn+1<0,kr,i,j−1/2,kn+1=kr,i,j−1,kn+1Vi,j−1/2,kn+1≥0,kr,i,j,kn+1Vi,j−1/2,kn+1<0,(28)kr,i,j,k+1/2n+1=kr,i,j,kn+1Vi,j,k+1/2n+1≥0,kr,i,j,k+1n+1Vi,j,k+1/2n+1<0,kr,i,j,k−1/2n+1=kr,i,j,k−1n+1Vi,j,k−1/2n+1≥0,kr,i,j,kn+1Vi,j,k−1/2n+1<0.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:(29)ϕi,j,kn+1ρwSwi,j,kn+1−ϕi,j,knρwSwi,j,kn+ϕi,j,kn+1ρgy1Sgi,j,kn+1−ϕi,j,knρgy1Sgi,j,knΔxΔyΔz+ρwVwxi+1/2,j,kn+1−ρwVwxi−1/2,j,kn+1+ρgy1Vgxi+1/2,j,kn+1−ρgy1Vgxi−1/2,j,kn+1ΔyΔzΔt+ρwVwyi,j+1/2,kn+1−ρwVwyi,j−1/2,kn+1+ρgy1Vgyi,j+1/2,kn+1−ρgy1Vgyi,j−1/2,kn+1ΔxΔzΔt+ρwVwzi,j,k+1/2n+1−ρwVwzi,j,k−1/2n+1+ρgy1Vgzi,j,k+1/2n+1−ρgy1Vgzi,j,k−1/2n+1ΔxΔyΔt=∑j=1Nρwqwj+ρgy1qgjn+1Δt.The discrete form of mass conservation equation for dissolved gas components is as follows:(30)ϕi,j,kn+1ρox1Soi,j,kn+1−ϕi,j,knρox1Soi,j,kn+ϕi,j,kn+1ρgy2Sgi,j,kn+1−ϕi,j,knρgy2Sgi,j,knΔxΔyΔz+ρox1Voxi+1/2,j,kn+1−ρox1Voxi−1/2,j,kn+1+ρgy2Vgxi+1/2,j,kn+1−ρgy2Vgxi−1/2,j,kn+1ΔyΔzΔt+ρox1Voyi,j+1/2,kn+1−ρox1Voyi,j−1/2,kn+1+ρgy2Vgyi,j+1/2,kn+1−ρgy2Vgyi,j−1/2,kn+1ΔxΔzΔt+ρox1Vozi,j,k+1/2n+1−ρox1Vozi,j,k−1/2n+1+ρgy2Vgzi,j,k+1/2n+1−ρgy2Vgzi,j,k−1/2n+1ΔxΔyΔt=∑j=1Nρox1qoj+ρgy2qgjn+1Δt.The difference discrete form of mass conservation equation for light oil components is as follows:(31)ϕi,j,kn+1ρox2Soi,j,kn+1−ϕi,j,knρox2Soi,j,kn+ϕi,j,kn+1ρgy3Sgi,j,kn+1−ϕi,j,knρgy3Sgi,j,knΔxΔyΔz+ρox2Voxi+1/2,j,kn+1−ρox2Voxi−1/2,j,kn+1+ρgy3Vgxi+1/2,j,kn+1−ρgy3Vgxi−1/2,j,kn+1ΔyΔzΔt+ρox2Voyi,j+1/2,kn+1−ρox2Voyi,j−1/2,kn+1+ρgy3Vgyi,j+1/2,kn+1−ρgy3Vgyi,j−1/2,kn+1ΔxΔzΔt+ρox2Vozi,j,k+1/2n+1−ρox2Vozi,j,k−1/2n+1+ρgy3Vgzi,j,k+1/2n+1−ρgy3Vgzi,j,k−1/2n+1ΔxΔyΔt=∑j=1Nρox2qoj+ρgy3qgjn+1Δt.The discrete form of mass conservation equation for heavy oil components is as follows:(32)ϕi,j,kn+1ρox3Soi,j,kn+1−ϕi,j,knρox3Soi,j,knΔxΔyΔz+ρox3Vxoi+1/2,j,kn+1−ρox3Vxoi−1/2,j,kn+1ΔyΔzΔt+ρox3Vyoi,j+1/2,kn+1−ρox3Vyoi,j−1/2,kn+1ΔxΔzΔt+ρox3Vzoi,j,k+1/2n+1−ρox3Vzoi,j,k−1/2n+1ΔxΔyΔt=∑j=1Nρox3qojn+1Δt.The discrete form of the energy conservation equation is as follows:(33)ϕi,j,kn+1ρwSwUwi,j,kn+1+ρgSgUgi,j,kn+1+ρoSoUoi,j,kn+1−ϕi,j,kn+1ρwSwUwi,j,kn+ρgSgUgi,j,kn+ρoSoUoi,j,kn+1−ϕi,j,kn+1ρRURi,j,kn+1−1−ϕi,j,kn+1ρRURi,j,knΔxΔyΔz+ρwVxwUwi+1/2,j,kn+1−ρwVxwUwi−1/2,j,kn+1+ρoVxoUoi+1/2,j,kn+1−ρoVxoUoi−1/2,j,kn+1+ρgVxgUgi+1/2,j,kn+1−ρgVxgUgi−1/2,j,kn+1ΔyΔzΔt+ρwVywUwi,j+1/2,kn+1−ρwVywUwi,j−1/2,kn+1+ρoVyoUoi,j+1/2,kn+1−ρoVyoUoi,j−1/2,kn+1+ρgVygUgi,j+1/2,kn+1−ρgVygUgi,j−1/2,kn+1ΔxΔzΔt+ρwVzwUwi,j,k+1/2n+1−ρwVzwUwi,j,k−1/2n+1+ρoVzoUoi,j,k+1/2n+1−ρoVzoUoi,j,k−1/2n+1+ρgVzgUgi,j,k+1/2n+1−ρgVzgUgi,j,k−1/2n+1ΔxΔyΔt=λi+1/2,j,kn+1Ti+1,j,kn+1−Ti,j,kn+1Δx−λi−1/2,j,kn+1Ti,j,kn+1−Ti−1,j,kn+1ΔxΔyΔzΔt+λi,j+1/2,kn+1Ti,j+1,kn+1−Ti,j,kn+1Δy−λi,j−1/2,kn+1Ti,j,kn+1−Ti,j−1,kn+1ΔyΔxΔzΔt+λi,j,k+1/2n+1Ti,j,k+1n+1−Ti,j,kn+1Δz−λi,j,k−1/2n+1Ti,j,kn+1−Ti,j,k−1n+1ΔzΔxΔyΔt+∑j=1NρwqwjUwn+1+∑j=1NρoqojUon+1+∑j=1NρgqgjUgn+1Δt.


### 3.3. 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:(34)FX=0.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 iteration method, Taylor expansion of the m-th equation Fm can be obtained:
(35)FmX+δX=FmX+∑l=14×Nx×Ny×Nz∂Fm∂xlδxl+oδX2.Equation (35) can be expressed as follows:
(36)FmX+δX=FmX+J⋅δX+oδX2.In Equation (36), J is the Jacobian determinant and can be calculated by the following equation:
(37)Jml=∂Fm∂xl=Fmx1,x2,⋯,xl+Δxl,⋯,x4×Nx×Ny×Nz−Fmx1,x2,⋯,xl,⋯,x4×Nx×Ny×NzΔxl.LetFmX+δX=0 and ignore oδX2, then you can get the following:
(38)J⋅δX=−F.DeterminantJ 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.
(39)Xnew=Xold+δX.Continuous iteration until the following convergence conditions are satisfied:(40)14×Nx×Ny×Nz∑l=14×Nx×Ny×NzδxLxl<ε.Through the above steps, the solutions of the discrete equations (Equations (29), (30), (31), (32), and (33)) can be obtained.


## 3.1. Basic Governing Equations for Compositional Numerical Simulation of Multicomponent Thermal Fluids

The compositional numerical simulation model consists of four kinds of equations: mass conservation equation, energy conservation equation, constraint equation, and boundary condition equation. Among them, the mass conservation equation is based on the principle of mass conservation. For any component, the principle of mass conservation can be expressed as follows:(2)ΔM=M1−M2+M3.According to the principle of mass conservation, the equation of mass conservation of water components is as follows:(3)∂∂tϕρwSw+ρgy1Sg+∇⋅ρwKt⋅krwμw∇Pw−ρwξ∇Z+ρgy1Kt⋅krgμg∇Pg−ρgξ∇Z=∑j=1Nρwqwj+ρgy1qgj.The mass conservation equation of gas components is as follows:(4)∂∂tϕρox1So+ρoy2Sg+∇⋅ρox1Kt⋅kroμo∇Po−ρoξ∇Z+ρgy2Kt⋅krgμg∇Pg−ρgξ∇Z=∑j=1Nρox1qoj+ρgy2qgj.The mass conservation equation of light oil components is as follows:(5)∂∂tϕρox2So+ρoy3Sg+∇⋅ρox2Kt⋅kroμo∇Po−ρoξ∇Z+ρgy3Kt⋅krgμg∇Pg−ρgξ∇Z=∑j=1Nρox2qoj+ρgy3qgj.The mass conservation equation of heavy oil components is as follows:(6)∂∂tϕρox3So+∇⋅ρox3K⋅kroμo∇po−ρoξ∇Z=∑j=1Nρox2qoj.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 theT−1 time step, the absolute permeability of the T time step is calculated. The specific calculation method is shown in the following equation:
(7)Kt=beta1×PV5+beta2×PV4+beta3×PV3+beta4×PV2+beta5×PV+beta6×T3+beta7×T2+beta8×T+beta9×Kt−1.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:(8)ΔQ=Q1+Q2+Q3−Q4.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:(9)∂∂tϕρwSwUw+ρoSoUo+ρgSgUg+1−ϕρRUR+∇⋅ρwUwVw+ρoUoVo+ρgUgVg=∇⋅λ∇T+∑j=1NρwUwqwj+ρoUoqoj+ρgUgqgj.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:(10)Sw+So+Sg=1.The molar fraction constraint equations are as follows:(11)x1+x2+x3=1,(12)y1+y2+y3=1.The phase equilibrium equations are as follows:(13)K1p,T=y2x1,(14)K2p,T=y3x2,(15)y1=PsatPg.The capillary force equations are as follows:(16)Pw=Po−PcowSw,(17)Pg=Po−PcogSg.When the production rate is constant, the internal boundary condition is expressed as follows:(18)qI=const.When the bottom hole pressure is constant, the internal boundary condition is expressed as follows:(19)pwf=const.Initial conditions are as follows:(20)P=Pit=0,T=Tit=0.


## 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 intoNx×Ny×Nz 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:(21)λi±1/2,j,kn+1=2λi,j,kn+1λi±1,j,kn+1λi,j,kn+1±λi±1,j,kn+1,λi,j±1/2,kn+1=2λi,j,kn+1λi,j±1,kn+1λi,j,kn+1±λi,j±1,kn+1,λi,j,k±1/2n+1=2λi,j,kn+1λi,j,k±1n+1λi,j,kn+1±λi,j,k±1n+1.The velocity of fluid in the grid is expressed by the following equations:(22)Vi±1/2,j,kn+1=±Kxμi±1/2,j,kn+1kr,i±1/2,j,kn+1Pi,j,kn+1−Pi±1,j,kn+1Δx,(23)Vi,j±1/2,kn+1=±Kxμi,j±1/2,kn+1kr,i,j±1/2,kn+1Pi,j,kn+1−Pi,j±1,kn+1Δy,(24)Vi,j,k±1/2n+1=±Kxμi,j,k±1/2n+1kr,i,j,k±1/2n+1Pi,j,kn+1−Pi,j,k±1n+1Δz.In the formula, the ratio of absolute permeability to viscosity is equal to the harmonic average of two adjacent grids.(25)Kxμi±1/2,j,kn+1=2Kx/μi±1,j,kn+1Kx/μi,j,kn+1Kx/μi±1,j,kn+1+Kx/μi,j,kn+1,Kyμi,j±1/2,kn+1=2Ky/μi,j±1,kn+1Ky/μi,j,kn+1Ky/μi,j±1,kn+1+Ky/μi,j,kn+1,Kzμi,j,k±1/2n+1=2Kz/μi,j,k±1n+1Kz/μi,j,kn+1Kz/μi,j,k±1n+1+Kz/μi,j,kn+1.Relative permeability takes the grid value in the direction of flow. The calculation method of relative permeability is shown in the following equations:(26)kr,i+1/2,j,kn+1=kr,i,j,kn+1Vi+1/2,j,kn+1≥0,kr,i+1,j,kn+1Vi+1/2,j,kn+1<0,kr,i−1/2,j,kn+1=kr,i−1,j,kn+1Vi−1/2,j,kn+1≥0,kr,i,j,kn+1Vi−1/2,j,kn+1<0,(27)kr,i,j+1/2,kn+1=kr,i,j,kn+1Vi,j+1/2,kn+1≥0,kr,i,j+1,kn+1Vi,j+1/2,kn+1<0,kr,i,j−1/2,kn+1=kr,i,j−1,kn+1Vi,j−1/2,kn+1≥0,kr,i,j,kn+1Vi,j−1/2,kn+1<0,(28)kr,i,j,k+1/2n+1=kr,i,j,kn+1Vi,j,k+1/2n+1≥0,kr,i,j,k+1n+1Vi,j,k+1/2n+1<0,kr,i,j,k−1/2n+1=kr,i,j,k−1n+1Vi,j,k−1/2n+1≥0,kr,i,j,kn+1Vi,j,k−1/2n+1<0.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:(29)ϕi,j,kn+1ρwSwi,j,kn+1−ϕi,j,knρwSwi,j,kn+ϕi,j,kn+1ρgy1Sgi,j,kn+1−ϕi,j,knρgy1Sgi,j,knΔxΔyΔz+ρwVwxi+1/2,j,kn+1−ρwVwxi−1/2,j,kn+1+ρgy1Vgxi+1/2,j,kn+1−ρgy1Vgxi−1/2,j,kn+1ΔyΔzΔt+ρwVwyi,j+1/2,kn+1−ρwVwyi,j−1/2,kn+1+ρgy1Vgyi,j+1/2,kn+1−ρgy1Vgyi,j−1/2,kn+1ΔxΔzΔt+ρwVwzi,j,k+1/2n+1−ρwVwzi,j,k−1/2n+1+ρgy1Vgzi,j,k+1/2n+1−ρgy1Vgzi,j,k−1/2n+1ΔxΔyΔt=∑j=1Nρwqwj+ρgy1qgjn+1Δt.The discrete form of mass conservation equation for dissolved gas components is as follows:(30)ϕi,j,kn+1ρox1Soi,j,kn+1−ϕi,j,knρox1Soi,j,kn+ϕi,j,kn+1ρgy2Sgi,j,kn+1−ϕi,j,knρgy2Sgi,j,knΔxΔyΔz+ρox1Voxi+1/2,j,kn+1−ρox1Voxi−1/2,j,kn+1+ρgy2Vgxi+1/2,j,kn+1−ρgy2Vgxi−1/2,j,kn+1ΔyΔzΔt+ρox1Voyi,j+1/2,kn+1−ρox1Voyi,j−1/2,kn+1+ρgy2Vgyi,j+1/2,kn+1−ρgy2Vgyi,j−1/2,kn+1ΔxΔzΔt+ρox1Vozi,j,k+1/2n+1−ρox1Vozi,j,k−1/2n+1+ρgy2Vgzi,j,k+1/2n+1−ρgy2Vgzi,j,k−1/2n+1ΔxΔyΔt=∑j=1Nρox1qoj+ρgy2qgjn+1Δt.The difference discrete form of mass conservation equation for light oil components is as follows:(31)ϕi,j,kn+1ρox2Soi,j,kn+1−ϕi,j,knρox2Soi,j,kn+ϕi,j,kn+1ρgy3Sgi,j,kn+1−ϕi,j,knρgy3Sgi,j,knΔxΔyΔz+ρox2Voxi+1/2,j,kn+1−ρox2Voxi−1/2,j,kn+1+ρgy3Vgxi+1/2,j,kn+1−ρgy3Vgxi−1/2,j,kn+1ΔyΔzΔt+ρox2Voyi,j+1/2,kn+1−ρox2Voyi,j−1/2,kn+1+ρgy3Vgyi,j+1/2,kn+1−ρgy3Vgyi,j−1/2,kn+1ΔxΔzΔt+ρox2Vozi,j,k+1/2n+1−ρox2Vozi,j,k−1/2n+1+ρgy3Vgzi,j,k+1/2n+1−ρgy3Vgzi,j,k−1/2n+1ΔxΔyΔt=∑j=1Nρox2qoj+ρgy3qgjn+1Δt.The discrete form of mass conservation equation for heavy oil components is as follows:(32)ϕi,j,kn+1ρox3Soi,j,kn+1−ϕi,j,knρox3Soi,j,knΔxΔyΔz+ρox3Vxoi+1/2,j,kn+1−ρox3Vxoi−1/2,j,kn+1ΔyΔzΔt+ρox3Vyoi,j+1/2,kn+1−ρox3Vyoi,j−1/2,kn+1ΔxΔzΔt+ρox3Vzoi,j,k+1/2n+1−ρox3Vzoi,j,k−1/2n+1ΔxΔyΔt=∑j=1Nρox3qojn+1Δt.The discrete form of the energy conservation equation is as follows:(33)ϕi,j,kn+1ρwSwUwi,j,kn+1+ρgSgUgi,j,kn+1+ρoSoUoi,j,kn+1−ϕi,j,kn+1ρwSwUwi,j,kn+ρgSgUgi,j,kn+ρoSoUoi,j,kn+1−ϕi,j,kn+1ρRURi,j,kn+1−1−ϕi,j,kn+1ρRURi,j,knΔxΔyΔz+ρwVxwUwi+1/2,j,kn+1−ρwVxwUwi−1/2,j,kn+1+ρoVxoUoi+1/2,j,kn+1−ρoVxoUoi−1/2,j,kn+1+ρgVxgUgi+1/2,j,kn+1−ρgVxgUgi−1/2,j,kn+1ΔyΔzΔt+ρwVywUwi,j+1/2,kn+1−ρwVywUwi,j−1/2,kn+1+ρoVyoUoi,j+1/2,kn+1−ρoVyoUoi,j−1/2,kn+1+ρgVygUgi,j+1/2,kn+1−ρgVygUgi,j−1/2,kn+1ΔxΔzΔt+ρwVzwUwi,j,k+1/2n+1−ρwVzwUwi,j,k−1/2n+1+ρoVzoUoi,j,k+1/2n+1−ρoVzoUoi,j,k−1/2n+1+ρgVzgUgi,j,k+1/2n+1−ρgVzgUgi,j,k−1/2n+1ΔxΔyΔt=λi+1/2,j,kn+1Ti+1,j,kn+1−Ti,j,kn+1Δx−λi−1/2,j,kn+1Ti,j,kn+1−Ti−1,j,kn+1ΔxΔyΔzΔt+λi,j+1/2,kn+1Ti,j+1,kn+1−Ti,j,kn+1Δy−λi,j−1/2,kn+1Ti,j,kn+1−Ti,j−1,kn+1ΔyΔxΔzΔt+λi,j,k+1/2n+1Ti,j,k+1n+1−Ti,j,kn+1Δz−λi,j,k−1/2n+1Ti,j,kn+1−Ti,j,k−1n+1ΔzΔxΔyΔt+∑j=1NρwqwjUwn+1+∑j=1NρoqojUon+1+∑j=1NρgqgjUgn+1Δt.


## 3.3. 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:(34)FX=0.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 iteration method, Taylor expansion of the m-th equation Fm can be obtained:
(35)FmX+δX=FmX+∑l=14×Nx×Ny×Nz∂Fm∂xlδxl+oδX2.Equation (35) can be expressed as follows:
(36)FmX+δX=FmX+J⋅δX+oδX2.In Equation (36), J is the Jacobian determinant and can be calculated by the following equation:
(37)Jml=∂Fm∂xl=Fmx1,x2,⋯,xl+Δxl,⋯,x4×Nx×Ny×Nz−Fmx1,x2,⋯,xl,⋯,x4×Nx×Ny×NzΔxl.LetFmX+δX=0 and ignore oδX2, then you can get the following:
(38)J⋅δX=−F.DeterminantJ 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.
(39)Xnew=Xold+δX.Continuous iteration until the following convergence conditions are satisfied:(40)14×Nx×Ny×Nz∑l=14×Nx×Ny×NzδxLxl<ε.Through the above steps, the solutions of the discrete equations (Equations (29), (30), (31), (32), and (33)) can be obtained.


## 4. 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 high-temperature and high-pressure experimental device is shown in Figure6. The experimental model is high-pressure sand-packed model (with temperature and pressure sensor, temperature control heating jacket, and thermal insulation layer), intermediate vessel (fluid buffer vessel for huff and puff), steam generator, CO2 intermediate vessel, N2 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.Figure 6
The flowchart of the multicomponent thermal fluid high-temperature and high-pressure experimental device.Experimental steps are shown as follows:(1)
Clean the actual oil sand(2)
Oil sand is filled into the model. The model volume is 5.1 L, the pore volume is 1.179 L, the irreducible water volume is 0.295 L, and the saturated oil volume is 0.884 L(3)
Saturate the model with formation water at the formation temperature (56°C) and then place the model for 12 hours(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 mLsteam+N2/CO2 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 pressureAccording 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 Table3. 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.Table 3
Basic parameters of numerical simulation model.
ParameterValueParameterValueMesh dimension (x×y×z)40×1×1X-direction mesh step (cm)1Y-direction mesh step (cm)2.21Z-direction mesh step (cm)2.21Porosity (%)33Initial oil saturation (%)75Permeability (10-3 μm3)1200Initial water saturation (%)2556°C oil viscosity (mPa·s)1681Water viscosity (mPa·s)1Pore volume (L)1.179Initial formation pressure (MPa)20Figure 7
Mesh of numerical simulation model.Figure 8
Water-oil relative permeability curves.Figure 9
Gas-oil relative permeability curves.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 Figure10.Figure 10
Physical simulation results and numerical simulation results of multicomponent thermal fluid huff and puff (scatter points represent experimental results, red curves represent numerical simulation results considering formation damage, and black curves represent numerical simulation results without considering formation damage).As can be seen from Figure10, 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.


## 5. 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(3)
Based on the formation damage experimental results, a new multicomponent thermal fluid formation damage mathematical model is established by using multielement nonlinear regression method. The calculated results of the new multicomponent thermal fluid reservoir damage mathematical model are basically consistent with the formation damage experimental results(4)
The new multicomponent thermal fluid formation damage model is embedded in the traditional compositional numerical simulator, and a new multicomponent thermal fluid reservoir damage numerical simulation model is established. Comparing with multicomponent thermal fluid huff and puff experimental results and traditional numerical simulator calculation results, it is found that the new multicomponent thermal fluid formation damage numerical simulator calculation results are closer to the experimental results



---

*Source: 8845896-2021-01-15.xml*