赵帅, 孙友宏, 杨秦川, 李强
吉林大学 建设工程学院, 吉林 长春 130026
收稿日期：2018-05-11
基金项目：国家自然科学基金青年基金资助项目(51404107)；吉林大学博士研究生交叉学科科研资助计划(10183201839)。
作者简介：赵帅(1990-)，男，山东临沂人，吉林大学博士研究生;
孙友宏(1965-)，男，江苏如皋人，吉林大学教授，博士生导师。

摘要：采用热-流耦合分析模式, 对水力压裂之后扶余油页岩储层的传热导流渗透能力进行数值模拟, 发现流体主要沿油页岩层理方向形成地裂隙流出, 但是随温度的增加, 孔隙度增大也会有少许流体从油页岩的原生孔隙流出, 渗流场压力在同一截面自裂隙垂直于油页岩层理向两端呈现下降趋势; 流体对油页岩地层热量的传导主要是沿裂隙方向进行.加热时间增长, 裂隙两侧油页岩裂解, 孔隙度增加, 氮气向油页岩储层的扩散速度也得到了提高, 加热至40 d之后, 裂隙周围油页岩首先达到裂解温度, 加热至60~100 d, 油页岩层的平均温度自500 K提升至650 K, 整个油页岩能够被有效热解.
关键词：油页岩水力压裂渗透率热传导数值模拟
Thermal-Fluid Coupling Simulation on In-situ Unsteady Pyrolysis of Fuyu Oil Shale
ZHAO Shuai, SUN You-hong, YANG Qin-chuan, LI Qiang
Construction Engineering College, Jilin University, Changchun 130026, China
Corresponding author: LI Qiang, E-mail: liqiangjlu @163.com
Abstract: Numerical simulation on heat transfer and permeability of Fuyu oil shale reservoir after hydraulic fracturing was carried out by using thermal-fluid coupling analysis model. It is found that the fluid mainly flows out of the fracture along the oil shale bedding direction. However, as the temperature increases, the porosity increases and a little fluid flows out of the primary pores of the oil shale. The seepage field pressure shows a downward trend in the same section from the fracture perpendicular to the oil shale bedding. The heat conduction of fluid to oil shale formation was mainly in the direction of fracture. As the heating time increases, the oil shale cracks on both sides of the fissure, and with the porosity increase, the diffusion rate of nitrogen to the oil shale reservoir is also improved. After 40 days of heating, the oil shale around the fracture first reaches the cracking temperature. Between 60 and 100 days of heating, the average temperature of oil shale rises from 500 K to 650 K, and the whole oil shale can be effectively cracked.
Key words: oil shalehydraulic fracturingpermeabilityheat conductionnumerical simulation
原始的油页岩储层结构致密, 孔隙度低, 渗透率低, 加热之后干酪根虽然能够被裂解, 但由于储层导流能力差, 干酪根裂解生成的部分产物不能被有效采出, 导致原本含油率低的油页岩在常规原位裂解工艺条件下采收率低下^[1-4].目前油页岩油的提取主要依靠地面低温干馏技术, 该方法能够获得较高的油收率, 但是工艺较复杂, 占地面积大, 容易导致环境污染, 生产成本高.油页岩原位热解工艺, 是通过储层改造的方式在加热井和开采井之间的油页岩层制造裂缝, 并随压裂液泵入陶粒支撑剂, 支撑裂缝, 提高油页岩层导流能力，增加传热介质与油页岩的接触面积，提高油页岩层导热速率^[5-6].
国内外对非常规油气藏水力压裂提高采收效率的模拟较多.Ma等对煤层气储层进行水利压裂模拟表明, 储层改造可以提高生产采收效率^[7].Zhang等建立了煤层气两相压裂三维模型, 分析了渗透率、原始体积密度、孔隙度等对生产的影响^[8].Tang等采用固定网格模型模拟非常规油藏多级水利压裂的裂缝扩展, 预测油井产能^[9].对于像油页岩一样的多孔介质在热-流耦合作用下的渗流传热特征研究较少, 本文以扶余油页岩的原位裂解先导试验为背景, 从传热学角度研究油页岩储层水力压裂之后, 提高油页岩层的孔隙率、渗透性能与导流能力的行为, 注入高温高压氮气加热油页岩层及氮气在油页岩层的渗流动力学特征.以多物理场耦合有限元分析软件进行数值模拟, 对氮气在多层裂隙内的流动与传热特性进行分析, 为油页岩原位注热开采提供指导.
1 模型与参数1.1 油页岩物性特征为探索扶余油页岩有效裂解温度, 进行了氮气气氛、不同升温速率下热失重区间的测试, 如图 1所示, 温度处于350~550 ℃为主要的失重过程, 即有机质(干酪根)的主要分解反应阶段, 也为主产油阶段.
图 1(Fig. 1)

