2. 中国石油天然气集团公司油气藏改造重点实验室二氧化碳压裂增产研究室;
3. 上海交通大学工程热物理研究所
2. CO2 Fracturing Stimulation Laboratory of Key Laboratory of Reservoir Stimulation, CNPC;
3. Institute of Engineering Thermophysics, Shanghai Jiaotong University
0 引言
CO2来源广泛，对储层伤害小[1-5]，是一种经济又环保的备选压裂液；同时超临界状态下的CO2具有与气体相近的扩散能力，密度与液体相近，能高效溶于原油而降低稠油的黏度[1-5]，从而易于驱替出水基压裂液难以作用的油气资源，提高采收率。CO2压裂技术正受到越来越多的关注。为评估CO2的压裂效果和人造裂缝质量，获得较准确的井底温度和压力很有必要。通常情况下CO2压裂井井深在2 000 m以上，井底的压力较高，直接测量对传感器元件和施工操作要求高，需要投入的成本高。借助地面方便获得的入口施工参数和前期勘测得到的地质参数，通过理论分析来获得井筒内温度和压力分布和具体的井底压力具有明显的优势。A.R.HASAN等[4, 6]基于H.S.RAMEY[7]在1962年提出的半稳态模型，考虑油井的倾斜角度等影响因素，提出了适用性更广的稳态计算模型，该模型可用于计算长时间注液后流体温度和压力沿井筒深度方向的分布，进而推算得到井底的压力。然而该模型未涉及时间项，并不适用于短时间注液过程，如压裂或酸化等非稳态过程。王鸿勋等[8]提出的非稳态流动模型可用于计算从压裂开始到趋于稳定过程中的井筒流体与周围环境热交换。不同于J.R.EICKMEIER等[9]的数值方法，王鸿勋等[8]考虑了井筒中原有积液对注入液换热的影响，并采用全隐式的离散方法达到数值解的无条件稳定。但是该方法主要针对的是水力压裂过程，计算中将压裂液物性参数处理成常数，并忽略了流动摩擦对流体温度和压力的影响。
不同于常规的水基压裂液，CO2热物性在压裂工况下变化大，不能简化为常数考虑，同时CO2压裂的入口排量大，导致井筒中的流动摩阻不可忽略。针对CO2的这2大主要特点，笔者提出了井筒内耦合流体热物性的非稳态流动传热模型。为获得较准确的变物性规律，笔者分别采用Vesovic模型[10-12]和Span-Wagner方法[13]计算黏度、比热容和密度等参数。耦合了变物性的非稳态流动传热模型难以得到精确的分析解，常用的方法是对井筒空间和施工时间各自取离散微元后数值求解控制方程。笔者综合考虑液体变物性和黏性生热对温度和压力的影响，通过数值分析得到了液体温度压力随时间和井深的变化关系，并重点探讨了井底的温度压力随时间的变化情况，所得结论可为地层裂缝的计算提供更准确的入口参数。
1 数学模型的建立井筒结构如图 1所示。图中，dti和dto分别表示油管内径和外径，dci、dco分别表示套管内径和外径，dcem为水泥环壁的厚度。主要研究对象为垂向上的一维流动，假设同一井深处的流体温度和压力相等。为方便求解，对模型中的流动传热还做了以下几个合理假设[8]。
|
| 图 1 井筒示意图 Fig.1 Schematic diagram of the wellbore |
(1) 忽略井筒流体和地层岩石在垂向上的热传导；
(2) 地面排量及注液温度和压力为定值，不随时间变化；
(3) 注液前油管和环空中的流体稳定流动；
(4) 油/套管的管径不随井深发生变化；
(5) 油管周围地层各向同性，为匀质地层；
(6) 假设地表到恒温点温度恒定为Teio，恒温点z1以下的岩石温度随井深z呈线性递增，其表达式为：
|
(1) |
式中：tg为地温梯度，Tei为各点的地层温度。
在上面假设的基础上，对dz流体微元应用动量定理可以得到，开口控制体的动量变化等于净流入的动量与合外力之和，表达式为：
|
(2) |
式中：p为各点的流体压力，ρ为流体密度，g为当地重力加速度，f为压裂液在井筒流动时的摩阻系数，采用N.H.CHEN[14]在1979年提出的关联式计算，v为流体速度，t为压裂时间。
同样对流体微元应用能量守恒定理，流体微元的能量变化率等于净流入的能量与黏性生热和外层壁面的导热量之和：
|
(3) |
式中：T为流体温度，cp为流体的定压比热容，μj为流体的节流系数，Q为油管单位长度上流体与地层的换热量，Q1为单位体积流体的黏性力产热功率。
一般情况下，可以认为产热功率Q1等于黏性力做功。于是流体微元存在下面的机械能守恒式：
|
(4) |
|
(5) |
式中：ek为单位质量流体的机械能，
控制体中的流体满足下面的连续性方程：
|
(6) |
联立式(3)~式(5)可以得到稳态和非稳态情况下流体的黏性力做功功率。
稳态时油管内单位体积的产热功率为：
|
(7) |
非稳态时油管内单位体积的产热功率为：
|
(8) |
基于总热阻的假设，单位长度油管内流体与地层的换热量Q可以写为：
|
(9) |
式中：L为基于油管内壁换热面积的总热阻。
油管壁、套管壁和水泥环之间由热传导传递热量，油管和套管内的流体与管壁之间对流换热，总热阻的计算式为：
|
(10) |
式中：hti、hto和hci分别为油管内壁、油管外壁和套管内壁的流体对流换热系数，kt、kc和kcem分别为油管壁、套管壁和水泥环壁的导热系数。
根据式(4)~式(10)，可将式(3)写成：
|
(11) |
结合给定的初始和边界条件，求解式(2)~式(10)可以得到井筒中流动换热的温度与压力随时间和井深的变化规律。对于液态CO2的热物性参数，采用P-R两参数方程[10]计算流体的密度，使用文献[13]提出的基于亥姆霍兹自由能的方法计算CO2流体跨临界流动中的定压比热容，利用Vesovic[11-12]方程计算CO2的黏度变化，实验关联式中用到的参数如各自文献中所示。
P-R两参数方程为：
|
(12) |
式中：R是通用气体常数，a、b是2个与气体种类有关的参数，V为气体摩尔体积。
Span-Wagner关联式为：
|
(13) |
|
(14) |
式中：A为亥姆霍兹自由能；Φ为无因次亥姆霍兹自由能；Φ0为理想部分亥姆霍兹自由能；Φr为残余部分亥姆霍兹自由能；δ为对比密度；δ=ρ/ρc；τ为逆对比温度，τ=Tc/T；Cp*和ΔCp分别为理想和残余比热容。
Vesovic模型为：
|
(15) |
式中：η是流体黏度，η0为理想部分黏度，Δη为残余部分黏度，Δηc为超临界误差黏度。
2.1 控制方程的离散将压力控制方程式(2)按空间向后差分，时间中间差分离散为：
|
(16) |
其中：
|
(17) |
式中：下标i和上标j分别为离散网格中的空间节点和时间节点，dz为井筒轴向上的空间步长，dt为时间步长；下同。
为保证解的稳定性，式(11)采用全隐式离散为：
|
(18) |
其中：
|
(19) |
将式(19)化为：
|
(20) |
|
(21) |
|
(22) |
对每一个流体微元分别求出参数A和B的值，带入式(20)后可以求得
基于式(17)~式(22)的离散方法，结合实际工况的参数可以计算得到井筒的温度和压力分布。为验证模型的可靠性，笔者以文献[15-17]中给出的CO2压裂参数作为输入条件，将计算得到的温度和压力分布与实验测量值进行对比，结果如图 2和图 3所示。
|
| 图 2 第1次注CO2 13 h后温度和压力沿井深的分布曲线 Fig.2 Distribution of temperature and pressure along the well after 13 h of the first injection of CO2 |
|
| 图 3 第2次注CO2 26.58 d后温度和压力沿井深的变化曲线 Fig.3 Distribution of temperature and pressure along the well after 26.58 d of the first injection of CO2 |
从图 2可以看到，笔者的模型计算得到的井底温度和压力分别为377.81 K(104.66 ℃)和53.184 MPa，与文献中给出的实测温度值374.15 K(101 ℃)和实测压力值52.02 MPa比较，相对误差分别为0.98%和2.24%。从图 2还可以看到，模拟结果与现场测量值接近，最大温度偏差为0.96%，出现在井底处。
目前缺乏超临界CO2压裂的实测数据，笔者提出的模型可在有限的已知参数下，计算得到CO2在井筒中的温度和压力场变化规律。下面以榆林气田的神62井油管注液为例，分析各参数对井筒内的温度和压力的影响。因为不同的初始条件和边界条件对结果的影响较大，求解过程中用到的参数如表 1所示。
| 参数 | 数值 |
| 排量/(m3·min-1) | 1.5/2.0 |
| 表层套管深度/m | 514.19 |
| 恒温点深度/m | 30 |
| 井筒流体入口温度/℃ | 0 |
| 井筒流体入口压力/MPa | 25/50 |
| 油/套管壁导热系数/(W·m-1·℃-1) | 40.66 |
| 油/套管壁比热容/(J·kg-1·℃-1) | 460.5 |
| 每100 m地温梯度/℃ | 3.0 |
| 钻头所在深度/m | 3 000 |
| 恒温点温度/℃ | 20 |
| 油管内径/mm | 76.00 |
| 油管外径/mm | 88.90 |
| 套管内径/mm | 121.36 |
| 套管外径/mm | 139.70 |
| 水泥壁导热系数/(W·m-1·℃-1) | 0.6 |
由式(11)可以看到，温度主要受井口流体温度、流动时黏性生热率和地层加热量的影响。图 4为是否考虑黏性产热量时井底温度对比曲线，由图可以清楚地看到黏性力对温度影响比较大，图示中井底流体的温度差值已经达到了10 ℃；因此在对井筒进行换热分析时，不能忽略摩擦生热的影响。对于排量较大的CO2压裂过程，流速快导致摩阻梯度偏大，换热分析时必须考虑黏性力对温度场的影响。
|
| 图 4 黏性生热对井底温度的影响 Fig.4 Effect of viscous heat generation on bottomhole temperature |
非稳态流动换热时的初始条件对最终结果的影响较大，王鸿勋等[9]的研究中选取了井筒积液与地层已经达到热平衡的温度分布作为非稳态计算的初始条件，笔者则选取了低排量、低压注液时，流体已经在井筒中达到稳态流动时的温度和压力分布作为非稳态计算的初始条件。初始条件的温度压力场如图 5所示。由图可以看到，温度随井深增加而呈非线性升高，压力随井深增加而逐渐增大。
|
| 图 5 稳态时的井筒内温度和压力沿程变化曲线 Fig.5 The temperature and pressure in the wellbore at steady state |
非稳态流动换热时井筒内的温度和压力分别如图 6和图 7所示。由图 6和图 7可知，流体温度随井深增加而升高，随低温流体持续注入而下降；在重力作用下压力随井深而上升，随压裂进行而增加。不同时刻流体压力随井深变化曲线如图 8所示。由图可以观察到不同时刻井筒内的压力变化趋势相差不大，随时间增加压力上升总量略有增加。流体沿井筒流动，与地层之间的温度差随着井深增加而增加，井筒内流体受到地层的加热使得温度随井深升高，变化趋势如图 9所示。
|
| 图 6 井筒内温度随时间和井深变化曲面图 Fig.6 The temperature variation in the wellbore with time and depth |
|
| 图 7 井筒内压力随时间和井深变化图 Fig.7 The pressure variation in the wellbore with time and depth |
|
| 图 8 不同时刻流体压力随井深变化曲线 Fig.8 Fluid pressure distribution in wellbore at different times |
|
| 图 9 不同时刻流体温度随井深变化曲线 Fig.9 Fluid temperature distribution in wellbore at different times |
井筒内流体的温度同时受井口低温CO2的影响，施工初期的温度沿井深变化剧烈，1 min时CO2温度沿3 000 m井深上升了45 K左右，而到10 min时上升量仅为10 K左右。这说明随着井口流体的泵入，井筒内的流体温度逐渐趋于稳定，并且越来越接近流体的入口温度。井口不断泵入的CO2的冷却作用和井筒原始积液逐渐流出井筒，导致了井底温度随注液时间增加而下降，最后趋于稳定，这与王鸿勋等[8]的试验测量结果的变化趋势吻合。图 10、图 11和图 12分别比较了井筒不同深度处的温度和压力随时间的变化情况。
|
| 图 10 井筒内不同位置处的温度随时间变化曲线 Fig.10 Changes of temperature at different locations versus time |
|
| 图 11 井筒内不同位置处的压力随时间变化曲线 Fig.11 Changes of pressure at different locations versus time |
|
| 图 12 施工前期井筒内不同位置处的压力随时间变化曲线 Fig.12 Pressure changes at different locations during the pre-construction period |
从对比的情况可以看到，离井口越近，温度和压力越易达到稳定。这是因为不管是压力还是温度，随着井深递增受原始积液、流体扰动和地层岩石的影响都在加剧，越难在流动中达到动态平衡。
3 结束语笔者在研究中考虑了摩阻压降和流体变物性对流体温度的影响，在控制方程中引入了与流体扰动对应的热功率项。数值模拟结果表明：往井筒内大排量注入CO2时，黏性生热量占主导地位，对井筒中温度、压力分布和流体热物性有显著的影响；CO2温度沿井深呈非线性升高；当重力的作用大于井筒壁面摩擦力的作用时，井筒内压力沿井深增加；油井越深，井底的温度和压力达到稳定状态需要的时间越长。
虽然笔者计算的井筒压裂液为低温CO2流体，但推导过程和所建立的数学物理模型同样适用于其他类型压裂液的注入过程；计算时考虑不同液体的物性参数和流动中的状态变化，可以得到不同压裂液在井筒中的温度压力分布；适当修正控制方程也可以使模型适用于生产井和水平井等其他油井的施工过程。
| [1] |
孙宝江, 孙小辉, 王志远, 等. 超临界CO2钻井井筒内流动参数变化规律[J].
中国石油大学学报(自然科学版), 2016, 40(3): 88-95.
SUN B J, SUN X H, WANG Z Y, et al. Flow behavior analysis during supercritical CO2 drilling:Consideration of varying thermo-physical properties of CO2 in wellbore[J]. Journal of China University of Petroleum, 2016, 40(3): 88-95. DOI: 10.3969/j.issn.1673-5005.2016.03.011 |
| [2] |
王在明.超临界CO2钻井液特性研究[D].青岛: 中国石油大学(华东), 2008. WANG Z M. Feature research of supercritical carbon dioxide drilling fluid[D]. Qingdao: China University of Petroleum(Huadong), 2008. |
| [3] | GUPTA A P, GUPTA A, LANGLINAIS J. Feasibility of supercritical carbon dioxide as a drilling fluid for deep underbalanced drilling operation[R]. SPE 96992, 2005. |
| [4] | HASAN A R, KABIR C S. Heat transfer during two-phase flow in wellbores: Part Ⅰ-formation temperature[R]. SPE 22866, 1991. https://www.onepetro.org/conference-paper/SPE-22866-MS |
| [5] |
王瑞和, 倪红坚. CO2连续管井筒流动传热规律研究[J].
中国石油大学学报(自然科学版), 2013, 37(5): 65-70.
WANG R H, NI H J. Wellbore heat transfer law of carbon dioxide coiled tubing drilling[J]. Journal of China University of Petroleum, 2013, 37(5): 65-70. DOI: 10.3969/j.issn.1673-5005.2013.05.009 |
| [6] | HASAN A R, KABIR C S, WANG X W. A robust steady-state model for flowing-fluid temperature in complex wells[J]. SPE Production & Operations, 2009, 24(2): 269-276. |
| [7] | RAMEY H S. Wellbore heat transmission[J]. Journal of Petroleum Technology, 1962, 14(4): 427-435. DOI: 10.2118/96-PA |
| [8] |
王鸿勋, 李平. 水力压裂过程中井筒温度的数值计算方法[J].
石油学报, 1987, 8(2): 91-99.
WANG H X, LI P. Numerical method for calculating wellbore temperature during hydraulic fracturing[J]. Acta Petrolei Sinica, 1987, 8(2): 91-99. |
| [9] | EICKMEIER J R, ERSOY D, RAMEY H J. Wellbore temperatures and heat losses during production or injection operations[J]. Journal of Canadian Petroleum Technology, 1970, 9(2): 115-121. |
| [10] | GASEM K A M, GAO W, PAN Z, et al. A modified temperature dependence for the Peng-Robinson equation of state[J]. Fluid Phase Equilibria, 2001, 181(1/2): 113-125. |
| [11] | CHUNG T H, AJLAN M, LEE L L, et al. Generalized multiparameter correlation for nonpolar and polar fluid transport properties[J]. Industrial & Engineering Chemistry Research, 1988, 27(4): 671-679. |
| [12] | FENGHOUR A, WAKEHAM W A, VESOVIC V. The viscosity of carbon dioxide[J]. Journal of Physical and Chemical Reference Data, 1998, 27(1): 31-44. DOI: 10.1063/1.556013 |
| [13] | SPAN R, WAGNERW. A new equation of state for carbon dioxide covering the fluid region from the triple-point temperature to 1100 K at pressures up to 800 MPa[J]. Journal of Physical and Chemical Reference Data, 1996, 25(6): 1509-1596. DOI: 10.1063/1.555991 |
| [14] | CHEN N H. An explicit equation for friction factor in pipe[J]. Industrial & Engineering Chemistry Fundamentals, 1979, 18(3): 296-297. |
| [15] | LYU X R, ZHANG S C, MA X F, et al. Numerical investigation of wellbore temperature and pressure fields in CO2fracturing[J]. Applied Thermal Engineering, 2018, 132: 760-768. DOI: 10.1016/j.applthermaleng.2017.12.095 |
| [16] |
窦亮彬, 李根生, 沈忠厚, 等. 注CO2井筒温度压力预测模型及影响因素研究[J].
石油钻探技术, 2013, 41(1): 76-81.
DOU L B, LI G S, SHEN Z H, et al. Wellbore pressure and temperature prediction model and its affecting factors for CO2 injection wells[J]. Petroleum Drilling Techniques, 2013, 41(1): 76-81. DOI: 10.3969/j.issn.1001-0890.2013.01.015 |
| [17] |
毕刚, 李根生, 沈忠厚, 等. 超临界CO2径向水平井钻井技术应用前景分析[J].
石油机械, 2013, 41(6): 14-19.
BI G, LI G S, SHEN Z H, et al. Analysis of the application prospect of supercritical CO2 radial horizontal drilling technology[J]. China Petroleum Machinery, 2013, 41(6): 14-19. DOI: 10.3969/j.issn.1001-4578.2013.06.004 |