岩体单裂隙流固热化学耦合作用数值模拟研究

第４６卷第９期　２　０　１　５年５月　人　民　长　江　

Ｙａｎｇｔｚｅ　Ｒｉｖｅｒ　Ｖｏ１．４６．ＮＯ．９　Ｍａｙ，　２０１５　

文章编号：１００１—４１７９（２０１５）０９—００８４—０７　

岩体单裂隙流固热化学耦合作用数值模拟研究　

田晓丹　，姜晓桢２　（１．河海大学文天学院水利工程系，安徽马鞍山２４３０００；　２．南京水利科学研究院岩土工程研究所，江苏南京２１００２４）　摘要：为全面研究渗流、应力、化学、温度影响下岩体裂隙的变化，建立了岩体单裂隙的流一固一热一化四场全　耦合数学模型，并以Ｃｏｍｓｏｌ　Ｍｕｈｉｐｈｙｓｉｃｓ软件为基础，将该数学模型转化为一个统一的偏微分方程组，实现了　四场全耦合数值求解。该模型给出了更接近真实物理过程的数值解答，避免用松散耦合法求解多场耦合问题　带来的误差。同时通过对两个实例的计算模拟，分析了在渗流、压应力、水化学溶蚀和温度等多因素耦合作用　下，岩体单裂隙渗透特性的变化规律，得出了裂隙岩体渗透特性的影响机理和关键影响因素。　关键词：岩体裂隙；渗透特性；多场耦合；数值模拟　中图法分类号：Ｐ６４２　文献标志码：Ａ　ＤＯＩ：１０．１６２３２／ｊ．ｃｎｋｉ．１００１—４１７９．２０１５．０９．０２０　

裂隙岩体的渗流及其破坏问题是岩体水力学研究　的热点和难点问题，涉及水利水电工程、隧道工程、采　矿工程等众多领域。裂隙岩体在多种因素作用下所产　生的物理、化学的或者力学的作用过程正是导致岩体　发生破坏的重要原因之一。由于裂隙岩体所处环境的　复杂性，使得数学模拟占有非常重要的地位，建立一个　能够反映渗流一应力一化学一温度作用及孔隙结构和　裂隙面几何形状变化之间相互影响的数学模型，将对　工程建设起到重要指导作用。　目前很多研究都采用多组份化学反应传输模　型¨Ｉ２　或者它们的简化模型　来描述岩石裂隙面矿　物反应的时空变化过程。Ｙｅｈ　Ｇ．Ｔ．，Ｔｒｉｐａｔｈｉ　Ｖ．Ｓ建立　了溶解和沉淀过程中裂隙渗透率与变化中的孔隙率之　间的关系，但很难反映孔隙空间的溶解、沉淀、闭塞、次　裂隙发育等作用引起的裂隙面几何形状的影响。裂隙　覆盖物的沉淀及表面形态的变化将改变裂隙壁面之间　和岩石内的有效弥散率，改变岩石组分和中和势　。　Ｄｒｅｙｂｒｏｄｔ　Ｗ等对石灰岩的溶解过程做了深入的研　究　Ｊ，并建立了石灰岩的表面溶解、离子的对流扩散、　ＣＯ，的转换和水流流动的耦合求解模型，模拟了石灰　岩在一转盘中的溶解过程。该模型较为全面地考虑了　石灰岩的溶解过程，但是该模型只考虑了石灰岩的表　面溶解，并没有考虑多场耦合条件下石灰岩的溶解过　程。Ｙａｓｕｈａｒａ和Ｅｌｓｗｏｒｔｈ等人利用规则的球型矿物颗　粒模型从压力溶解的角度对天然裂隙渗透特性进行了　理论模型研究　，但其理论是基于岩石颗粒之间的　压力溶解理论建立起来的，文中的有些参数取值比理　论值大上千倍，是否适合裂隙还有待进一步的验证。　Ｌｉｕ和Ｓｈｅｎｇ等人利用理想化锯齿裂隙面模型，从裂隙　面化学溶蚀和接触面压力溶解相互竞争的角度建立了　二维裂隙渗透特性应力一渗流一化学耦合模型　，对　Ｐｏｌａｋ的试验进行了理论和数值解析”…，但并未对温　度场的作用加以考虑。由以上综述可见，目前国内外　的研究工作主要集中在耦合作用下裂隙岩体变形力学　特性机理、力学本构关系以及耦合模型建立上，而关于　裂隙岩石的流固热化学耦合模型建立的研究还较为薄　弱，缺乏系统研究渗流一应力一化学一温度耦合作用　下裂隙岩石的渗透特性变化规律，尤其是从细观内部　结构模型出发来考虑矿物化学溶解、溶质迁移和沉淀　过程中渗透特性的变化规律。　