图 1 扶余油页岩热重分析的TG和DTG曲线Fig.1 TG and DTG curve of Fuyu oil shale

如图 2所示, 该温度区间内, 油页岩内赋存的干酪根发生裂解, 产生的油气会沿着油页岩内部的孔隙排出, 但是部分干酪根富集在油页岩内部一个密闭的空间内, 受热裂解产生的油气不能及时排出, 因此密闭空间的压力迅速增加, 直到温度继续升高, 达到油页岩的应力极限, 密闭空间发生破坏, 产生裂隙甚至裂缝, 因此随温度的升高油页岩的孔隙度是不断变大的, 如表 1所示.
图 2(Fig. 2)

图 2 油页岩裂解过程中孔隙的变化Fig.2 Changes in porosity during the pyrolysis of oil shale (a)—350 ℃；(b)—550 ℃.

表 1(Table 1)

表 1 扶余油页岩热物理性能Table 1 Thermophysical properties of Fuyu oil shale 温度 导热系数/(W·m-1·K-1) 比热容 孔隙度 K 平行层理 垂直层理 J·kg-1·K-1 298 0.615 5 0.509 5 2 206.0 0.024 5 423 0.695 1 0.550 6 2 004.0 0.047 6 523 0.406 8 0.328 2 1 733.0 0.062 3 623 0.363 7 0.336 9 970.7 0.187 9 723 0.356 8 0.332 0 1 187.0 0.405 3 773 0.232 4 0.217 7 1 033.0 0.558 2

表 1 扶余油页岩热物理性能 Table 1 Thermophysical properties of Fuyu oil shale

1.2 模型参数扶余油页岩原位裂解试验工程位于吉林省扶余市永平乡, 油页岩层上覆地层主要为泥岩和页岩, 油页岩的埋深在477~486 m, 储层厚度为9 m, 在水力压裂过程中, 实时进行微震数据的监测, 微震监测数据如图 3所示.经过水力压裂改造形成2条垂直于油页岩层理方向的裂隙把油页岩储层分为3层.注热井作为高温高压氮气的输入井, 氮气沿射孔段流入油页岩地层, 并且通过地层的裂隙流至开采井.
图 3(Fig. 3)

图 3 扶余油页岩水力压裂微震监测与模型建立图Fig.3 Microseismic monitoring and modeling map of hydraulic fracturing in Fuyu oil shale

模型假设:1)裂隙是沿油页岩层理方向扩展的; 2)由于氮气主要沿裂隙进行传热, 因此假设沿x，y方向油页岩的导热系数是相同的, 并且大于沿z方向的导热系数; 3)氮气沿裂隙各方向的流速是均匀的.氮气的计算参数如表 2所示.
表 2(Table 2)

表 2 氮气的计算参数Table 2 Calculation parameters of nitrogen 参数 单位 数值 定压比热容 kJ/(kg·K) 1.038 定容比热容 kJ/(kg·K) 0.741 密度ρ g/cm3 1.16 气体黏度μ Pa·s 175.44×10－7 导热系数λ W/(m·K) 0.024 75 输入温度T K 773 气体常数R J/(mol·K) 8.314 4 流量Q m3/min 11 输入压力P MPa 8.5 比热容比γ — 1.401 平均热膨胀系数β 1/K 0.007 53 压缩系数z — 0.292

