基于分形理论的赫巴流体在多孔介质中的渗流模型_杨仙.pdf

第 48 卷 第 3 期 煤田地质与勘探 Vol. 48 No.3 2020 年 6 月 COAL GEOLOGY 2. Xi’an Research Institute Co. Ltd., China Coal Technology and Engineering Group Corp., Xi’an 710077, China; 3. School of Geosciences and Info-physics, Central South University, Changsha 410000, China; 4. School of Ination and Electrical Engineering, Hunan University of Science and Technology, Xiangtan 411201, China; 5. Chongqing Engineering Co. Ltd., China Coal Technology Hershel-Bulkley fluid; porous medium; fractal theory; seepage model 泥浆是一种用途广泛的工程浆液，其在钻井工 程、顶管工程、钻孔灌注桩和地下连续墙等施工中 都是必不可少的[1-2]。工程施工过程中，泥浆会渗流 进入地层，给不同的工程带来有利或不利影响。建 立泥浆在地层中的渗流模型，深入研究其渗流机 理，有利于在施工过程中更好地控制渗流过程， ChaoXing 第 3 期 杨仙等 基于分形理论的赫巴流体在多孔介质中的渗流模型 123 增加渗流的有利影响，减小其不利影响。要建立 泥浆渗流模型，首先需要选择泥浆流变模型。当 前，研究者们常采用宾汉流体和幂律流体描述泥 浆的流变性能[3-4]。宾汉流体多用于表征泥浆的塑 性特征，幂律流体多用于表征泥浆的假塑性特征， 宾汉流体和幂律流体均可看作赫巴流体的一种特 殊形式。赫巴模型既能反映泥浆的塑性特征，又 能反映其假塑性特征[5]，比宾汉模式和幂律模式 更全面，对于塑性和假塑性均较为明显的泥浆， 采取赫巴流体进行模拟更为合适，准确率最高[6]。 分形理论是数学家本华曼德博首先提出的，即用 分形分维的数学工具来描述研究客观事物[7]。众 多研究者指出岩土介质适宜于使用分形理论来进 行研究C. Sparrow[8]指出天然多孔岩土体往往具 有良好的分形特征；D. Avnir 等[9]提出了岩土多 孔介质具有分形特征的概念；J. P. Hansen 等[10]用 盒子法测得了砂岩的分形维数；A. J. Katz 等[11] 和 C. E. Krohn [12]基于 SEM 实验，结合分形理论 表征出砂岩的孔隙结构；郁伯铭团队较全面地论 述了多孔介质输运性质的分形分析理论[13-17]；员 美娟等 [18]基于服从分形分布的弯曲毛细管束模 型，运用分形几何理论推导出卡森流体在多孔介 质中流动的流量、流速、启动压力梯度和有效渗 透率的分形解析解；万秀峰[19]在研究砂卵石地层 水泥–水玻璃复合注浆加固时，基于分形理论，把 浆液简化为宾汉流体，推导其在地层中的渗透模 型；周子龙等[20]根据分形理论推导出孔隙通道的 曲折效应方程，并通过幂律流体本构方程导出考 虑孔隙曲折效应的浆液扩散模型。由上述研究可 以看出，已有研究基于分形理论，把工程浆液模 拟为宾汉流体和幂律流体，研究其在岩土体中的 渗流机理。考虑到对于某些塑性和假塑性均较为 明显的泥浆，采取赫巴流体描述流体特性更全面、 更准确，本研究中把地层岩土体简化为多孔介质， 泥浆模拟为赫巴流体，基于分形理论建立了赫巴 流体在多孔介质中的渗流模型，并对其计算结果 进行详细分析。 1 赫巴流体在单根圆直管中的渗流模型 要研究赫巴流体在多孔介质中的渗流模型，首 先，需要研究赫巴流体在圆直管中的渗流模型。赫 巴流体在圆直管中的渗流情况如图1所示。 图 1 赫巴流体在圆直管中的渗流示意 Fig.1 Schematic seepage diagram of Hershel-Bulkley fluid in round straight tube 在图1中， 针对流体单元柱元素， 有如下表达式 d 2 d rp L   1 式中p为注浆压力；L为流体流经的距离；τ为剪切 应力；r为径向距离；dp/dL为压力梯度。 赫巴流体本构方程 0 n k 2 由式1和式2得到 1 1 0 0 d1d d2 d n n vrp rkkL        3 式中为速度梯度；n为流性指数；k为稠度系数； τ0为屈服应力；v为流速。 当0≤r≤rp时，每一层流体相对于邻层，保持 静止状态，速度均相等。当rpr≤r0时，每一层流体 相对于邻层，保持运动状态，速度不相等。 当0时，即可求出 0 p 2d d L r p    4 假设赫巴流体的流动方式为层流，则边界条件 为当r r0 时，v 0。 当rpr≤r0时 0 1 0 1d d 2 d rn r rp vr kL       5 结合边界条件及式4和式5，可得 111 0pp d 1 2 d nn nnn np vrrr r nk L        - ） - ） 6 当0≤r≤rp时 11 p0p d 1 2 d n nn np vrr nk L     - ） 7 假设 1 d 1 2 d n np A nk L    8 圆管中，赫巴流体的流量是通过相对静止区0≤ r≤rp和相对运动区rpr≤r0两部分的流量总和 ChaoXing 124 煤田地质与勘探 第 48 卷 0 p 0 p 2 pp 111 2 0pp0pp 13121 2 0p00pp0p π2π d π2πd 2 π2 π π 3121 r r nnn r nnn r nnn nnn qr vrv r A rrrrArrr rr n An A A rrrrrr rr nn            - ） - ） - ） - ）- - ） - ） 9 式中r0为管道半径；rp为塞流层半径；vp为塞流 层流速；q 为单管流量。 2 多孔介质分形理论基础 根据分形理论，多孔介质孔隙大小分布在一定 范围内，且在统计上满足分形标度关系[21] f D M LL（ ） 10 式中L表示尺度；ML表示物体的体积、面积、质 量或曲线长度等指标；Df为分形维数，取值02。分 形物体累计数N与颗粒大小的分布也满足如下的标 度关系[21] f max D R N LR R     ≥ 11 式中Rmax为最大孔隙尺寸；R为孔隙尺寸。 对式11进行微分，得到 ff 1 fmax dd DD ND RRR    12 多孔介质的随机弯曲通道长 度LT可根据分形理 论得出[17] TT 1 T0 2 DD LLR   13 式中L0为通道的直线长度；DT为迂曲度分形维数； DT 1时表示毛细通道是直线；DT取值范围为12。 从式13可以看出，LT和半径R有关，R越大，则LT 越小。 对式13微分可得到 TT 11 T0T0 d2 d DD LLRDL   14 孔隙分形维数计算Df公式如下[17] f minmax ln 2 ln/ D RR   15 多孔介质迂曲度分形维数DT计算公式如下[17-20] T 0av ln 1 ln/2 T D LR 16 流线迂曲度的表达式为[17-20] 2 11 1 41 11 111 2211 T               17 式中T为流线迂曲度；为多孔介质孔隙率；Rav 为多孔介质中平均孔隙半径；Rav的计算公式为[17] 1 f 4 avmax f 4 D RR D 18 多孔介质的结构参数计算公式为[17] 0 2π 31 LR  ， 19 max 22π 2 4131 R R         ， 20 式中R ，为多孔介质的平均粒径。 多孔介质的结构参数对于流体在多孔介质中的 渗流影响很大。如，单个孔隙大小直接影响渗透率 和渗流速度，小孔体积占岩样孔隙总体积的比例越 大，其渗透率和渗流速度就越小。根据毛管模型理 论与达西定律的关系可知，渗透率的平方根与平均 孔隙半径成正比[22]。在分形理论中，充分考虑了多 孔介质的结构参数问题，因此，新建立的渗流模型 能够更好地表征实际渗流情况。 3 赫巴流体在多孔介质中的渗流模型 在多孔介质中，流体实际渗流路径有一定的随 机性，可看作是弯曲管道中的渗流，如图2所示。 图 2 赫巴流体在多孔介质中的渗流示意 Fig.2 Schematic seepage diagram of Hershel-Bulkley fluid in porous media 因此，在式8中用LT代替L，并假设 1 T d 1 2 d n np B nk L    21 则 13121 2 p0ppp 2 π2 π π 3121 nnn nnn n Bn B q RB RrrRrrRr nn    22 3.1 赫巴流体在多孔介质中的流量 由于多孔介质中弯曲的管道大小不一，且在统 计学上服从分形幂规律，因此，多孔介质横截面的 总流量可以通过积分计算得到 max min d R R Qq RN  23 ChaoXing 第3期 杨仙等 基于分形理论的赫巴流体在多孔介质中的渗流模型 125 式中Q为多孔介质横截面的总流量。 f Tf TT 1 fmax 2 0T π π dd 1 2 D DD n DD n D Rp BNRR n kLD      24 假设 f TT 1 fmax 2 0T π 1 2 D n DD n D Rp C n kLD    25 则 maxmax TfTf minmin max Tf min 131 2 pp 21 p p 2 dd 31 2 d 21 nn RR DDDD nn RR n R DD n R nC QC R rRRR rRR n nCr R rRR n          - ）- - ）- - ） 26 该积分没有解析解，可使用MATLAB计算其数 值解。 3.2 赫巴流体在多孔介质中的流速 多孔介质横截面积的计算公式如下 2 fmax f π1 2 D R S D     27 式中S为多孔介质的横截面积。 用式26除以式27即可得到赫巴流体在多孔 介质中的平均速度。 Q v S ， 28 式中v，为赫巴流体在多孔介质中的平均速度。 4 模型计算结果与分析 4.1 流速与压力梯度之间的关系 图3为根据式28计算出来的瞬时平均流速随 压力梯度的变化关系。计算过程中，选取赫巴流体 本构方程相关的参数为n0.6，k0.9，τ010 Pa；选 取多孔介质的相关参数有0.6，Rmax0.05 m， Rmin0.000 5 m。为了便于更直观地了解本次所采用 泥浆的流变性能，采用六速旋转黏度仪对浆液的剪 切应力进行了实测。实验室配置出泥浆后，用六速 旋转黏度仪进行测试，当旋转速度为600 r/min 1 022 s–1时，赫巴流体的测试剪切应力为66 Pa，根 据式2计算的理论剪切应力值为67.5 Pa，与测定值 吻合良好。因式2可直接计算出泥浆，且该流体配 置至完全符合所选择的计算参数较为困难，因此， 本次实验室只对参数n0.6，k0.9，τ010 Pa时的配 置泥浆进行测定，其他泥浆的剪切应力值均直接由 式2计算得出。 图 3 流速随压力梯度的变化关系 Fig.3 The relationship between velocity and pressure gradient 由图3可以看出，压力梯度与流速之间的关系符 合幂指数关系， 随着压力梯度的增加， 流速逐步增加， 且流速增加的趋势越来越明显，这与实际工程中的状 况是一致的。本模型计算的是单位时间的流量和瞬时 流速，因此，把压力梯度作为一个常量代入式28。 但在实际工程中，如注浆等施工时，随着浆液从注浆 口往多孔介质弯曲通道中的扩散时间和扩散距离的增 加，压力梯度是逐步降低的。因此，如果涉及到注浆 时间累积的计算时， 压力梯度不能当作常量进行计算。 4.2 流速与 n、、k 值之间的关系 当赫巴流体本构方程的相关参数k0.9，τ010 Pa时，变化n值，通过式2计算其剪切应力。当旋 转速度为600 r/min1 022 s–1时，不同n值对应的赫 巴流体剪切应力见表1。 表 1 n 值变化时赫巴流体剪切应力计算值 Table 1 Calculation values of apparent viscosity of Hershel-Bulkley fluid with the change of n n 0.550.60 0.65 0.70 0.75 剪切应力/ Pa 50.767.5 91.4 125.0 172.7 图4为根据式28计算出来的流速随n值的变化 关系。除赫巴流体本构的相关参数外，其余选取的 计算参数有压力梯度dp/dL200 kPa/m，多孔介质 相关参数有0.6，rmax0.05 m，rmin0.000 5 m。 当赫巴流体本构方程的相关参数n0.6，τ010 Pa时，变化k值，通过式2计算其剪切应力。当旋转 速度为600 r/min1 022 s–1时，不同k值对应的赫巴 流体剪切应力见表2。 图5为根据式28计算出来的流速随k值的变化 关系。除赫巴流体本构方程的相关参数外，其余选 ChaoXing 126 煤田地质与勘探 第48卷 取的计算参数有压力梯度dp/dL200 kPa/m，多孔 介质相关参数有0.6，Rmax0.05 m，Rmin0.000 5 m。 图 4 流速随流性指数 n 值的变化关系 Fig.4 The relationship between velocity and n 表 2 k 值变化时赫巴流体剪切应力计算值 Table 2 Calculation values of apparent viscosity of Hershel-Bulkley fluid with the change of k k 0.50 0.60 0.70 0.80 0.90 剪切应力/ Pa 41.9 48.3 54.7 61.1 67.5 图 5 流速随稠度系数 k 值的变化关系 Fig.5 The relationship between velocity and k 由图4、图5可以看出，流速与赫巴流体本构参 数n、k值之间的关系均呈幂指数关系。结合表1、表 2可知，n、k值越大，剪切应力越大，流速越小。n 值变化比k值变化对剪切应力影响更大，n值变化对 流速的影响也更大。结合图4、图5中曲线趋势及 式2可知当n值不变、k值变化，或k值不变、n值 变化，使得两种赫巴流体的剪切应力相同时，n值小 的流体， 其流速会远大于k值更小的流体。 由此可见， 对于赫巴流体流速而言，n值的变化是一个极敏感因 素。由此可知，要更好地控制泥浆在多孔介质中的 渗流，简单地控制黏度不一定能达到预期的目的， 更为科学的方法是综合控制n、k值，特别是选取合 适的n值范围，对流体渗流速度影响很大。 4.3 流速与孔隙率之间的关系 图6为根据式28计算出来的流速随孔隙率的 变化关系。选取的计算参数有赫巴流体本构方程 的相关参数n0.6，k0.9，τ010 Pa；计算选取的多 孔介质相关参数有Rmax0.05 m，Rmin0.000 5 m；压 力梯度dp/dL200 kPa/m。 图 6 流速随孔隙率的变化关系 Fig.6 The relationship between velocity and porosity 由图6可以看出，流速与多孔介质孔隙率之间 的关系呈二项式关系，孔隙率越大，流速越大，且 随着孔隙率的增加，流速增加的趋势越来越明显。 当孔隙率由30增加到45时，流速增加2.3倍；当 孔隙率由30增加到60时，流速增加7倍。这也解 释了在实际工程中孔隙率对渗流影响极大的深层 原因。 本节中根据所建立的赫巴流体在多孔介质中的 渗流模型，着重对施工参数压力梯度、流体特性参 数n、k值以及岩土体特性参数孔隙率等对瞬时渗流 平均流速的影响进行了计算与分析，计算结果与实际 工作状况相吻合。基于建立的渗流模型，还可以进一 步分析流体的启动压力梯度、有效渗透率等与各类参 数之间的关系；引入时间因子，还能进一步建立渗流 扩散模型，研究泥浆在岩土体中的扩散半径、扩散过 程中压力的衰减以及渗流全程平均流速变化规律等。 5 结 论 a. 把地层简化为多孔介质，泥浆模拟为赫巴流 体，基于分形理论，建立了赫巴流体在多孔介质中 的渗流模型，并对渗流模型的计算结果进行了分析 研究。研究结果显示，压力梯度、赫巴流体本构参 数及地层孔隙率均会对渗流流速产生影响。 b. 压力梯度与流速之间的关系符合幂指数关系， 随着压力梯度的增加，流速逐步增加，且流速增加的 ChaoXing 第3期 杨仙等 基于分形理论的赫巴流体在多孔介质中的渗流模型 127 趋势越来越明显；流速与赫巴流体本构参数n、k值之 间的关系均呈幂指数关系，n、k值越大，剪切应力越 大，流速越小，且n为一个极敏感因素；流速与多孔介 质孔隙率之间的关系呈二项式关系，孔隙率越大，流 速越大，且流速增加的趋势越来越明显。 c. 为了更好地控制泥浆在地层中的渗流，可以 基于建立的渗流模型，在施工参数及泥浆配置方面 提出一些措施。且泥浆配置过程中，需综合控制n、 k值，特别是选取合适的n值范围，简单地控制流体 黏度可能达不到预期目的。 d. 基于建立的渗流模型，引入时间因子，则可 进一步建立渗流扩散模型，对渗流扩散半径、渗流 扩散过程中的压力衰减、渗流全程平均流速等进行 更为全面深入的研究。 请听作者语音介绍创新技术成果 等信息，欢迎与作者进行交流 参考文献References OSID 码 [1] 莫世扬，杨晓伟，洪元堂，等. 非开挖顶管工艺在公路污水管 线下穿工程中应用分析[J]. 公路工程，2019，442151–155. MO Shiyang，YANG Xiaowei，HONG Yuantang，et al. Ap- plication analysis of trenchless pipe jacking technology in high- way sewage pipeline project[J]. Highway Engineering，2019， 442151–155. [2] 邓利蓉，刘振兴，曹钧，等. 硬岩破碎地层高温钻孔钻进的钻 井液实验研究[J]. 煤田地质与勘探，2015，434117–119. DENG Lirong，LIU Zhenxing，CAO Jun，et al. Experiment of drilling fluid used to drill hard fractured complex ation[J]. Coal Geology Exploration，2015，434117–119. [3] 侯璐瑶， 王奎升， 侯世全. 基于宾汉模式的桥梁桩基钻孔泥浆 流变特性初步研究[J]. 铁道建筑，2010641–43. HOU Luyao，WANG Kuisheng，HOU Shiquan. Preliminary study on rheological properties of drilling mud from Bingham model to bridge pile foundation[J]. Railway Engineering， 2010641–43. [4] 蔡亮学. 水平定向钻管道穿越回拖过程动态力学特性研究[D]. 东营中国石油大学华东，2011. CAI Liangxue. Investigation on dynamics of pipe during pull- back in horizontal directional drilling installations[D]. Dongy- ingChina University of PetroleumEast China，2011. [5] 汪友平，郑秀华，夏柏如，等. 圆管中赫巴流体层流流动规律 的研究[J]. 钻井液与完井液，2008，25127–28. WANG Youping，ZHENG Xiuhua，XIA Bairu，et al. Studies of the laminar flow of Herschel-Bulkley fluids in pipes[J]. Drilling Fluid and Completion Fluid，2008，25127–28. [6] 李鑫，高德利，刁斌斌，等. 基于赫巴流体的页岩气大位移水 平井裸眼延伸极限分析[J]. 天然气工业， 2016， 3610 85–92. LI Xin，GAO Deli，DIAO Binbin，et al. Analysis on the open-hole extension limit of a shale-gas extended-reach horizon- tal well based on Hershel-Bulkley fluids[J]. Natural Gas Indus- try，2016，361085–92. [7] MANDELBROT B B. How long is the coast of britain statisti- cal self-similarity and fractional dimension[J]. Science，1967， 1563775636–638. [8] SPARROW C. The fractal geometry of nature by B. Mandel- brot[J]. Journal of the Royal Statistical Society，1984，1474 616–618. [9] AVNIR D，FARIN D，PFEIFER P. Chemistry in non-integer dimensions between two and three. II. Fractal surfaces of ad- sorbents[J]. The Journal of Chemical Physics，1983，797 3566–3571. [10] HANSEN J P，SKJELTORP A T. Fractal pore space and rock permeability implications[J]. Physical Review B Condensed Matter，1988，3842635–2638. [11] KATZ A J，THOMPSON A H. Fractal sandstone pores Implications for conductivity and pore ation[J]. Physical Review Letters，1985，54121325–1328. [12] KROHN C E. Sandstone fractal and euclidean pore volume dis- tributions[J]. Journal of Geophysical ResearchSolid Earth， 1988，93B43286–3296. [13] 吴金随. 多孔介质里流动阻力分析[D]. 武汉华中科技大学， 2007. WU Jinsui. Analysis of resistance for flow through porous me- dia[D]. WuhanHuazhong University of Science and Technol- ogy，2007. [14] 邹明清. 分形理论的若干应用[D]. 武汉 华中科技大学， 2007. ZOU Mingqing. Fractal theory and its applications to porous media, rough surface and thermal contact conductance[D]. Wu- hanHuazhong University of Science and Technology，2007. [15] YU B M. Fractal character for tortuous stream tubes in porous media[J]. Chinese Physics Letters，2005，221158–160. [16] XU P，YU B M. Developing a new of permeability and Kozeny-Carman constant for homogeneous porous media by means of fractal geometry[J]. Advances in Water Resources， 2008，31174–81. [17] LI J H，YU B M. Tortuosity of flow paths through a Sierpinski carpet[J]. Chinese Physics Letters，2011，283117–119. [18] 员美娟，郁伯铭，郑伟，等. 多孔介质中卡森流体的分形分 析[J]. 物理学报，2011，602410–415. YUAN Meijuan，YU Boming，ZHENG Wei，et al. Fractal analysis of Casson fluid flow in porous media[J]. Acta Physica Sinica，2011，602410–415. [19] 万秀峰. 砂卵石地层水泥–水玻璃复合注浆加固关键技术研 究[D]. 长沙中南大学，2014. WAN Xiufeng. Gordian technique of cement and sodium silicate composite grouting reinforcement in sand gravel strata[D]. ChangshaCentral South University，2014. [20] 周子龙，杜雪明，陈钊，等. 考虑孔隙曲折效应的浆液扩散压 力[J]. 中国有色金属学报，2016，2681721–1727. ZHOU Zilong，DU Xueming，CHEN Zhao，et al. Grout dis- persion considering effect of pore tortuosity[J]. The Chinese Journal of Nonferrous Metals，2016，2681721–1727. [21] CRAWFORD J W，MATSUI N. Heterogeneity of the pore and solid volume of soilDistinguishing a fractal space from its non–fractal complement[J]. Geoderma， 1996， 733/4 183–195. [22] 黄延章. 低渗透油层渗流机理[M]. 北京石油工业出版社， 1998. HUANG Yanzhang. Seepage mechanism of the low permeability reservoir[M]. BeijingPetroleum Industry Press，1998. 责任编辑 聂爱兰 ChaoXing