１流固热化耦合模型　本文在前人工作的基础上，给出了岩体单裂隙流　固热化四场耦合的数学模型，其基本假定如下：①岩　体为饱和、非沸腾多孔弹性介质，孔隙被水充满，并假　设ＲＥＶ存在；②在应力场中，只考虑压力溶解产生的　

收稿日期：２０１４—０８—１５　基金项目：国家自然科学基金资助项目（５０７７９０１２）　作者简介：田晓丹，女，讲师，硕士，主要从事水力学及河流动力学研究。Ｅ—ｍａｉｌ：ｔｉａｎｘｉａｏｄａｎｈｈｕ＠１６３．ｅｏｍ　第９期　田晓丹，等：岩体单裂隙流固热化学耦合作用数值模拟研究　８５　影响；③岩体中水的流动符合达西定律；④不考虑相　变，即不考虑因温度引起的气相与液相之间的转化；⑤　热量变化包括热传导和热对流，热传导过程符合Ｆｏｕ—　ｒｉｅｒ定律；⑥水的密度不随时间变化，只随温度变化；　⑦溶质运移过程中考虑吸附作用，并假定为等温线性　吸附；⑧不考虑岩体受热传导影响产生的温度应力；　⑨温度场只单向地影响溶质在流体中的扩散速率以　及溶质的溶解速度。　１．１渗流方程　在多孔介质中，渗流符合达西定律。考虑裂隙开　度变化对渗流场的影响，需在达西公式等号的左侧添　加裂隙开度随时间变化的量，如式（１）所示　旦（ｋｂ　）＋旦（ｋｂ　）：　（１）　８ｘ　ｄｘ　８　ｄ，ｄｔ　式中，ｋ为渗透系数，ｂ为裂隙开度，　为渗透水头。　１．２压力溶解控制方程　如图１所示，参考文献［７］建立了一个力学模型　来描述岩石颗粒的压力溶解过程，该过程导致了渗透　特性的改变和压力的增加。　图１颗粒压力溶解示惫　在应力作用下，溶解将在颗粒与颗粒接触部分的　薄膜水中发展，物质将通过溶解被输送。物质的溶解　将随着接触面积的增大而增大，使接触应力减小，这个　过程将一直持续到接触应力减小到一定程度才停止。　这个临界应力可以通过一定应力和温度条件下的能量　平衡得到，可表示为　Ｅ　（１一÷）　ｃ　—■　一　Ｌ　式中，ｏｒ　为临界应力；　为熔点；　为固体的摩尔体　积；Ｅ　为溶解热；Ｔ是系统温度。　在颗粒接触边界的溶解应变率可以表示为　Ｏ￣＇ｄｉｓｓ＝　（一　）（３）　Ｊ　Ｌｊ　式中，　为溶解应变率；ｋ＋为溶解率常数；Ｒ为气体　常数；ｄ为颗粒直径；　为受压力作用的颗粒部分和　不受压力作用的空隙部分的化学势差，这是压力溶解　的驱动力。化学势可以表示为　△　＝（ｏｒ。＋Ｐ）Ｖ　＋Ａｆ一２ＨＴ　（４）　式中，ｏｒ。为分离压力，等于接触部分的压力超过空隙　水压力的数值；Ｐ为空隙中的压力；ｆ为摩尔赫尔霍兹　自由能，它与弹性和塑性应变能有关；式中第３项代表　了颗粒与颗粒接触部分的表面能，日是流固接触面的　平均曲率，　是颗粒与颗粒接触部分的具体界面能。　由于孔隙中的压力为孔隙水压力，可以用ａｆ表示　为颗粒与颗粒接触部分与空隙自由表面的摩尔赫尔姆　霍兹自由能之差。当整个系统达到平衡时，化学势差　将变为零，接触部分的应力Ｏｒ　将变为临界应力ｏｒ　，因　此可以得到　Ａ　ｑ　：　＝ｏｒ　＝２ＴＨ　一气　（５）　ｙ

ｍ　式中，上标ｅｑ表示达到饱和状态下相应的量。把式　

（５）代人式（４）中，同时由于空隙中的压力Ｐ较小可忽　略，可以得到　

一△　＝（ｏｒ。一ｏｒ　）Ｖ　（６）　将式（６）带人式（３）中，可以得到颗粒与颗粒接触　部分的溶解应变为　３　ｋ＋，　、　…　ｓ　面Ｌ　ｏｒ—ｏｒｃ，　Ｌ　