表 2 氮气的计算参数 Table 2 Calculation parameters of nitrogen

2 模拟分析2.1 渗流场压力损失根据油页岩原位裂解工艺需求, 高温高压氮气自射孔端流入沿油页岩层水平层理方向的裂隙, 具体输入参数如表 1所示.沿油页岩层X，Y方向, 气体在单个裂隙中压力损失模型^[10]满足:

(1)

(2)

(3)

f(ε)是与形状有关的因数, 由于开采井与加热井之间水力压裂产生的裂缝有效宽度仅为5 mm, 裂缝的长度15 m, 长宽比大于3 000, 根据经验f(ε)取值0.994^[10], 无滑移状态下b=0.
在油页岩的微观孔隙裂隙中, 根据达西-韦斯拉科方程^[11]，多孔介质中压力的损失为

(4)

氮气在多孔介质内的流速为

(5)

在热解过程中, 氮气向油页岩原生孔隙的运移满足多孔介质的哈根-波肃叶方程^[11]：

(6)

压力的衰减主要沿裂隙进行, 由于裂隙并非直线型且表面非光滑, 沿油页岩层理方向裂隙的高度仅为5 mm, 因此压力的沿程损失与局部损失均比较大, 随注热井注入氮气压力的升高, 开采井压力也会逐渐升高.在氮气输入压力为8.5 MPa情况下，自裂隙出口输出的氮气压力仅为2.78 MPa, 如图 4所示.
图 4(Fig. 4)

图 4 压力分布云图(单位：MPa)Fig.4 Pressure distribution nephogram (a)—射孔端输入压力；(b)—裂隙压力分布.

渗流场的压力分布如图 5所示, 高压驱替流体渗流主要发生在沿油页岩层理方向形成地裂隙中, 渗流场压力在同一截面自裂隙垂直于油页岩层理向两端呈现下降趋势, 沿平行于油页岩层理方向, 压裂形成地裂隙非直线并且壁面粗糙导致沿程阻力较大, 因此压力损失也会较快, 随着压力增加, 微裂缝提高致密介质渗透率的作用增强.
图 5(Fig. 5)

图 5 裂隙及两井的压力分布拟合曲线Fig.5 Fitting curves of pressure distribution in fractures and two wells (a)—裂隙压力；(b)—井压力.

高压驱替作用下, 流体主要从沿油页岩层理方向形成地裂隙流出, 但是也会有流体从油页岩的原生孔隙流出.根据流线分布图 6可以看出, 高压气体在沿裂隙扩散的同时, 也会向油页岩原生储层扩散, 这有利于原位注热过程中, 提高热氮气对地层的热传导能力.
图 6(Fig. 6)

图 6 渗流场流线分布图(单位：MPa)Fig.6 Streamline distribution of seepage field (a)—地层流线分布；(b)—裂隙流线分布.

2.2 温度场分布根据Brinkman-Forchheimer扩展的Darcy渗流模式^[12], 氮气在通过裂隙对油页岩层的热传导过程中质量、动量与能量守恒.

(7)

(8)

(9)

(10)

F为多孔介质基体对氮气的作用力:

(11)

(12)

(13)

(14)

在上述模型中加入高温氮气对油页岩的热传导的温度场, 可以获得热-流耦合作用下, 氮气对油页岩的传热效果.按照表 3所示油页岩层的计算参数, 分别设置加热时间为10，20，40，60，100，150 d, 沿裂隙每米设置温度监测点, 获得氮气对油页岩层的加热温度云图和裂隙轨迹温度分布随时间的变化情况.
表 3(Table 3)

