（固体力学专业论文）均匀岩层水力压裂的数值模拟及机理研究.pdf

中国科学技术大学硕士学位论文 摘要 在中国目前探明的油气储量中，低渗透油田占了很大的比例。目前制约低渗 透储层开发的关键问题是已掌握的开发技术难以使此类油田得到经济高效地开 发。虽然对低渗透油田的开发采取了一些与之相适应的配套工艺技术，也取得了 一定成绩，但是，由于对影响开发效果的一些关键因素认识不清，且缺乏系统综 合的研究，因而导致此类油藏的开发具有一定的风险性和盲目性，同时由于低渗 透油田天然裂缝的存在，其开发特征不同于常规油藏，因此研究低渗透油藏的开 发技术至关重要。 低渗透油藏的开发主要通过水力压裂提高单井产能。为了研究岩石内部微裂 纹在水力压裂过程中的扩张及相关渗流行为，需要对岩石的水力压裂过程进行数 值模拟。本文中的数值模拟采用渗流应力耦合模型来模拟岩石内部的力学和渗流 行为，并采用粘结单元来模拟由水力损伤造成的预设扩展裂缝，在此基础上提出 了模拟三维水力压裂的有限元模型和模拟主裂纹和天然裂纹相互作用的有限元 模型。为了保证在裂纹扩张过程中，流体压力能随裂缝的扩张动态地加载到裂纹 面上，我们给出了流体压力在裂纹内传递的模型。 基于以上建立的有限元数值模型，我们对垂直裂纹和水平裂纹进行了三维数 值模拟，对主裂纹和天然裂纹的相互作用进行了二维数值模拟，并得到了整个水 力压裂过程中岩石内部的应力分布，孔隙压力分布，压裂液的滤失以及裂缝的几 何尺寸。随后我们进一步分析了获得的结果来研究水力压裂过程的力学机理。分 析揭示了裂纹的最终几何尺寸与压裂液的压力梯度大小和隔层影响的关系，主裂 纹与天然裂纹的相互作用和两者的角度有很大的关系。本次数值模拟的结果对解 释和研究岩石水力压裂过程的力学机理有很重要的作用并且对石油工程中压裂 作业有现实的指导意义。 中国科学技术大学硕士学位论文 a b a s t r a c t t h er e s e r v e so fl o w - p e r m e a b l eo i lw e l lt a k e sag r e a tp r o p o r t i o ni nt h ew h o l e r e s e r v e sd i s c o v e r e du pt on o w a tp r e s e n tt h ec r i t i c a lp r o b l e mi st h a tp r e s e n t d e v e l o p m e n tt e c h n i q u e f o ro i lw e l li s i n c a p a b l eo fm a k i n go i lw e l ld e v e l o p e d e c o n o m i c a l l ya n de f f i c i e n t l y , w h i c hl i m i tt h ed e v e l o p m e n to fl o w - p e r m e a b l eo i lw e l l a l t h o u g hs o m er e l a t e dt e c h n i q u eh a sb e e na d o p t e da n ds o m ea c h i e v e m e n tw a sm a d e ， t h ed e v e l o p m e n to fl o w - p e r m e a b l eo i lw e l li sf u l lo fr i s k sa n db l i n d n e s sb e c a u s eo f l a c ko fu n d e r s t a n d i n go ft h ek e yf a c t o r sw h i c hh a v eg r e a ti n f l u e n c eo nt h e d e v e l o p m e n te f f e c t sa n ds y s t e m i ca n ds y n t h e t i cr e s e a r c h a tt h es a m et i m e ，n a t l l r a l 丘a c 眦se x i s ti nt h eo i lw e l la n dt h ed e v e l o p m e n tt e c h n i q u ei sd i f f e r e n tf r o mn o r m a l o i lw e l ls ot h ed e v e l o p m e n tt e c h n i q u ef o rr e s e a r c ho nl o w - p e r m e a b l e 胡w e l li so f g r e a ti m p o r t a n c e h y d r a u l i cf r a c t u r i i l gi su s e dc o m m o n l y t oi m p r o v et h ep r o d u c t i o no f t h eo i lw e l l f i n i t ee l e m e n ts i m u l a t i o n sw e r ec a r r i e do u tt oi n v e s t i g a t et h eh y d r a u l i cp r o g r e s s i v e d a m a g ea n da s s o c i a t e df l o wb e h a v i o ri nr o c k t h ec o u p l e dp o r ep r e s s u r ea n ds t r e s s m o d e lw a su s e dt os i m u l a t et h er o c km e d i aa n dc o h e s i v ee l e m e n t sw e r eu s e dt o s i m u l a t et h ep r e - d e f i n e dh y d r a u l i c 丘a e t u r e t h e nn u m e r i c a lm o d e l sf o rc o u p l e dp o r e f l u i df l o wa n ds t r e s sm o d e lw a sp r o p o s e d af l u i dp r o p a g a t i o nm o d e li nf l a c t u r ew a s p r o p o s e dt ou p d a t ea n dr e d i s t r i b u t e dt h em o v i n gf l u i dl o a d i n g si nf r a c t u r e w i t ht h e p r o p o s e dc o u p l i n gm o d e l ，t h et h r e ed i m e n s i o n a lh y d r a u l i cf r a c t u r i n gp r o c e s sa n d t h e b a t u r a lf r a c t u r e s i n f l u e n c e0 nt h em a i nf r a c t u r ew e r es i m u l a t e db yc o m m e r c i a l e n g i n e e r i n gs o f t w a r ea b a q u s t h e nt h es t r e s sd i s t r i b u t i o n , t h ep o r ep r e s s u r e d i s t r i b u t i o n , t h ef l u i dl o s sa n dt h eg e o m e t r yo ft h ep r o g r e s s i v ed a m a g ei n f o r m a t i o n w e r eg a t h e r e da n ds t u d i e d t h er e s u l t ss h o wt h a tt h el e n g t ha n dt h ew i d t ho ft h ep a t h o ft h ep r o g r e s s i v ed a m a g ea r es t r o n g l yi n f l u e n c e db yb o t ht h eh y d r a u l i cp r e s s u r ea n d t h ei n j e c t i o nt i m e t h er e s u l t sp r o v i d eg o o di n t e r p r e t a t i o na n du n d e r s t a n d i n go ft h e m e c h a n i s mo fh y d r a u l i cp r o g r e s s i v ed a m a g ei nr o c k t h i ss t u d yi sv e r yu s e f u la n d i m p o r t a n tt ot h eo i le n g i n e e r i n ga n ds o m eo t h e rr o c ke n g i n e e r i n gf i e l d s 中国科学技术大学学位论文相关声明 本人声明所呈交的学位论文，是本人在导师指导下进行研究 工作所取得的成果。除己特别加以标注和致谢的地方外，论文中 不包含任何他人已经发表或撰写过的研究成果。与我一同工作的 同志对本研究所做的贡献均已在论文中作了明确的说明。 本人授权中国科学技术大学拥有学位论文的部分使用权， 即：学校有权按有关规定向国家有关部门或机构送交论文的复印 件和电子版，允许论文被查阅和借阅，可以将学位论文编入有关 数据库进行检索，可以采用影印、缩e 乜或扫描等复制手段保存、 汇编学位论文。 保密的学位论文在解密后也遵守此规定。 作者签名：蕴烟 年月日 中国科学技术大学硕士学位论文 1 1 水力压裂技术介绍 第一章绪论 岩体是由多种矿物晶粒、胶结物及孔隙缺陷等组成的混合体，是自然界中经 过亿万年的地质演变和地质构造运动所形成的最为复杂的固体材料之一。岩石在 承受一定外界载荷下就会发生裂纹的起裂和扩张过程。由于岩石中本身内部原本 就存在一些随机分布的天然孔隙缺陷，所以岩石的裂纹起裂和扩张既有新孔隙缺 陷的萌生、增大、扩张、贯通的过程，也有原有缺陷发展、扩张、贯通过程。但 是深层地底情况复杂，影响岩石裂纹扩张的因素也是多样的。 目前国内众多油田的油井压裂施工普遍采用水力压裂这种方式。水力压裂是 人们利用地面高压泵组，以超过地层吸收能力的排量将高粘压裂液泵入井内而在 井底产生高压，当该压力克服井壁附近地应力并达到岩石抗张强度，就在地层产 生裂缝并向前延伸一定长度1 ”。为了使产生的裂缝在高压液体返排后不至于闭 合，在形成裂缝以后，接着泵注混有支撑剂的携砂液，携砂液将继续延伸裂缝并 将支撑剂输送到裂缝内。当支撑剂泵送完毕后，将粘性压裂液破胶降为低粘度的 液体流回井内，然后返排到地面。这样目标层内留下一条高导流的裂缝，以利地 层远处的油和气流向井底。 水力压裂作为一种应用广泛的现代工程技术，在石油开采、煤层注水、防治 冲击地压、地应力测量、水库诱发地震、水坝及边坡失稳控制、地热开发等众多 领域也具有十分广阔的应用前景删。 1 2 压裂开发的目的与意义 中国低渗透油气资源十分丰富，截至1 9 9 6 年底已探明低渗透油田共2 8 5 个， 地质储量约4 0 x1 0 ”t ，占全部探明地质储量的2 4 5 ，广泛分布于全国勘探开 发的2 1 油区。但到目前为止，陆上已开发地质储量中低渗透油气储量仅占1 7 5 ， 在已探明未动用储量中低渗透储量占6 0 8 “。近年来，新增的油气储量中，低 渗透储量的比例更大，以中国石油天然气公司为例，虽然历年探明储量幅度较大， i - 中国科学技术大学硕士学位论文 但低渗透、特低渗透比例也很大( 表卜1 ) 。因此，对己开发的低渗透油气田如 何进一步提高开发效益，对未动用的低渗透油气储量如何尽快有效的投入开发， 对保持中国石油工业持续稳定发展有着十分重要的意义。 表1 - 11 9 9 7 - 2 0 0 1 年中国石油天然气股份有限公司探明储量表 时间 1 9 9 7 1 9 9 8 1 9 9 9 2 0 0 0 2 0 0 1 探明储量( 1 0 8 t ) 5 2 4 8 3 7 4 2 4 5 低渗透、特低渗透储量所占比例( ) 7 5 7 0 6 5 4 6 6 9 注：资料来源于中国石油天然气集团公司信息和经济研究中心 目前制约低渗透储层开发的关键问题是已掌握的开发技术难以使此类油田 得到经济高效地开发，表现为投入开发的大部分低渗透油田采油速度低、采收率 低，随着含水的上升，采油指数和采液指数下降等。虽然对低渗透油田的开发采 取了一些与之相适应的配套工艺技术，也取得了一定成绩，但是，由于对影响开 发效果的一些关键因素认识不清，且缺乏系统综合的研究，因而导致此类油藏的 开发具有一定的风险性和盲目性，同时由于低渗透油田天然裂缝的存在，其开发 特征不同于常规油藏，因此研究低渗透油藏的开发技术至关重要。 低渗透油田储层物性差，非均质性强，开发难度大，经济效益差，因此，加 强低渗透油气藏渗流机理研究，加强开发与开采新工艺技术研究，不断提高低渗 透油气田开发水平，获得最佳经济效益，乃是当今石油工业的当务之急。 1 3 国内外研究现状 低渗透油藏的开发主要通过水力压裂提高单井产能，以注水保持地层能量， 油水井压裂后油藏内的渗流机理和动态变化在很大程度上取决于水力裂缝的参 数和方位，油井产量除与水力裂缝参数有关外，还与注水强度、井网格局等因素 有关”。低渗透油井的开发不同于常规油井的开发，关键在于水力裂缝作用的双 重性，一方面水力裂缝的存在提高了油水井的增产增注能力，另一方面却加剧了 中田科学技术大学硕士学位论文 地层的非均匀性，当裂缝方位不利时可能会导致油水井过早水淹，降低开采效益。 而水力裂缝的方位是受最小主应力控制的，对某一具体区块，最小主应力的方向 是一定的，因而水力裂缝的方位也是一定的，但是在不同的井网方位条件下，水 力裂缝的方位可能处于有利或不利方位，所产生的开发效果却是截然不同的。因 此，在低渗透油田开发中，首先要根据裂缝的方位确定出合理的井网方位，在此 基础上优化井网类型和井网密度，然后优化裂缝长度，导流能力和相应的注采压 力差，最终以经济效益为目标进行整体方案优化设计，提出合理的开发方案，指 导油田生产。 1 3 1 国内研究现状 国内外对有关低渗透油藏整体压裂改造的研究最早可以追溯到上个世纪5 0 年代，但是在理论和应用的突破还是在9 0 年代以后。 1 9 9 3 年，石油大学( 北京) 杨能宇在硕士论文“注水开发油田整体压裂改 造水力裂缝参数对采收率的影响研究h 跏中，利用电解模型和数值模拟方法研究 了五点井网整体压裂改造后，水力裂缝参数和裂缝方位对产量和采收率的影响。 在建模中采用把裂缝放大，导流能力减小的等效渗流阻力法处理裂缝系统，把油 藏和裂缝作为同一体系来考虑，使用二维两相油藏数值模拟方法。 1 9 9 7 年，石油大学( 北京) 陈志海在硕士论文“低渗透油田压注采水力裂 缝参数优化设计跏中，建立了三维两相非均匀油藏数值模拟模型和二维两相水 力裂缝模拟模型，利用裂缝和地界交界处的压力和流量相等原则实现了两个模型 之间的连接，计算时分别求解裂缝和地层内的压力和饱和度分布，更真实的反映 了裂缝的存在对低渗透油藏开采动态的影响。 从国内发表的文献看，有整体压裂改造的研究与应用研究经历了以下的发展 过程：地层内的流体渗流由仅考虑单相的纯油流发展到考虑油水两相的流动；压 裂方式由单井压裂增产发展到区块整体压裂改造。 1 3 2 国外研究现状 国外对压裂后人工裂缝对油藏动态影响的研究比较早。方法主要是实验和数 值模拟计算，从研究方法上可以明显地分为两个阶段两种类型。2 0 世纪6 0 年代 以前主要用电解模型实验方法，通过测定电位值，确定渗流场的变化，再计算扫 中国科学技术大学硕士学位论文 油效率。6 0 年代以后，随着计算机技术的发展，多采用数值模拟计算方法。具 有代表性的研究成果如下： 1 9 7 8 年，h o l d i t c h 等人l ”发表了“低渗透气藏开发中的井距和缝长的优化”， 这是最早公开讨论裂缝缝长和井距优化问题的论文之一。本文突出的优化模型主 要包括g e e r t s m a 和d ck l e r k 提出的垂直缝缝长计算模型，g r i n g g a r t e n 等人提 出的油藏动态解析油藏模型的经济模型，利用此优化模型分别对低渗透、中渗透 和高渗透气藏进行了计算，较为详细的提供了一套在低渗透气藏中对缝长和井距 的优化方法。 1 9 8 5 年，b r i t tn 1 1 在“中渗透油藏油井压裂的最优化”一文中用二维三相 模型分析了人工裂缝对五点井网单元水驱采收率的影响。结果表明：裂缝处于不 利于方位时，其半长超过井距的2 5 将对水驱采收率产生不利影响。考虑采油、 注水井都压裂时，总缝长超过井距的5 0 就会影响最终采收率。 1 9 8 8 年，d n m e e h a n 等人n ”发表了“储集层非均匀性与裂缝方位对裂缝 长度和井距优化的影响”，该文指出若裂缝长度与井距相比很小，裂缝方位就不 会产生井间干扰；并给出了考虑裂缝方位的均匀流裂缝解，证明当井间干扰很大 时，均匀流裂缝模型就不适用：文中最后用数值模拟方法证实当井距和裂缝长度 之比小于2 时，裂缝方位是很重要的。 1 9 9 1 年，k o n e p l y o v 和z a z o v s k y 1 ”在“面积注水油藏油井压裂后水驱油的 数值模拟”一文中研究了多井压裂时裂缝方位、长度、导流能力等对水驱油的影 响。在研究中考虑了压裂后投产初期油井排液和压裂液侵入带的影响。本文采用 二维两相油藏数值模拟，计算了线性布井注采中一口井和多井压裂时，裂缝处于 不同方位，对油井生产动态的影响。 1 9 9 5 年，d n m e e n h a n “1 发表的“在非均质油藏中缝长和井距的优化” 中，利用了数值模拟方法研究了地层有限封闭中心压裂气井非均质性对采气井生 产动态的影响，指出井距和水力裂缝的长度是影响低渗透非均质气藏开发的主要 因素。得出油藏的非均质性、渗透率的各向异性对非均质气藏开发具有重要影响， 并进行了定量描述。通过经济模型对井距和缝长进行了优化，对具体的区块进行 模拟计算，得出了最优的布井和压裂施工方案。 由此可见，国外研究成果可以归纳如下： 中国科学技术大学硕士学位论文 ( 1 ) 考虑了油气藏的非均质性，但以单相气藏研究居多； ( 2 ) 油藏边界考虑五点法布井或有限边界中心一口并； ( 3 ) 研究集中于高渗流气藏或低渗透气藏。 总结国内外的研究现状不难看出：在数值模拟方面，使用有限单元法模拟三 维水力压裂过程中的力学渗流耦合行为并以此研究水力压裂机理的文献很少。这 也是本文的主要工作。 1 4 渗流场与应力场的耦合分析研究进展 流固耦合现象和流固耦合问题越来越受到许多领域学者和专家的重视。1 9 2 5 年，t e r z a g h i 提出了饱和士的一维固结理论，他首先将可变形的饱和多孔介质 固结视作流动一变形的耦合问题，提出了有效应力概念，并建立了一维固结模型， 它在土力学中得到了广泛应用。t e r z a g h i ( 1 9 4 3 ) “”将他的一维固结理论推广到 三维。b i o t ( 1 9 4 1 “”，1 9 5 6 “7 1 ) 进一步研究了三维变形材料与孔隙压力的相互作 用，并在一些假设如材料为各向同性、线弹性小变形、孔隙流体不可压缩且充满 固体骨架的孔隙空间、流体流动满足d a r c y 定律的基础上，以饱和土体的总应力 和孔隙流体压力为状态变量，建立了比较完善的三维固结理论，奠定了流固耦合 理论研究的基础。b l o t “”( 1 9 5 5 ) 将其理论推广到各向异性孔隙固体的弹性固结 理论，假设孔隙中流体是可压缩的粘性流体，对固体骨架应力与孔隙流体压力、 固体骨架变形与流体运动作了区别，并将d a r c y 定律表示为流体相对于固体骨架 的运动，该理论可以预测流体渗流、各向异性介质中应力一应变随时间的变化历 史。z i e n k i e w i c z j ( 1 9 8 4 ) 考虑了几何和材料的非线性，并在b l o t 的三维固结 理论的基础上提出了广义b i o t 公式。 国内李锡夔噼1 等讨论了考虑饱和土固结效应的结构一土壤相互作用问题，张 洪武和钟万勰”1 等利用z i e n k i e w i c z 和s h i o m i 建立的广义b i o t 公式对饱和土 固结非线性问题的理论和算法进行了研究。陈平、张有天嘲，黎水泉和徐秉业伽 以及仵彦卿、柴军瑞嘲等对裂隙岩体渗流与应力耦合问题进行了研究。冉启全 建立了弹塑性交形油藏中多相渗流的流固耦合模型，并用有限差分法求解渗流方 程，有限元法求解应力方程。董平川建立了描述变形饱和储层流体流动的数学 模型，并用有限元法进行求解。薛世峰d ”建立了非混溶饱和两相渗流与孔隙介质 中国科学技术大学硕士学位论文 耦合作用的数学模型，并用流体和固体方程解耦分析方法进行求解。范学平进行 了低渗透压裂气藏、油藏压裂后流固耦合数学模拟研究，用有限差分法分别对流 体和固体方程进行求解。 近年来，随着试验测试和计算技术的发展，流固耦合问题从理论到应用都有 了长足发展。其主要特点是：由简单的单相孔隙介质模型向更复杂的双相( 孔隙 一裂隙) 连续介质及拟连续或非连续的裂隙网络介质模型发展：新的数学理论方 法进一步改善了非线性有限变形体本构模型( 粘弹性、弹塑性、蠕变) 的求解方 法。 1 5 本文的主要工作 本文采用有限元方法模拟并分析了了水力压力作用下裂缝的扩展过程。论文 的主要内容有： 在第二章中，从有限元数值模拟这一角度出发，详细阐述了水力压裂数值模 拟涉及到的相关技术细节。包括孔睬介质的有效应力原理、孔隙介质应力平衡方 程、流体连续性方程、耦合方程的求解方法、岩石孔隙介质的基本力学行为以及 用c o h e s i v e 单元模拟裂纹扩展时涉及到的相关损伤力学计算方法。 在第三章中，根据已有参考文献，获得裂缝面内水力压力分布的表达式，并 利用a b a q u s 用户子程序接口f l o w 和d i s p 编写用户子程序嵌入到a b a q u s 中，模 拟了垂直裂纹和水平裂纹的真三维水力瞬态压裂全过程，并对其结果进行了分 析。 在第四章中，模拟了平面内主裂缝遇到天然裂缝时的扩张情况，并分析和总 结了模拟结果。 在第五章中，对全文的工作进行了总结，并展望了进一步的研究内容和方向。 中国科学技术大学硕士学位论文 第二章水力压裂数值模拟中的力学理论基础 2 1 引言 水力压裂过程是一个十分复杂的力学、物理、化学过程，很难直接观察到， 只能借助于建立在种种假设和简化条件的基础上进行间接分析。数值模拟作为认 识水力压裂机制的一种重要手段，通过模拟各种地质条件下裂缝的启裂和延伸过 程，可以获得水力压裂裂缝扩展的基本规律。然而水力压裂数值模拟技术是一项 综合性的、跨学科的边缘学科，涉及到流体力学、渗流力学、岩石力学、断裂力 学、计算力学等多学科知识，本章将介绍水力压裂数值模拟涉及到的力学基础知 识。 2 2 岩石的孔隙介质模型 2 2 1 孔隙介质模型的基本概念 岩石是由固体骨架和相互连通的孔隙以及储存于骨架孔隙中的流体( 油、气 和水) 三者组成的多孔介质。岩石中固体骨架、液相和气相在数量上的比例关系 不仅可以描述岩石的物理性质和它所处的状态，而且可在一定程度上反映岩石的 力学性质。岩石的三相图是表示岩石体内三相相对含量的直方图，如图2 1 所示， 图中d 矿表示单元体积：d 巧表示固相颗粒体积；d 圪表示液相体积；k 表示孔隙体 积，它可以完全饱和( 即孔隙体积厢内完全充满液体) 或者部分饱和；某些岩 石( 如含粘土矿物颗粒的岩石) 中包含的颗粒在吸收水后会发生膨胀，有时吸收 的湿液体积相当大，不可忽略，因此可定义该部分湿液体积为k 。 中国科学技术大学硕士学位论文 d 上 岩石的孔隙率力定义为岩石中孔隙的体积与岩石颗粒的体积之比，即 疗：盟= id v g d r ,( 2 2 1 ) d矿dvd v 使用上标0 表示参照构形值，因而当前构形的孔隙率又可用下式表示： 州一冬贮善一韭 d 硭d vd v 9 d v ( 2 2 2 ) = l - j s j 。( 1 一n 。一n ? ) 一n 。 所以： 上等乌；生 ( 2 2 3 ) i - 以。一门? j 、。7 其中： j 竺l 别是介质当前构形的体积与其参考构形体积的比； 以= i 砑d v g 产 介质相体颗粒当前构形的体积与其参考构形体积的比； 鸭竺芳是单位介质当前构形中吸收的液体体积； 岩石的孔隙比e 定义为岩石中孔隙的体积与岩石颗粒的体积之比，即 8 竺黑 ( 2 2 4 )8 = ：一 ( 224 ) d r + d k 7 孔隙率盯和孔隙率p 的转换公式为： p = 三( 2 2 5 ) 1 。 一8 - 中国科学技术大学硕士学位论文 p 胛= 一l + e ( 2 2 6 ) 饱和度s 定义为岩石中孔隙水的体积与孔隙的体积之比，即 j 筌盟 d v , ( 2 2 6 ) 当完全饱和时s 的值为1 0 。而对于不含水岩石s 为0 。 某一点的自由湿液体积比为湿液体积比总体积。 竺告= 册 ( 2 2 7 ) 总的含水率等于单位体积中的自由湿液和吸收的液体，其表达式为： 吩2 册+ 吩 ( 2 2 8 ) 2 2 2 孔隙介质的有效应力原理 饱和岩石体中孔隙水是连续的，它与通常的静水一样，能够承担或传递压力。 人们把饱和岩石体中由孔隙水来承担或传递的应力定义为孔隙水应力甜。，通过 岩石颗粒间的接触面传递的应力称为有效应力厅。 如图2 2 所示，假定岩石中一点的总应力盯由三部分组成，它们分别是湿液 中的平均压应力，其他不湿润流体的平均压应力以及“有效应力”矿。其 中湿液的平均压应力又称为“湿液压力”。有效应力由下式定义： 矿= o + ( x u w + ( 1 - - x ) u o ) i ( 2 2 9 ) 应力分量如此表达是为了保证拉应力为正，而和t a 为压应力值，这就解 释了上述方程式中的符号。其中工是个依赖于饱和度和流固系统表面张力的因 子嘲，当岩石完全饱和时，工等于1 0 ；在不饱和岩石系统中，x 取值范围为从 0 0 到1 0 。由于缺少实验数据证明因子x 对饱和度的依赖程度，我们简单假设x 的取值等于岩石介质的饱和度。 为简化模型计算，我们假设整个模拟区域内施加到不湿润流体上的压力保持 中国科学技术大学硕士学位论文 不变且不随时间变化，并且它足够小可被忽略不计。这就要求不湿润流体在整个 介质内可以自由扩散，它的压力不能超过介质边界上施加到该不湿润流体上的 压力。介质边界上施加到不湿润流体上的压力在整个模拟期间保持不变。岩土工 程问题就是这种简化应用的典型例子，在岩土孔隙介质中空气就是不湿润流体， 它直接暴露于大气压中，具有极高的渗透能力。此外，要求被模拟区域不能大到 大气压力的重力梯度会显著改变空气压力，另外也不能发生促使大气压力发生瞬 间变化的外部事件。根据以上这些假设，足够小，以至于它在介质的变形中不 产生作用。因此可以删除公式( 2 2 9 ) 中的“。项，与它相关的加载项也可以从平 衡方程中删除。这样有效应力就可以改写成： 矿= 盯- i - 肋。i ( 2 2 1 0 ) 图2 2 有效应力示意图 在吸液膨胀的岩石系统中，我们假设有效应力由吸收的流体平均压应力和岩 石骨架的有效应力2 部分组成，它们的成分权重与吸收液体的相对体积和孔隙材 料有关，其有效应力的表达式如下： 矿= ( 1 一吩) 厅一吩历i ( 2 2 n ) 其中，厅是岩石骨架的有效应力； 磊是吸收液体内部的平均压应力； m 是吸收液体体积与总体积的比值。 1 0 - 中国科学技术大学硕士学位论文 我们假设岩土孔隙介质的基本力学响应由简单固体骨架弹性和液体的关系 组成。岩石骨架的有效应力厅定义成土壤应变历史、温度和其他状态变量的函数。 厅= 矛( s 扭a j nh i s t d i ) ，t e m p e r a t u r e ，s t a t e v a r i a b l e s ) ( 2 2 1 2 ) 2 2 3 孔隙渗流的基本行为 孔隙流体流动的基本行为不是服从d a r c y 定理，就是服从f o r c h h e i m e r 定理。 d a r c y 定律一般用于描述低速流体流动，而f o r c h h e i m e r 定律主要用于模拟涉及 更高流动速度的情况。d a r c y 定律可表达为层流条件下通过多孔介质单位面积的 湿液体积流率肼p 。与负的水力压头的梯度成正比。 s f i v w ；- k j 挲 ( 2 2 1 3 )2 -懈z j j ， 其中：毒是多孔介质的渗透系数。 口是水头，它的由下式定义： 中竺：+ 旦 ( 2 2 1 4 ) g p ， 兵中；二是某指定参考面之上的高度； g 是重力加速度的大小，它的作用方向与z 相反。 f o r c h h e i m e r 定律可表达为紊流条件下通过多孔介质单位面积的湿液体积流 率与负的水力梯度成二次相关性。 册，l + 厄) 警 ( 2 2 1 5 ) 其中：p g ，e ) 是速度系数嘲，非线性渗透性可以通过对材料孔隙比依赖性而 定义。 当渗流速度趋近0 时，f o r c h h e i m e r 定律就和d a r c y 定律接近。同样，当= 0 时，这两种定律是等价的。 渗透系数痞可以是各向同性，也可以是多孔介质饱和度和孔隙比的函数。毒 称为水力传导系数或渗透系数，它具有速度的量纲( 长度时间) 。渗透率置只与 多孔介质本身的结构特性有关删，而与单相牛顿体的特性无关，它的定义如下： 中国科学技术大学硕士学位论文 童2 i v 丽i 露 2 2 1 6 其中：v 是流体的动力粘度( 流体的动力粘度与质量的比率) 。 我们假设g 的大小和方向固定，则有： 警= 去陪) ( 2 2 1 7 ， 其中g 竺一g a z 锄是重力加速度( 我们假设密度成随位置缓慢变化) 。 在多相流动系统中渗透性主要受研究对象的饱和度和介质的孔隙比支配。我 们假设这种依赖性是可分的，因而： j ；c = 后，k ( 2 2 1 8 ) 其中：毛0 ) 提供了对饱和度的依赖性，函数七，可由用户定义，n g u y e n 和 d u r s o 嗍( 1 9 8 3 ) 观测到，在部分饱和介质稳态流动中，渗透性随j 3 改变，从而我 们取七，- - $ 3 作为默认值当七( 1 ) = 1 0 时，则七 p ) 代表了介质完全饱和时的渗 透性。 引入基本渗流定律后，质量连续性方程又可改写成下式： 工m 爱( 蜊训 一号( 怠( 砌训，) ) 础南警以p 亿2 1 9 心6 怠一v w 抓o 2 3 渗流应力耦合模型 2 3 1 应力平衡方程及其离散化 根据虚功原理，问题域中任一时刻t 考虑体积当前构型的应力平衡方程可用 下式表达： 口：6 s d v = i f 6 v d s + 侈8 v d v ( 2 3 1 ) 中雷科学技术大学硕士学位论文 其中：枷是虚速度场； 6 占竺。y m f 娑1 是虚应变率； 、盘， 仃是真( 柯西) 应力； t 是单位面积的表面外力7 夕是单位体积的体力，它包含了湿液的重力，即： 兀= ( s n + n t ) p 。g ( 2 3 。2 ) 其中；p 一是湿液的密度； g 是重力加速度，可以把它视作指向竖直方向的常量。 为了简化处理，我们考虑显式写出湿液的重力载荷，这样一来夕中的重力项 将仅仅与干孔隙介质相关。从而，我们可以把虚功方程改写为： 工口：6 占d 矿= f 8 v d s + ，6 v d 矿+ i ( s 以+ 珥) p ，g 6 v d 矿( 2 3 2 ) 其中：是除湿液之外的所有介质体力： 在有限元模型中，通过引入插值函数( 又称为形函数) ，平衡方程可近似表 达为一组有限方程式。用于表示这些变量离散方程的符号可写成变量带大写上 标，例如v ”表示节点变量v 的各插值函数之和。插值函数的值取决于材料坐标， 即由拉格朗日公式决定。 为简单起见，我们暂时仅考虑问题无内部约束( 例如不可压缩性) 的情况， 此时离散化方程完全仅由近似平衡方程组成。虚速度场可按下式进行插值： 6p=n”5v”(233) 其中：n ”( 墨) 是插值函数，它取决于材料坐标墨a 虚应变率的插值形式为： 6 占= 夕”6 v ” ( 2 3 4 ) 其中夕”是应变矩阵，最简单的形式为： p n 啊( 警 亿。s ， 中国科学技术大学硕士学位论文 从而虚功方程司被禹散化为： 6 v 工：f r d v 曲h 舢+ 工芦矿+ 工( 册+ q ) p 。鲥叼2 矗6 其中：6 v ”被假定为任意的。 上式方程等号左边与却共轭的项是内力矩阵，： ，筌：仃d 矿 ( 2 3 7 ) 同理，位于方程式等号右边的外力矩阵p ”可定义为： p ”竺”t d s + f d v + 工( s n + n t ) p n ”g d v ( 2 3 8 ) p ”项中包含了达朗贝尔力。 任意取函。非0 值，则意味着内力和外力达到平衡状态，亦即： ，”一p ”= o( 2 3 9 ) 这些离散平衡方程连同流体连续性方程共同确定了岩石孔隙介质的状态。当 使用隐式积分求解时，平衡方程最终被写成时间增量形式，哪怕是最简单的情况， 它们都是非线性问题。一般都使用n e w t o n 法求解。这就意味着系统需要j a c o b i a n 矩阵，它用于定义离散问题方程每一项基本变量的变化，对本问题而言，这些基 本变量就是节点坐标( 或者节点坐标的等价变量一位移一x ) 以及节点湿 液压力值甜g 。为此定义变量妒，它包含了节点位置和节点湿液压力值彰。 矽芒善矾舌血 ( 2 3 1 0 ) 由离散平衡( 2 3 9 ) 的变量可知，项d p x 导致j a e o b i a n 矩阵产生质量矩阵 ( 对达朗贝尔力而言) 和载荷刚度矩阵。 带湿液重力的载荷刚度项的表达式如下： 一争斛) g d 矿 ( 2 3 1 1 ) 其中：，竺l 芳l 是当前构形体积与参照构形体积比。 中国科学技术大学硕士学位论文 ；f ；：d ( 打) + 仃：妒” 矿 q 3 1 2 上式第一项中包含d ( ，口) ，它表示节点位置和孔隙流体压力变化产生的应力 变化，该项由有效应力原理和材料本构定义，它的表达式如下： d ( j o ) = d 7 ( ，( 1 一q ) 于) 一d ( ，珥露j r ) + ，( d 7 c r + 盯d 了7 ) 一d ( 加，) j ( 2 3 1 3 ) 其中：矿( ，矛) 用于表示有效应力随材料响应的变化，也就是说它的变化与 应变或者其他状态变量相关； d 口竺a s y m f 警1 是材料的自旋。 、掰 根据这种假设，岩石孔隙介质对j a c o b i a n 的贡献量为： d ，= 工i 专p ”： d 9 p ( 1 一一) 厅) 一d ( 以磊，) + 口：( d 夕”+ 2 p ”d 口) ：士+ 毒) 咄：帆川s 卜 亿3 蚴 其中：d 占竺研m f 譬1 是应变率，因此又有如下关系： 废， d j f ：n 搬( d s ) ：上：d f ( 2 3 1 5 ) 2 3 2 流体连续性方程 考虑某一体积岩石中包含固定数量的骨架固体材料，这个体积的当前构形占 用空间k 其表面积为s 它的参照构形占用的空间为矿。任一时间流过该体积 空间的湿液体积为k ，被吸附而成为胶体的湿液体积为k 。则该体积岩石中湿液 的总质量为： l p d v w + d v ，】= 工风( + 氇) d 矿 ( 2 3 t 6 ) 湿液质量随时间的变化率为： 昙( 风( + 吩) ) d 矿= 工了l 面d ( 纸( + 啊) ) d 矿( 2 3 1 7 ) 中国科学技术大学硕士学位论文 单位时l 司内j 恿过表曲进入该望l 司的湿_ 便质量为： 一上风玎。尼l ，。d s ( 2 3 1 8 ) 其中： v 。是湿液相对固相的平均速度( 即渗流速度) ； 打是表面s 的外法线。 让通过表面s 的湿液质量等于体积域矿内的液体质量变化，就可以得到湿液 质量的连续性方程： f 了l 磊d 慨( + 珥m 矿= 一l 矶甩一弋舔 ( 2 3 1 9 ) 根据散度定理，又e h 于体积可以是任意形状，这就得到了点态连续性方程： 了l 面d ( ，岛( + ) ) + 昙( 岛) = o ( 2 3 2 0 ) 上式的等价弱形式为： f 抛。；鲁( 地( + 珥) ) d 矿+ 帆丢( 风) d 矿= o ( 2 3 2 1 ) 其中：觑w 是任意的连续变量场。 式( 2 3 2 1 ) 又可以写成参照体积形式： l 砌。专昙( 哦( + 啊) ) d 矿。+ l 8 u w ( 风) d r 。= o ( 2 3 2 2 ) 这个连续性方程可以写成欧拉公式的后差分积分形式： l 砜 ( 纸( + 惕) ) t + ，t - - ( 1 ，“( + 珥) ) ， + 出l 6 j 丢( 成) ，+ 。d 矿= 。 q _ 名3 用当前体积构形表示，它又可写成： 甄 ( ，以( + 珥) ) ，。一( ，成( + 惕) ) ， d 矿 + a ，工a “。 j 丢( 风甩。忍。) l + 。d 矿= 。 2 3 2 4 这样我们可以就根据惯例认为任何与时间f 非显式相关的变量都取自时间 根据散度定理。方程( 2 3 2 4 ) 叉可再写为： 中国科学技术大学硕士学位论文 i m 龛c k + 哆( 舯卅缸象警以k3 瞄， + & 6 舅以乒v 声= o 其中为了表达方便，我们正则化了方程湿液密度的参照构形z 。 由于栉，= 肼，所以式( 2 3 2 6 ) 又可写成： f m 象( ) 卜专睦州) 一嗡职警 k 。嘲， + 址i s 舅湖v 拶= o 2 4c o h e s i v e 单元的损伤力学模型 c o h e s i v e 单元用于绑定两个间距可以忽略的两个界面。这种情况下 c o h e s i v e 单元的本构采用分离张开模式，即只考虑法向和两个横向剪切方向的 应力应变关系。 在c o h e s i v e 单元发生破坏之前，其力学响应表现为线弹性，其应力应变关 系如下； 仁仨 = 誊差冬 仨 = 趾 其中毛= 鲁，岛= 砉，局= 砉，e 一，e - 和e t 分别为法向和两个横向剪切方向 的应变。6 。6 。和5 。为法向和两个横向剪切方向的张开位移，t 口为c o h e s i v e 单元的初始厚度。 当c o h e s i v e 单元的应力到达一定的l 晦界值的时候，c o h e s i v e 单元发生破坏。 c o h e s i v e 单元的破坏模型主要包括损伤判断准则和损伤演化方程。 在本数值模拟中我们采用最大应力判断准则，其表达式如下： 中国科学技术大学硕士学位论文 临界应力损伤判断准则： m a x 悟，毒，豺乩 c r n 为法向应力，吒，吼( 二维的情况下不存在) 为切向应力 弹性模量线性损伤演化准则：e = ( 1 一d ) e o ， 。= 黼， e o 是无损伤单元的杨氏模量，e 是损伤单元的杨氏模量， d 是损伤因子，点是加载过程中，单元的最大应变， 毹是单元完全破坏时的应变， a 0 是达到损伤准则时的张开应变 2 5 本章小结 本章主要是阐述了水力压裂有限元数值模拟涉及到的力学原理和实现方法， 它是研究水力压裂规律的有力工具，并简单介绍了水力压裂数值模拟涉及到的多 们力学学科。 第二节阐述了岩石孔隙介质模型。 第三节介绍了渗流应力耦合模型的数学模型和算法。 第四节介绍了c o h e s i v e 单元的损伤力学模型。 中国科学技术大学硕士学位论文 第三章水力压裂的三维数值模拟 3 1 垂直裂纹扩张的瞬态模拟 3 1 1 问题描述 本数值模拟是采用c o h e s i v e 单元预设裂纹扩张的路径再加上a b a q u s 的应力 渗流耦合模型来模拟水力压力导致的垂直裂纹在三维空间中在长宽高三个方向 的扩张。在图l 的几何模型中：轴2 为高度方向即重力方向，是垂直裂纹的高度 扩张方向；轴l 、3 为水平方向，分别是裂纹扩张的长度和宽度扩张方向。由于 裂缝扩张在水平和高度扩张方向( 轴l ，轴2 ) 多层，所以本模型取垂直裂纹几 何模型的四分之一。模型的基本尺寸为l 方向长1 0 0 m ，3 方向宽8 0 m ，2 方向高 4 0 m 。考虑到油层的高度一般在1 5 m 左右，模型在高度方向上分为两层，在图l 中模型上层用橙色表示的单元就是岩层单元，下层用绿色表示的就是油层单元。 其中c o h e s i v e 单元预设的裂纹扩张平面为几何模型中两条红线所在的垂直于轴 3 的平面，初始裂纹位于平面内，用蓝色的线在图上标出。 图3 1 几何模型以及有限元网格划分 中国科学技术大学硕士学位论文 本模型中裂纹扩张的方向在水平方向上不受限制，但在高度方向上有岩层的 限制，所以需要