根据图１所描述的力学模型可知，在压力ｏｒ。作用　下，颗粒之间接触部位发生溶解，导致两个颗粒之间的　距离变近，此时这个变形可理解为该裂隙面在此处开　度的变化量，即随着裂隙面接触部位的溶解，相互接触　的两个裂隙面的距离变近了，所以溶解应变率在此处　的含义就变为，单位时间内两个相互接触的裂隙面在　垂直裂隙面方向上的应变率。裂隙面在初始时刻的开　度等于两个颗粒中心之间的距离，即颗粒的直径ｄ；经　过一段时间，颗粒接触部分因压力而发生溶解，导致其　在垂直裂隙面上的应变率为溶解应变率　。，即应变　率８　乘以初始开度ｄ（颗粒直径）可得压力作用下裂　隙岩体的溶解量（即裂隙开度的改变率），即式（７）两　边同乘以颗粒直径ｄ可得式（８）　ｄｓ１　３　ｋ＋（ｏｒ。一　）　（８）　

式中，ｄ　为裂隙开度，其与单位时间ｄ　的比值即为裂　隙开度改变率。　虽然颗粒接触和裂隙面接触的不同在于接触部分　的几何模型不一样，但可以将裂隙面接触的几何模型　概化为颗粒接触的几何模型，如图２所示。　以接触直径ｄ　代替颗粒直径ｄ可得　

ｄ　＝　／　（９）　Ｖ　丌　

式中，Ａ　为宏观接触面积。　对于单轴压缩，法向压力作用在裂隙的接触部分　８６　人　民　长　江　上，可得　Ａ‘　。Ａ　Ａ　８　

（１０）　

（１１）　式中，　为平均的宏观应力；Ａ　为总的接触面积；　为接触比，即接触面积与总裂隙面积的比。　以上步骤得出了在平均宏观应力下，裂隙接触部　位压力对裂隙岩体溶解作用的控制方程。　裂隙的接触部分　局部接触面积　图２裂隙接触概化　１．３溶质运移方程　溶质在多孔介质中运移有两个主要过程，即对流　和扩散过程，但在迁移中还有吸附和解吸附（溶解）过　程。本文中的模型除了考虑吸附、对流、扩散和溶解４　种基本作用外，还必须考虑１．２节中的压力作用下，裂　隙接触部位产生的溶质，所以综合上述５种作用得出　如式（１２）的控制方程　３，　１　Ｒｄ　＝　（Ｄ　Ｃ）一ｕＶ　Ｃ＋÷［Ｄｄ（Ｃ　一Ｃ）＋　Ｏ＇ｂ　ｕ　Ｐ　３　ｋ＋＋（　。一　。）］　，，＾、　———　—一Ｊ　式（１２）中等号左侧为溶质浓度随时间的变化量，　等号右侧从左至右分别代表了溶质的扩散作用、对流　作用、吸附解吸以及压力溶解５个过程。其中，Ｃ为裂　隙中流体的溶质浓度，由于裂隙流体中一般都有一个　初始的溶质浓度（一般为岩体中溶质的浓度，即溶质　的饱和浓度），所以这里的溶质浓度ｃ从更准确的意　义上可理解为相较于初始溶质浓度的变化量；Ｄ为扩　散系数；／．ｔ为裂隙中流体运动速度；Ｄ　＝Ｄ／ｇ，　为裂　隙岩体介于溶解与未溶解之间的扩散边界层厚度；Ｃ　为岩体中溶质的浓度；Ｐ　为固体溶质的密度；Ｍ为溶　质摩尔质量；ｋ＋为溶解率常数；Ｒ为气体常数；Ｔ为系　统温度；　为临界应力，即物质的溶解伴随着接触面　积的增大，导致了接触应力的减小，当溶解停止时的接　触应力；　为分离压力，等于接触部分的压力超过空　隙水压力的数值。　１．４热动平衡方程　地下水在岩体裂隙中的热量运移主要有３种形　式，即传导、对流和辐射，其中含水介质中骨架的温度　和地下水的温度并不相同，也存在热量的交换，可以把　含水介质中的水温和骨架温度作为两个未知函数，建　立如式（１３）的热动平衡方程。　Ｔ　ｃ　＝　（ｎＡ　Ｖ　Ｔ　一ｃ　）＋Ｋ　（　—　）　（，　（１３）　式（１３）中等号左侧为系统热量随时间变化的改　变量，等号右侧第一项为热量的传导、对流和辐射项，　第二项为含水介质中骨架与水之间的热交换过程。其　中，ｃ　为水的热容量；Ｔ为水温；Ｔ为骨架温度；Ａ　为热动力弥散系数；／Ｉ，为岩体孔隙率；Ｋ　为热交换系　数，　和含水介质的性质有关，当　＞　时，热量从　骨架进入水中，Ｋ　取负值，当Ｔ＜Ｔｗ时，热量从水进　入骨架，故Ｋ　取正值。