表 3 油页岩层的计算参数表Table 3 Calculation parameters of oil shale 参数 单位 数值 导热系数 W/(m·K) 0.65 岩层温度T0 K 288 颗粒平均直径dp μm 110 裂隙长度L0 m 15 裂隙宽度L m 9 支撑剂尺寸w mm 0.5 密度ρo g/cm3 1.80 地层压力P0 MPa 9.2 重力加速度g m/s2 9.8

表 3 油页岩层的计算参数表 Table 3 Calculation parameters of oil shale

随加热时间的增长, 油页岩地层裂隙内的温度最先升高, 但由于裂隙长达15 m, 沿裂隙热损失较大.平行于层理的方向氮气对油页岩地层的传热与裂隙的形状相关, 裂隙与油页岩层中轴线夹角呈锐角的一侧, 温度会优先升高, 如图 7所示.这是因为高温高压氮气沿裂隙输入, 由于输入角度的不同, 会优先冲击裂隙的一侧, 此外高温高压会在裂隙拐角的一侧形成较强的湍流, 强化换热过程.氮气加热油页岩层的速度在沿油页岩层理方向要高于垂直于油页岩层理方向, 这一方面是因为, 水力压裂产生的裂隙主要沿油页岩层理方向, 氮气的流动通道主要是沿裂隙进行, 因此主要的换热面是沿油页岩层理方向裂隙的表面; 另一方面是因为, 在沿层理方向油页岩的导热系数大于垂直于层理方向油页岩的导热系数.
图 7(Fig. 7)

图 7 平行于层理方向油页岩温度变化Fig.7 Temperature change of oil shale parallel to the bedding direction (a)—10 d；(b)—20 d；(c)—40 d；(d)—60 d；(e)—100 d；(f)—150 d.

2.3 升温速率的变化油页岩地层孔隙度的变化会影响传热传质的进程, 油页岩地层孔隙度的增加, 会导致载热流体扩散距离的增加, 对远离裂隙处的油页岩地层进行热量的传导.对不同温度下油页岩孔隙度变化的测试结果如图 8所示.
图 8(Fig. 8)

图 8 温度与孔隙度、出口流量的关系Fig.8 Temperature versus outlet flow and porosity

由于裂隙两侧油页岩层温度均已提升, 孔隙度提升, 热损失减少, 导流能力提高, 有一部分高温高压氮气沿裂隙进入油页岩地层, 并进行热传导的过程.此外, 由于驱替流压的存在, 高温高压氮气会被吸附在有机质热解产生的连通性较差的孔隙内; 而对于孔隙连通性较好的油页岩地层, 高温高压氮气对油页岩地层除了进行热量的传导之外, 还会渗流进入油页岩层驱替油页岩热解产生的油气产物, 有利于提高油气采收效率.
由于加热前期裂隙周围温度较低, 高温高压氮气与裂隙的温差较大, 因此随加热时间的延长裂隙温度上升迅速.如图 9所示, 当加热至100 d以后, 裂隙温度达到700 K, 氮气温度773 K, 温差减小, 因此传热过程变得缓慢; 氮气对油页岩地层的传热在80 d之前传热缓慢, 这是因为储层裂隙温度较低, 不能达到油页岩裂解温度, 因此其孔隙度相对较低, 氮气渗透扩散到油页岩储层的进程缓慢, 但是80 d之后, 裂隙温度升至650 K达到油页岩裂解的温度, 裂隙周围的油页岩首先裂解, 产生孔隙, 随时间的延长孔隙增加, 因此氮气渗透扩散的速度增加, 进而油页岩储层温度迅速增加.
图 9(Fig. 9)

图 9 裂隙与地层温度变化及升温速率Fig.9 Temperature change and temperature rising rate in fractures and stratum (a)—平均温度变化；(b)—升温速率.

3 结论1) 加热以后, 氮气首先沿水力压裂产生的裂隙对油页岩层进行热传导, 当油页岩层温度升高之后, 孔隙度增加, 氮气扩散进入油页岩层的速度增加, 热传导的速率也加快.
2) 在高压驱替作用下, 流体主要从沿油页岩层理方向形成地裂隙流出, 但是也会有少许流体从油页岩的原生孔隙流出; 在同一截面垂直于油页岩层理方向压力损失显著, 这是由于油页岩的原生孔隙度较低, 渗流阻力较大.
3) 由于存在传热损失, 在注热初期裂隙升温缓慢, 油页岩难以达到裂解温度, 孔隙度低, 导流能力差, 导致地层升温缓慢, 因此可以尝试增加氮气流量, 强化传热, 提高地层升温速率.
参考文献

[1]	Martins M F, Salvador S, Thovert J F, et al. Co-current combustion of oil shale--Part 1.Characterization of the solid and gaseous products[J].Fuel, 2010, 89: 144–151.DOI:10.1016/j.fuel.2009.06.036
[2]	Na J G, Im C H, Chung S H, et al. Effect of oil shale retorting temperature on shale oil yield and properties[J].Fuel, 2012, 95: 131–135.DOI:10.1016/j.fuel.2011.11.029
[3]	Guo H F, Peng S Y, Lin J D, et al. Retorting oil shale by a self-heating route[J].Energy & Fuels, 2013, 27: 2445–2451.
[4]	Jiang X M, Han X X, Cui Z G. New technology for the comprehensive utilization of Chinese oil shale resources[J].Energy, 2007, 32: 772–777.DOI:10.1016/j.energy.2006.05.001
[5]	姜鹏飞, 孙友宏, 郭威, 等. 压裂-注氮原位裂解油页岩加热工艺及传热模拟[J].东北大学学报(自然科学版), 2015, 36(9): 1353–1357. ( Jiang Peng-fei, Sun You-hong, Guo Wei, et al. Heating technology and heat transfer simulation for oil shale of in-situ pyrolysis by fracturing and nitrogen injection[J].Journal of Northeastern University(Natural Science), 2015, 36(9): 1353–1357.DOI:10.3969/j.issn.1005-3026.2015.09.029)
[6]	王永辉, 卢拥军, 李永平, 等. 非常规储层压裂改造技术进展及应用[J].石油学报, 2012, 33(sup 1): 149–158. ( Wang Yong-hui, Lu Yong-jun, Li Yong-ping, et al. Progress and application of hydraulic fracturing technology in unconventional reservoir[J].Acta Petrolei Sinica, 2012, 33(sup 1): 149–158.)
[7]	Ma C C, Jiang Y P, Xing H L, et al. Numerical modelling of fracturing effect stimulated by pulsating hydraulic fracturing in coal seam gas reservoir[J].Journal of Natural Gas Science and Engineering, 2017, 46: 651–663.DOI:10.1016/j.jngse.2017.08.016
[8]	Zhang J C, Bian X B. Numerical simulation of hydraulic fracturing coalbed methane reservoir with independent fracture grid[J].Fuel, 2015, 143: 543–546.DOI:10.1016/j.fuel.2014.11.070
[9]	Tang H Y, Winterfeld P H, Wu Y S, et al. Integrated simulation of multi-stage hydraulic fracturing in unconventional reservoirs[J].Journal of Natural Gas Science and Engineering, 2016, 36: 875–892.DOI:10.1016/j.jngse.2016.11.018
[10]	Beskok A, Karniadakis G E. Report:a model for flows in channels, pipes, and ducts at micro and nano scales[J].Microscale Thermophysical Engineering, 1999, 3(1): 43–77.DOI:10.1080/108939599199864
[11]	Bikic S, Bukurov M, Ta?in S, et al. Determination of sucrose solution viscosity on evaporator by Hagen-Poiseuille equation[J].Journal on Processing and Agriculture, 2010, 14(4): 186–190.
[12]	王补宣. 多孔介质中单相对流换热分析的流体渗流模式[J].上海交通大学学报, 1999, 38(8): 966–969. ( Wang Bu-xuan. On the modelling of fluid flow in porous media[J].Journal of Shanghai Jiaotong University, 1999, 38(8): 966–969.DOI:10.3321/j.issn:1006-2467.1999.08.012)