含裂隙多孔介质非饱和渗流初步分析

《含裂隙多孔介质非饱和渗流初步分析》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

三峡大学硕士学位论文1绪论1.1问题的提出长期的物理和化学风化、地质构造以及卸荷等作用，使得天然岩石中存在着大量的裂隙；而气候的变化、干湿循环等作用也使得某些类型的土体形成裂隙。裂隙岩土体中地下水位以上部分是未被水充满的非饱和带，降雨入渗和地面水体的下渗都是通过该带到达稳定地下水面的非饱和渗流过程．可以说，在自然界中裂隙岩土体非饱和[1]渗流是普遍、客观存在的。虽然岩石和土体中的裂隙成因不尽相同，但水在裂隙及其赋存的岩土介质中渗流特性有相似性，大量试验结果也表明水在裂隙和岩土介质中[2][3]的饱和、非饱和都满足达西定律。所以，可以将裂隙岩土体渗流统一起来研究，不同之处主要表现在介质渗透性质的不同。工程建设中，对于岩基上的大坝、地下洞室、岩石边坡等工程，水在裂隙岩土体中的流动而形成的荷载对工程的安全是至关重要的。如1976年，美国的Teton土坝溃坝，经反复查证，确认坝基岩石节理发育，库水流经岩石裂隙使心墙齿槽土体发生管涌而最终遭致溃坝。还比如，我国的梅山连拱坝，1958年开始蓄水，1962年9月到11月在高水位持续了40天，右岸坝座错动，岩石裂隙大量渗水，坝体出现不正常[4]位移，支墩地基测到很大的扬压力后，紧急放空水库，才避免了溃坝灾害。另外，各类工程建设也常被某种工程性质极差的特殊土——裂土（裂隙粘土、膨[5-9]胀土、胀缩土的统称）所困扰，造成巨大经济损失。这类土危害工程建设的一个重要表现形式就是裂土边坡失稳。因为裂土具有吸水膨胀、脱水收缩的特性，裂土边坡在蒸发条件下含水率降低，从而导致土体形成许多不规则的裂隙，裂隙的存在破坏了土体的完整性，同时为水分的入渗开辟了通道。在降雨条件下，土体裂隙的存在使得水分渗入土体更加容易，从而引起土体渗流场、应力场的变化，极易导致边坡失稳[11]。而在裂土地区修筑的铁路、公路路堤与堑坡及渠坡也常在雨季产生溜坍、坍塌，甚至滑坡，引起边坡破坏，妨碍工程设施的正常使用和运行，且这种破坏往往具有多[6-10]发性、反复性及长期潜在性。近年来已经充分认识到土体裂隙对土体渗流和边坡[11]稳定的重要影响。总之，随着社会生产的不断发展，深入开展裂隙岩土体非饱和渗流问题研究不仅具有突出的理论意义，而且具有重大工程实用价值。1.2研究现状流体通过多孔介质或裂隙介质的流动，即渗流，是一种与人类的一些工程活动密切相关的现象，对于渗流的研究最早开始于土体中。1856年，法国工程师Darcy总1 三峡大学硕士学位论文结出线性渗流定律，即著名的达西定律，标志着经典渗流理论的诞生。100多年来，建立在达西定律基础上的经典渗流理论在工程实践中得到了广泛的成功的应用。经典渗流理论是以连续介质假定为基础的。而近代连续介质力学是建立在质点或“表征体积单元”的基础上的。引入了表征体积单元，就可以把实际上为离散型的多孔介质抽象为由表征体积单元组成的连续介质，只要表征体积单元与研究域的尺寸相比很小，就可以用连续介质渗流理论进行分析。对于土体，表征体积单元很小，连续介质渗流理论可以直接进行应用。岩体的渗流有着不同的特点。岩体中一般分布有大小不等的裂隙。在流体通过岩石流动的过程中，尽管孔隙也可以作为渗流通道，但裂隙扮演着更重要的角色。裂隙和孔隙的关键性的差别有两点。首先，孔隙在三维方向上延伸，其各个方向的尺寸的差别不显著，而裂隙在二维方向上延伸，其两个方向上的尺寸比第三个方向上的大得多。其次，孔隙的延伸尺度比研究域小得多。例如孔隙的尺寸一般为几毫米或更小，而研究域一般大于几十米。然而，裂隙的延伸尺度可以达到整个研究域，其大小可以从几厘米，或更小，变化到几千米。这些差别则导致裂隙岩土体中的渗流规律不同于多孔介质。因为岩石与土体中裂隙存在着差异，所以在岩石裂隙渗流方面，国外有大批学者相继开始了岩石裂隙的试验和理论研究。WangJ.S.Y.于1991年对1987到1990年期间内有代表性的论文进行了评述，发表了具有代表性的综述性论文“Flowandtransport[12]infracturedrocks”，此时西方国家在该领域的研究达到了高潮。在国内，田开铭、张有天、张武功等人最早开展了对岩石渗流的研究，给出了关[13][14-17]于裂隙岩石渗透性的初步研究，揭示水在裂隙中的运动状态，以及裂隙岩体渗[18][19]流特性、数学模型及系数量测等，拉开了国内裂隙岩体渗流问题研究的序幕。目前，裂隙土体非饱和区的渗流问题还存在一些值得关注的现象。[20]首先，在裂隙土体非饱和区，水流常常以毛细管流、片流和间歇流的形式进行。毛细管流假定流体与孔隙或者裂隙两侧壁均接触，从而产生表面张力，形成负孔隙水压力，流体渗流由孔隙水压力和重力梯度之和，即水力梯度来驱动。它常见于多孔介质渗流。片流与毛细管流不同，它只与裂隙一面侧壁接触，裂隙两侧的片流中间由空气隔开，互不连续。片流主要受重力梯度驱动，几乎不受裂隙毛细管力影响。后有学[21]者在实验中观测到裂隙岩体中间歇流，他们发现小股流动通过裂隙岩体时，流动先是在某一个孔隙附近被吸附，然后流体在吸附点前面累积，当毛细吸着水累积到一定程度时，就接着向下移动，又形成新的小股流动，如此反复向前推进。这种流动形式目前还没有很好的理论来解释。2 三峡大学硕士学位论文其次，土体渗流时存在这样一种现象，即在非均匀流中水分和溶质可通过少部分[20]土体快速运移，通常定义为优势流（Preferentialflow）。国外已有大批学者针对土壤中的优势流问题展开了初步探索，按其形成原因将优势流分为大孔隙流(macroporeflow)、指流(fingerflow)、绕流(bypassing)、漏斗流(funnelflow)、沟槽流(channeling)、短路流(shortcircuiting)、分驱动流(partialdisplacement)和地下风暴流(subsurfacestormflow)等，并认为影响土壤中优势流的主要因素有：土壤结构、土壤质地、土壤水分含量、土壤初始水分含量、水与溶质的施加速率及溶质的施加方法。目前，优势流的[22][23]研究主要限于农业土壤、地下水污染等领域，直接针对岩土工程中裂隙岩土体[24]渗流问题的研究还不多见，值得去探索。另外，水流在饱和裂隙中流动时，有一部分水被岩体基质所吸入，这种吸入作用的大小对水流运动的影响及其相互作用[20]机理也尚未研究清楚。上述这些问题都表明裂隙岩土体非饱和流动性质的复杂性。为了研究裂隙岩土体的渗流特性，国内外的学者纷纷建立了或基于等效连续介质理论，或基于离散裂隙网络的适用于各种工程条件下的数值模[25]型。1996年，Altman等提出了一套裂隙岩土体非饱和渗流概念模型（如图1）。其中图1a是等效连续多孔介质模型，对于饱和裂隙渗流，需要确定等效渗透张量；对于非饱和裂隙渗流，需要给出土水特征曲线（吸力或孔隙水压力，或毛细管压力与含水率之间的关系）和渗透性函数（吸力或孔隙水压力，或毛细管压力与渗透系数之间的关系）。图1b是复合孔隙等效连续介质模型，它假定在含水量低的时候，流体只在岩土介质中流动，在高含水量的时候，流体同时在岩土介质和裂隙中流动，图1裂隙岩土体渗流概念模型相应渗透性函数为岩土体和裂隙渗透性函数的复合。图1c为双重孔隙模型，它假定裂隙和岩土体是相互作用的等效连续介质，3 三峡大学硕士学位论文它们之间存在流体交换，流体只在渗透性高的裂隙中流动。该模型常用于饱和裂隙岩土体渗流分析，而很少用于非饱和裂隙渗流。图1d为双重渗透性模型，它是双重孔隙模型的推广，它假定岩土体和裂隙中均存在流体流动，两者流体交换量与两种介质中水头差成正比。该模型常用于模拟裂隙岩土体中的非饱和优势流，它需要给出岩土体和裂隙各自的土水特征曲线和渗透性函数。图1e和图1f是离散裂隙网络模型，其中图1e不考虑岩土体渗流，而图1f则同时考虑岩土体渗流。图1e主要用于分析饱和渗流，它需要给出各个裂隙的几何形状和渗透性。随后它也被推广到非饱和裂隙渗流中，如图1f，此时它需要给出各个裂隙的几何形状，岩土体和各个裂隙渗透[20]性函数和土水特征曲线。Altman等提出的概念模型得到许多研究者认同，他在等效连续介质模型和离散介质模型之间建立了分级子模型，每个子模型描述了一类裂隙分布特性岩土体渗流特征。[26]目前国内描述裂隙岩土体渗透特性的数值模型主要归纳为三类。第一类是等效连续介质模型。原理是将裂隙岩土体概化为等效连续介质，其渗透[27][28]性用等效各向异性渗透张量来描述。它包括了图1a、b、c、d四种模型。它从宏观上将离散随机裂隙渗透性等效为连续介质渗透张量，利用成熟的连续介质理论分析问题；但它难于确定裂隙岩体等效渗透张量，以及等效连续介质模型有效性和适用[26]性有限。应用等效连续介质模型求解裂隙岩体非饱和渗流场的主要代表有：[30][31][32][33]Dykhuizen；Peters和Klavetter；Pruess,Narasimhan&Tsang以及Finsterle等。第二类是离散裂隙网络模型。即建立裂隙网络，然后研究流体在裂隙网络与岩土体介质中的渗流。对于饱和裂隙渗流，通过裂隙的流量与等效裂隙水力隙宽的立方成[34][35]正比，即渗流遵循所谓立方定律。而对于非饱和裂隙渗流，裂隙渗透系数是孔[20][36]隙水压力或者基质吸力的函数。这类模型也即图1e所示的模型。国内从事饱和[37][38][1][39]-非饱和裂隙网络渗流研究的学者有张有天、王恩志、胡云进和宋晓晨等。第三类是裂隙孔隙介质模型。也即图1f所示的模型。它能同时反映裂隙及孔隙的水力特性，如孔隙压力在裂隙中传播速度很快，而在孔隙介质中传播速度则很慢的特性。这是一种较为理想的模型，可弥补以上两种模型的不足。该模型的建立与算例分析将是本课题的研究重点之一。1.3本文的主要研究内容（1）裂隙和孔隙介质的非饱和土水特征曲线试验和渗流性函数预测。（2）介绍岩土体中裂隙的几何参数、测量方法；裂隙面产状和尺寸的统计方法；裂隙的密度和频率等裂隙的基本特征等。基于拉丁方抽样（LHS）方法，实现随机裂4 三峡大学硕士学位论文隙网络的生成。（3）等效连续介质模型的控制方程、求解渗流场的有限单元法以及边界条件。渗透张量的性质。分析该将模型用于裂隙岩土体所必要的前提条件，以及渗透张量的求解。编程实现等效连续介质模型的建立，并通过算例来模拟含裂隙岩土体的基本渗流情况。（4）裂隙孔隙介质模型及其建立的过程、适用条件以及渗流的相关特性等。通过算例来展示该模型的特征，并讨论渗流参数的改变对裂隙岩土体渗流的影响以及不同边界条件的改变对裂隙岩土体渗流的影响。5 三峡大学硕士学位论文2孔隙裂隙介质的渗流特性2.1孔隙介质的渗流特性土体作为一多孔介质是由无数碎散的、形状不规则且排列错综复杂的固体颗粒组成的。在对土体进行研究时，通常用宏观的方法，也就是在较大尺度范围内研究多孔介质的孔隙大小及其中水流的平均状况。如图2.1所示，取一表征单元体（REV）对土体进行研究。表征单元体的选取很关键。体积V1取得过小，那么它空隙率离散性太大，逐渐增大体积至V2，孔隙率有波动但逐渐减小，继续增大体积，直到出现某一体积，使得继续增大体积对空隙率）%（孔隙率体积图2.1非饱和土典型单元体(REV)的改变不产生较大影响时，也就是说该体积的孔隙特性基本能代表整个土体的孔隙时，该体积则可认为是土体的表征单元体。2.1.1非饱和土中水的存在形式与特征水气相在土单元中的存在形式多样。当土单元含水量特别小时，空气基本填充了整个单元孔隙，气相在孔隙中是完全连通的，如图2.2（a）；当含水量增加，空气体积逐渐减小，成为部分连通或内部单元的局部连通，如图2.2（b，c）；当土趋于饱和时，，气泡则被封闭在孔隙中或溶解在水里，如图2.2（d）。6 三峡大学硕士学位论文图2.2水气相存在形式示意图（包承刚，1979）天然土中水气相的关系也可以表示为图2.3所示关系图：水相和气相的性质饱和土非饱和土孔隙被水封闭气泡干土充满气相不连续无自由水存在气泡溶解在水中水相连续水相不连续气相连续水本质上不水、气各自流动移动图2.3天然土中水气相的关系从宏观上看，多孔介质中的空隙主要包括基质孔隙与裂隙，从而，水在其中的存在形式有以下几种形式：a）结合水；因为岩石颗粒表面及岩石空隙表面带有电荷，具有吸附水的能力，当固体表面的引力大于水分子自身重力，会在颗粒表面形成一部分水，即结合水。其不能在自身重力影响下运动。b）重力水；重力水是指能在自身重力影响下运动的这部分水，它距离颗粒表面较结合水远，所以重力对它的影响大于固体表面对它的引力。7 三峡大学硕士学位论文c）毛细水；水气界面下，因表面张力作用的存在，水会承受一定的吸持力，承受这种毛细管力的水被认为是毛细水。一般而言，毛管的直径在0.5~0.0002mm之间，裂缝宽度在0.25~0.0001mm之间的孔隙称为毛细管孔隙。多孔介质中的细小空隙通道以及裂隙会构成毛细管。在毛管孔隙中，液体质点之间以及液体与孔隙壁之间均处于分子引力的作用下，因毛细管力的作用，液体不能自由流动。若要使液体流动，需要施加外力去克服毛细管阻力。毛细管孔隙越小，产生的毛细管力越大。2.1.2孔隙介质的土水特征曲线2.1.2.1土水特征曲线定义土水特征曲线定义为土的含水量与土的吸力之间的关系曲线。土的含水量可以是重力含水量（ω），体积含水量（θ），也可以是饱和度（S）。吸力可以是基质吸力（也）%体积含水量（基质吸力（kPa）图2.4典型土水特征曲线（Fredlund和Xing，1994）叫做毛细压力（uu−），u是孔隙气压力，u是孔隙水压力），也可以是总的吸力awaa（基质吸力与渗透吸力的总和）。当吸力值很高时（一般大于1500kPa)，基质吸力和总吸力一般可假定为相等的。土水特征曲线经常绘制在半对数坐标系上（如图2.4）[42]在描述非饱和土性质的方程中还用到另外两种形式的含水量，我们称其为归一化重力含水量或归一化体积含水量，它们通常用残余含水量来表达：ww−rΘ=（2.1），nww−sr8 三峡大学硕士学位论文θ−θrΘ=（2.2）θ−θsr其中w,θ——重力含水量，体积含水量；w,θ——残余重力含水量，残余体积含水量；rrw,θ——饱和重力含水量，饱和体积含水量。ssw无量纲含水量定义为：Θ=（2.3）dws在土体不发生变形的前提下，从含水量的各种不同描述中都等得到相同的信息，如进气值和残余含水量等。如图2.4所示，曲线上有两个特征点：一个对应于土的进气值uu−的点，土的aw进气值是指空气开始进入土体边界的土颗粒或颗粒集合体之间的最大孔隙时对应的基质吸力值；另一个特征点是对应于残余含水量θ的点。当土体中含水量随着吸力的r增加而降低到一定值时，含水量的继续减少需要增加很大的吸力，含水量的这一临界值称为残余含水量。一般通过作图可获得残余含水量的值。如图2.4所示，过曲线的拐点作一条切线，高吸力范围内的曲线也可以表示为另一条直线。两条直线的交点所对应的纵坐标即为残余含水量。当土的含水量为零时，所有土类的总吸力大致相等，6其值略小于10kPa。同一土体，对应于吸湿过程和脱湿过程有两个不同的土水特征曲线，如图2.4所示。这种现象称为滞后现象。产生滞后现象的原因主要是由于土中孔隙的几何形状所引起的。由于吸湿和脱湿两个过程存在的瓶颈现象而导致相同的吸力值对应于不同的含水量。另外，脱湿过程和吸湿过程中的液、固表面接触角不同也是引起滞后现象的原因。2.1.2.2土水特征曲线实验方法土水特征曲线描述的是吸力与体积含水量、重力含水量或者饱和度之间的关系，也被称为持水曲线，测量土水特征曲线的方法有很多，这里简要介绍四种：Hangingcolumn法、压力板仪与体积控制或质量控制法、冷静式湿度计法以及离心机法。其中Hangingcolumn法与离心机法主要在低吸力范围内量测，以得到进气值。而压力板仪与体积控制（或质量控制）注重土水特征曲线中间段的测量1、Hangingcolumn法Hangingcolumn法，适用的吸力范围是0-80kpa，一般用于细粒含量少的粗粒土的脱湿过程。装置如图2.5所示，将土样连同环刀放置在多孔陶土板上，为保证土样与板接触良好，将其旋转约45°，并在土样上表面放置一块刚性板，板上放置1kg重物。然后9 三峡大学硕士学位论文在漏斗表面放一块带小孔的塑料薄片以防止蒸发且能和大气相通。土样达到平衡至少48小时，直到水气界面不再改变，记下平衡时毛细管中的水气界面位置。在一次平衡之后，调整水槽的高度来改变吸力，其大小可通过压力计来观测。注意第一次的压力增量不能超过预期的进气值。监测水气界面的位置，直到其24小时内都不发生变化，记下该位置。重复上述过程，实施下一级吸力。图2.5Hangingcolumn法测量装置2、压力板仪与体积控制（或质量控制）图2.6压力板仪压力板仪与体积控制（或质量控制），适用的吸力范围是0-1500kpa，一般用于细粒土。10 三峡大学硕士学位论文装置如图2.6所示，体积控制法中，将饱和土样与饱和高进气值陶土板充分接触，先将陶土板下面的空气排净，并充满水和量管连接，以保证试样水量变化可友量管指示读出，容器内从零逐级加压，待吸力平衡时，读出量管指示数。逐级加压并记录相应含水量变化，以便反算含水量。施加最后一级空气后，记录吸力值，并将试样称量，烘干后计算重力含水量。再利用两次量管体积读书差反算不同吸力值下试样含水量，绘出土水特征曲线。质量控制法中，将饱和试样与高进气值陶土板充分接触，陶土板下的底板上有一排水管，供土样图2.7冷镜测试仪排水之用。通过施加气压控制吸力，待吸力平衡后称量试件和仪器重，以便测定含水量变化。施加更高气压后，重复这一步骤。施加最后一级气压后，取出试样，烘干测定含水量，并根据前面测定的试样和仪器重，反算不同吸力值下对应的重力含水量值，绘制干燥段土水特征曲线。3、冷静式湿度计法（Chilledmirrorhygrometer法）冷静式湿度计法，适用于500kpa-100Mpa，一般用于高吸力范围，侧重得到干燥段的结束时刻的值。装置如图2.7所示，将土样放进托盘，并插入到冷镜湿度计内。通过该湿度计测得相对湿度，并精确到0.001。对不同体积含水量的试样进行相对湿度的测试，然后将试样称重，并精确到0.001g，得到质量含水量。4、离心机法离心机法，一般用于0-120kpa，仅对粗粒土实验。装置如图2.8所示，将土样置于封闭容器内，记下旋转轴中心到试样底部、顶部的距离（外、内旋转半径），精确到0.1mm。将该容器放入离心机内，以某个吸力相对应的离心速度对土样离心120分钟，记下在该过程中土样所失去的水分，精确到0.1ml。直到不再有水出来，记下失去水分的总量，精确到0.1ml。改变速度，图2.8离心机重复这一过程，即可得到土水特征曲线。11 三峡大学硕士学位论文2.1.2.3土水特征曲线估算1、土水特征曲线估算的理论基础目前，除了通过实验测定，还可以根据土颗粒和孔隙大小分布来估算土水特征曲线。Fredlund和Xing在1994年提出了根据土体空隙孔径分布曲线得到土水特征曲线[43]的理论，用方程（2.5）来近似模拟土水特征曲线：mΘ=ae−()bsψ（2.5）s式中，a、b、m为拟合参数。ss但此方程并非普遍适用，它受到吸力范围的限制。为建立土水特征曲线的估算理论，先考虑土中的孔隙大小分布曲线。土体被认为是一种内部孔隙随机分布的孔隙介[44]质。用孔隙半径r和半径分布函数f(r)来描述裂隙特征。这样，半径由r到((r+dr)的孔隙中充满水时的含量为f(r)dr，水的饱和度可表示为：Rθ()R=∫frdr()（2.6）Rmin含水率表示为：Rθθ()()R=Rf∫()rdr（2.7）Rmin式中：θ()R——半径小于或等于R的空隙中充满水时的水含量；R——土中最小的孔隙半径。min若R表示最大孔隙的半径，则达到饱和时的情况为：maxθ()R=θ（2.8）maxs根据毛管理论，毛细吸力与空气—水界面收缩膜的曲率半径成反比关系。有：CT2cosϕr==（2.9）ψψ式((2.9)即为Kelvin毛细模型方程。式中：T——孔隙水的表面张力，ϕ——水土间的接触角。曲率较小（毛细压力较小）时，可用孔隙半径来近似收缩膜的曲率半径。两个极端吸力条件可表示为：Cψ=（2.10）maxRminCψ=（2.11）minRmax12 三峡大学硕士学位论文式中，ψ——与最小孔隙半径对应的吸力，minψ——与最大孔隙半径对应的吸力。max根据该毛管理论，方程（2.6）可以写为：ψψCCmaxCCθψ()==f()()df()dh（2.12）∫∫ψψ2maxhhhh式中：h——积分变量，代表吸力。方程（2.12）是描述含水量与吸力关系的一般形式。如果土的孔径分布f(r)已知，那么土水特征曲线由方程（2.12）唯一确定。2、不同孔径分布形式下的含水率—吸力关系方程为使土水特征曲线能在整个吸力范围内都适用，先将含水率写成饱和度的形式，即θ写成θ/θ，则方程（2.12）可写为：s∞θψ()=θf()hdh（2.13）s∫ψ式中，f(h)——孔径分布函数不同孔径分布形式下的含水率和吸力关系曲线如下：（1）服从均匀分布即：f()hA=(2.14)式中，A是一个不变量。则带入方程（2.12）得到：ψAC11Bθψ()==dhAC(−)=−D（2.15）∫ψ2maxhψψmaxψ式中：B=AC，是个定值，DA=C/ψ也是一个定值。max（2）服从幂指数分布m+即：f()hAr=/（2.16）式中：m为某一整数。代入方程（2.12），则含水率与吸力的关系为：m+1ψmaxAhCmθ()ψψ==dhBD−（2.17）∫ψChm+12mmm其中：B==Am(ψ)/(C),DA/(mC)，为定量。max（3）服从正态分布1−−()hμ22/2σ即：fh()=e(2.18)2πσ式中：μ为孔径大小的平均值，σ为标准差。代入方程（2.13），则含水率与吸13 三峡大学硕士学位论文力的关系表达为：∞θψ()=θf()hdhs∫ψ∞θs2−y2=∫edy（3.19）2π()ψμ−/2σθψμ−s=erfc()22σ（4）服从对数正态分布211⎡⎤⎛⎞lnh−μ即：fh()=−exp⎢⎥⎜⎟（2.20）2πσh⎢⎥⎣⎦2⎝⎠σ式中：μ为孔径大小的平均值，σ为标准差。代入方程（2.7），则含水率与吸力的关系表达为：2R11⎡⎤⎛⎞lnh−μθψ()=−θexp⎢⎥⎜⎟dhs∫Rmin2πσh⎢⎥⎣⎦2⎝⎠σ⎡⎤⎛⎞ln()C−μ⎛⎞ln(C)−μ（2.21）⎢⎥⎜⎟⎜⎟ψψ=−θφ⎢⎥⎜⎟φ⎜⎟maxs⎢⎥⎜⎟σσ⎜⎟⎢⎥⎜⎟⎜⎟⎣⎦⎝⎠⎝⎠式中：Φ为正态分布函数。当ψ达到无穷大时，则max⎛⎞Cln()−μ⎜⎟ψθψ()=θφ⎜⎟（2.22）s⎜⎟σ⎜⎟⎝⎠（5）服从gamma分布αβ−−1/h⎧he⎪,,αβ>≤≤0,0h∞αfh()=⎨βαΓ()即：（2.23）⎪⎩0,其他∞α−−1h式中：Γ（）=α∫hedh。0当α为整数时，代入方程（2.13）得到含水率与吸力的关系表达为：14 三峡大学硕士学位论文∞αβ−−1/hheθψ()=θdhs∫αβαΓ()ψ∞θsα−−1h=∫hedh（2.24）Γ()αψβ/α−1i−ψβ/ψe=θs∑it=0i!β（6）服从负指数分布−λh即：fh()=>λeh,0（2.25）代入方程（2.13）得到含水率与吸力的关系表达为：∞−λhθψ()=θλedhs∫ψC∞−λCψ=θλed（2.26）s∫ψψ−ψλ/=θes（7）服从β分布αβ−−11⎧hh(1−)⎪,,αβ>≤≤0,0h1即：fh()=⎨B(,)αβ（2.27）⎪⎩0,其他1α−1Γ()()αΓββ−1式中：Bh(,)αβ=−=∫(1)hdh0Γ+()αβ当α,β为整数时，代入方程（2.13）得到含水率与吸力的关系表达为：αβ−−11∞hh(1−)θψ()=θdhs∫ψB(,)αβ（2.28）αβ+−1⎛⎞αβ+−1iiαβ+−−1=−θψs(1∑⎜⎟(1−ψ))i=α⎝⎠i2.1.2.4土水特征曲线模型通过实验或者估算可以得到一系列表示吸力与含水量之间关系的离散的点。从而得到以吸力值为横坐标，饱和度为纵坐标的散点图。然运用这些数据来测量流量、压力和变形现象一般都要求是连续曲线的数学表达形式。所以，许多学者提出了一系列具有一定实用性的土水特征曲线模型，如表2.1。15 三峡大学硕士学位论文表2.1土水特征曲线经验模型编提出者方程参数说明号θ−θrΘ=,θ、θ分别为饱和含水率和1sr1Gardner(1958)Θ=θθsr−n1+qψ残余含水率，q、n为拟合参数λ⎛⎞ψ2Brooks和Corey(1964)Θ=⎜⎟bΘ同上，ψ为进气值，λ为拟合参数b⎝⎠ψm⎡⎤13vanGenuchten(1980)Θ=⎢⎥Θ同上，amn、、为拟合参数n⎣⎦1()+αψ4Williams(1983)lnψ=+ablnθab、为拟合参数5Mckee和Bumb(1984)−−()ψab/Θ=eΘ同上，ab、为拟合参数16Mckee和Bumb(1987)Θ=Θ同上，ab、为拟合参数()ψ−ab/1+em⎡1⎤7Fredlund和Xing(1994)θθ=s⎢n⎥θb为饱和含水率，amn、、为拟合参数⎢ln⎡ea+(/)ψ⎤⎥⎣⎣⎦⎦Feng和Fredlund（1999）abc+ψd8w()ψ=w体积含水量，abcd、、、为拟合参数dhysteresismodelb+ψ2.1.3孔隙介质的渗透性函数非饱和土中水相和气相的渗透系数均为孔隙比（e），饱和度（S）或体积含水量（θw）的函数，所以也常被称为渗透系数函数或渗透性函数。由于饱和度或体积含水量与基质吸力之间的关系可以由土-水特征曲线来描述，因此渗透性函数也常常表达为基质吸力的函数。常用的渗透系数模型有：a）Green-Corey模型Green-Corey模型预测非饱和渗透系数具有足够的精确度，公式为：2pmks30Tε−2kj()θ=⋅⋅⋅⎡⎤(212)+−ih，(1im=,2,…)（2.30）ij2∑⎣⎦kgnscρμwji=式中：k()θ——第i段的含水率（或吸力值）对应的渗透系数i16 三峡大学硕士学位论文kk/——匹配系数，其中k为实测的饱和渗透系数sscsk——计算得到的饱和渗透系数scT——水的表面张力ε——孔隙比p——孔隙特征经验参数，一般12≤p≤，可假定p=2ρg——水的容重wμ——水的粘滞系数n——nm=−[θ/(θθ)]，饱和含水率θ与最小含水率θ之间的计算间段总数，ssLsLi为间段编号，i=1表示第一个间段，对应与饱和含水率θ，im=表示最后一个间段，s对应于最小含水率θ；j为从i到j的计数；m为θ与θ之间等分的总数；h为吸力LsLj对应的水头高度。在拟合时，先计算饱和时的渗透系数（即吸力趋于零时），即：2pm30Tε−2kj=⋅⋅+⎡⎤(212)−ih，i=1（2.31）sj2∑⎣⎦ρμwgnj=1然后通过有限个数值点差值得到渗透系数随吸力（或饱和度）变化的关系曲线。b）Gardner模型Gardner模型是一个使用非常普遍的非饱和渗透系数的计算模型，其计算式为：ksk=（2.32）wβ⎛⎞uu−aw1+α⎜⎟γ⎝⎠w式中：k——随吸力变化的非饱和渗透系数wα、β——拟合参数，为不变量γ——水的容重wk——饱和渗透系数，可由渗透试验得到suu−——吸力awc）Brooks-Corey模型(2.52.0/)+λ⎛⎞θθ−rkk()θ=⋅⎜⎟（2.33）s⎝⎠θθ−sr式中：k——饱和时的渗透系数sλ——拟合参数θ——含水率，θ、θ分别表示饱和含水率和残余含水率。sr17 三峡大学硕士学位论文d）vanGenuchten模型2−mnn−1kasa⋅−{1([]u−uw)1(⎡⎤+[]aua−uw)}⎣⎦km==,(11/)−n（2.34）m/2n⎡⎤1(+−[]auuaw)⎣⎦式中：k——饱和时的渗透系数samn、、——拟合参数uu−——吸力aw2.2裂隙介质的渗流特性2.2.1裂隙岩体的渗流特点与相对成熟的多孔介质饱和渗流理论相比，非饱和带裂隙岩体渗流还处于研究阶段，一些概念模型相继被提出，但对其中一些现象所发生的机理还不完全清楚。如美国科学院等机构合作成立的一个非饱和带裂隙岩体渗流研究小组提出的一些重要的，[45]但在模拟上存在困难而需要研究的现象，包括毛细竹流、薄膜流、优先流和裂隙一基质相互作用等。a）毛细管流通常假定多孔介质中的空隙为一束束的毛细管。毛细管中的水受到表面张力的作用，一般都高于自由水面，需要施加外力才能使其流动。毛管中的水是处于负压状态的，负压的大小与孔隙大小和饱和度有关。孔隙裂隙介质主要由基质孔隙和裂隙组成。用此毛管理论来研究多孔介质毛细通道中的水的渗流规律的渗流就被称为毛细管流。其基本假定为，孔隙中水与裂隙的两壁都是接触的。通常，裂隙的开度会比基质孔隙的孔径大很多，根据毛管理论，因为孔隙的毛细作用比微裂隙的毛细作用大很多，当到达一定的饱和度后，就可忽略微裂隙的毛细作用。这样，吸水时，裂隙中的水首先被吸入基质孔隙中，当孔隙中充满水后，水才能被吸入微裂隙中。而排水时，因为裂隙中产生的毛细吸力比基质吸力小得多，所以，裂隙中的水会先排出，然后基质孔隙中的水才会排出。当然，该假设成立的条件是微裂隙的毛细作用远小于孔隙的毛细作用，并且发生在孔隙与微裂隙都存在的多孔介质中。若微裂隙比较发育，那么其毛细作用会增大。由此可见，对裂隙孔隙岩体持水曲线起主导作用的是岩体基质孔隙特征。b）薄膜流薄膜流是裂隙介质中特有的一种流体流动现象。当非饱和裂隙中水的饱和度较18 三峡大学硕士学位论文小，且裂隙中充满空气时，水会沿着裂隙某一边侧壁，以薄膜的形式向下流动，即称其为薄膜流。薄膜流只与裂隙壁某一面接触，在水流和裂隙壁另一面之间是气相，薄膜流在自身重力作用下可以流动，不受裂隙的毛细管力束缚。薄膜流在不同开度大小裂隙中均可以出现，且流动不直接受裂隙开度控制。c）优先流优先流是指溶质沿着某一渗流路径的速度明显高于其它方向，从而快速运移的非均匀流。优先流存在的原因主要包括：①大孔隙(包括裂隙)的存在，②不稳定性水流的发展，③由于倾斜地层面的存在而改变入渗水的方向产生渗流漏斗。根据优先流理论，土体中的水在大小分布不均的孔隙中流动时，孔隙小的土体基质会形成透镜体，这些透镜体会成为水流的障碍，从而引起水流的集中。水流会绕过土体基质而优先进入大孔隙中。同样的，在孔隙裂隙岩体中，裂隙相对于基质孔隙大得多，当流量大于基质孔隙的水力传导能力或局部的非均质性使水流转移到裂隙中时，水就会绕过较小的岩石基质孔隙优先进入裂隙中。且这些水流会受到重力不稳定性或裂隙开度的变化性的影响。当然，流过裂隙的水流也会被吸入干燥的或非饱和基[46]质孔隙中去，但吸入率比较低。2.2.2单裂隙渗流基本理论裂隙的导水能力比其周围岩石基质通常大几个数量级，它所构成的裂隙网络是岩体中水的主要通道，极大的影响着裂隙岩体的渗流性质。所以，单裂隙的水力特性是[26]研究裂隙岩体渗流的起点。2.2.2.1立方定理平行板模型假定岩石裂隙是由两片光滑的平行板构成的裂隙，隙宽a为常数（如图a/22.9）。缝隙中的水流运动符合Navier-Stokes方程，即a/2∂u1i=−Fpvii,,+uijj（2.35）图2.9平行板缝隙水流流速分布图∂tρ式中，u——流速分量；iF——作用力，对缝隙水流F=0；iiρ——水的密度；p——水的压力；v——水的运动粘滞系数，其值与温度有关。19 三峡大学硕士学位论文当流速很小时，光滑平行板缝隙中的水流为层流，有u=0，也即z方向（垂直3于缝隙面）流速为零。令水力梯度最大方向为x向，则问题可简化为研究x方向流速的一维问题。因隙宽a为定值，则沿x方向各点u也为定值。对于恒定流，式（2.35）x写为：2du1dpx=（2.36）2dzρvdx缝隙内流速一般很小，其流速水头常可忽略。水力势（水头）即为位置水头与压力水头之和，即phz=+（2.37）ρg则式（2.36）可化为2duρdhx=（2.38）2dzvdx由于隙宽a为定值，因而水力梯度为dhJ==常值dx当za=±/2时，u=0，由式（2.38）求积分可得流速按抛物线分布，即x22ga(4)−zuJ=（2.39）x8v通过缝隙的流量为3a/2gaqu==2dzJ（2.40）∫x012v由式（2.40）可知，通过等宽缝隙的流量q与隙宽a的3次方成正比，即著名的立方定理。将其写成达西定律形式为qk=−aJ（2.41）f式中，k为缝隙的水力传导系数。f2gak=（2.42）f12v将式（2.42）写成如下形式ρgkk=()（2.43）finfμ2a其中，()k=（2.44）inf12式中：()k——缝隙内在水力传导系数，它只和缝隙的几何要素——隙宽a的平方成inf20 三峡大学硕士学位论文线性关系。2.2.2.2缝隙水流的水力损失由式（2.38）可求得缝隙水流的平均流速为2qgauJ==（2.45）av12对于管道水力学，长度l的水力损失的一般表达式为2luh=λ（2.46）mdg2式中，λ为阻力系数；d为管内径。仿照式（2.46），对于缝隙水力学，水力损失表达式可写为2luh=λ（2.47）mag2选定一无量纲数，即著名的雷诺数Re，设Re=uav/（2.48）则由式（2.45）、式（2.47）、式（2.48）可求得缝隙流为层流流态的阻力系数，即λ=6/Re（2.49）管道水力学水流一般为紊流流态，其水力损失与平均流速u的平方成正比[式（2.46）]。在层流流态区，式（2.49）中阻力系数含有u的负一次方，因而水力损失与平均流速u的一次方成正比例。无论是管道水力学还是缝隙水力学，如为层流流态，其水力损失均与u的一次方成比例。2.2.2.3立方定理适用范围Ломизе1951年对缝隙水力学进行了系统研究，并进行了光滑和粗糙平行板缝隙水力学实验，将9种不同隙宽的光滑平行板缝隙（am=0.50~5.23m）水力学实验成果整理后得出：当Re<500时，不同隙宽缝隙的试验成果与式（2.49）符合得较好；当Re>600时，λ与Re的关系为0.25λ=0.056/Re（2.50）1.752上式表明，当Re>600时，水力损失与u成比例，是阻力由u向u的过渡区。并且试验表明，当Re100000>时，式（2.50）已不再适用（Агроскин,1950）,水流进入平方阻力区。立方定理是按光滑平行板缝隙为层流推到出来的。据Ломизе的试验，其适用范围为Re<500。将式（2.45）的平均流速代入式（2.49），得21 三峡大学硕士学位论文13Re=gaJ（2.51）212v2设岩石中水的温度为15℃，相应的vc=0.0114m/s。若水力梯度为1，则由Re=500可求得隙宽am=0.927m，即隙宽不大于1mm，水流为层流状态。岩石中的裂隙隙宽一般较小，通常小于0.1mm，因此对于岩石中的大多数的小裂隙，立方定理是成立的。对于隙宽大于1mm的大裂隙，当水力梯度较大时，水流进入由层流向紊流的过渡区，立方定理则不再适用。2.2.2.4立方定理的修正现实中裂隙壁是非光滑的，如果将基于光滑平行板缝隙层流导出的立方定理用于实际裂隙，则必须对立方定理进行修正。实际的裂隙壁很粗糙且凹凸不平，会致使水流曲线不再是一条直线，与光滑平行板模型中的水流相比，流线加长了，那么同一水力梯度下，因为阻力系数加大了，同一隙宽的粗糙裂隙比光滑缝隙通过的流量要小。Ломизе经过粗糙裂隙水力试验得出，裂隙面的凸起度∆与隙宽的比值对裂隙过流能力影响较大。故将立方定理加入了一个粗糙度修正系数C，则3gaqJ=（2.52）12vC根据试验结果Ломизе认为C值与相对起伏差∆/a有关，可用下式表示1.5⎛⎞ΔC=+1⎜⎟（2.53）⎝⎠aLouis（1967）根据自己的试验资料，建议C值用下式计算1.5⎛⎞ΔC=+18.8⎜⎟（2.54）⎝⎠2a1975年，Жцленков给出了粗糙度修正系数的新公式5(λ−1)2C=+1（2.55）a式中：λ为裂隙面曲面面积与其投影面积之比。2除上述各修正公式以外，Neuzil、Tsang、Walsh、Barton和ElsworthandGoodman等也通过试验或理论研究提出了各自的修正公式。2.2.3裂隙岩体土水特征曲线在裂隙岩体中，由于很难通过试验来测定准确度较高的非饱和裂隙岩体土水特征曲线，所以一般的非大型的岩体工程并不适合采取试验的方法，而是采用更为简单的方法来得到较精确的结果。在上节中已经提到，Fredlund和Xing提出了根据土体孔隙孔径分布曲线来得到土水特征曲线的理论，后Zhang和Fredlund又将这一理论延伸到非饱和裂隙岩体中，22 三峡大学硕士学位论文[47]以得到非饱和裂隙岩体的土水特征曲线。自然界中，岩体中裂隙的发育程度不一。当裂隙开度小于0.25mm时就会产生毛细水，而当开度大于0.25mm时则可不考虑毛细水的存在。裂隙中的毛细作用机理与基质孔隙中的毛细作用机理是一样的。根据毛管理论，吸力与空气—水界面的曲面半径成反比关系，在裂隙中有：2Th=cosα（2.56）cγxw式中：h——毛细吸力水柱高度；cγ——水的容重；wx——裂隙开度的半径；T——孔隙水表面张力；α——接触角。当裂隙中存在进气压力u时，上式可表示为：a2Tψ=−=()cuuosα（2.57）awx由此，基于基质孔隙孔径分布的土水特征曲线估算方法同样适用于裂隙。含水率可写为：Rθθ()R=frdr()（2.58）sf∫Rmin式中：θ()R——张开度小于或等于R的裂隙中充满水时的水含量；θ——裂隙的饱和含水率；sfR——最小的裂隙张开度。min孔隙裂隙岩体的土水特征曲线定义为裂隙岩体饱和度或含水率与裂隙岩体中毛细吸力之间的关系曲线，它是研究非饱和岩体水力特征参数的基础。方程(2.58)对于任何能用毛管理论描述的多孔介质都有效(Zhang和Fredlund)。因此也同样可用于含微裂隙的岩体中。孔隙裂隙岩体是由岩体基质孔隙和裂隙组成的双重介质，且一般情况下微裂隙的张开度比基质孔径大很多。根据毛细管理论，当裂隙张开度远大于孔隙孔径时，孔隙介质的毛细作用比裂隙介质大得多，而在达到一定饱和度以后，可基本忽略裂隙的毛细作用。根据上述理论，可用下述两种方法来模拟孔隙裂隙介质的土水特征水曲线：a）将基质孔隙和裂隙的土水特征曲线作叠加考虑，既考虑孔隙介质的毛细作用，23 三峡大学硕士学位论文又考虑裂隙的毛细作用。计算时，先分别得到裂隙和孔隙的土水特征曲线，然后将两条土水特征曲线迭加，即为综合土水特征曲线。由综合土水特征曲线得到的水力特征参数即为综合水力特征参数。b）考虑到裂隙张开度比孔隙孔径大得多，则只考虑孔隙介质的毛细作用，忽略裂隙的毛细作用。那么，裂隙孔隙岩体的综合土水特征曲线被认为等同于孔隙介质的土水特征曲线。Zhang和Fredlund(2003)将方程(2.58)应用到裂隙孔隙岩体中，先对岩体的基质孔隙和微裂隙分别计算，再将两者叠加考虑，得到了最终的非饱和裂隙岩体土水特征曲线。叠加考虑后的土水特征曲线估算方程可表示为：RRmfθψ()=+θf()rdrθfrdr()（2.59）sm∫∫msfRRminmin式中：θ——基质孔隙的饱和含水率；smθ——裂隙的饱和含水率；sff()r——基质孔隙孔径的分布函数；mf()r——裂隙张开度的分布函数；fR——基质孔隙孔径，mR——裂隙张开度。fZhang和Fredlund的研究表明：由于岩体中的裂隙尺寸比基质孔隙大很多，裂隙比基质孔隙先达到非饱和状态，当吸力超过裂隙中的残余吸力后，饱和度—吸力关系曲线主要受岩体基质孔隙的影响，只有当岩体接近饱和时，岩体中裂隙起到的作用才会有明显增大。所以裂隙孔隙岩体介质与单孔隙介质具有相似的饱和度—吸力关系曲线。在土水特征曲线的估算方程中，vanGenuchten方程和Fredlund和Xing方程能很好的适应各种类型的岩土介质。mθθ−⎡1⎤wra）VanGenuchten模型：Θ==⎢⎥nθθ−+⎣1()αψ⎦srm⎡1⎤b）Fredlund和Xing模型：θθ=⎢⎥sn⎢ln⎡ea+(/)ψ⎤⎥⎣⎣⎦⎦两个方程均为三参数拟合方程，vanGenuchten模型与Fredlund和Xing模型都包括a,m,n三个拟合参数。根据Fredlund和Xing(1994)的研究：①当m,n固定不变时，参数a与进气值紧密相关；②参数n反映了曲线的坡度；③参数m控制着曲线的平24 三峡大学硕士学位论文滑度，m越小曲线越平滑。这两个方程的主要不同之处在于vanGenuchten方程在吸力值达到很高的时候趋于一条水平线，而Fredlund和Xing方程能满足从0~1GPa的吸力范围。2.2.4裂隙岩体渗透性函数在岩体中，非饱和现象及毛细吸力主要是由于岩体基质孔隙和岩体微裂隙的存在而共同作用产生的，孔隙裂隙岩体可以看作是由岩体颗粒的空隙和裂隙组成的双重介质。而用试验方法较准确地测定具有基质孔隙和裂隙双重介质的岩体的渗透性函数是一件较困难的事，但目前也没有成熟的确定方法。根据裂隙非饱和持水机理与多孔介质的相似性原理，孔隙裂隙岩体非饱和渗透系数可以按土体非饱和渗透系数确定方法进行估算（见2.1.3节）。2.3本章小结在本章中，分别就孔隙介质和裂隙介质的渗流特性做了总结。孔隙介质的渗流特性与其内部水的存在形式与特征有关，并主要反映在土水特征曲线与渗透性函数的变化上。这里也介绍了土水特征曲线及渗透性函数的实验方法和理论预测。而关于裂隙介质，重点介绍了单裂隙渗流的基本理论，以及裂隙多孔介质的土水特征曲线和渗透性函数的预测。25 三峡大学硕士学位论文3裂隙网络的随机模拟离散裂隙网络模型假定，裂隙的渗流可用裂隙几何形状和导水率进行预测。而生成离散裂隙网络模型，必须对裂隙的几何参数予以测量，且对这些数据进行统计分析，以获得裂隙的统计参数及其所服从的概率分布函数。在此基础上，用计算机生成在统计意义上等效的二维或三维裂隙网络。本文中只限于二维网络的生成。且假定每个单裂隙有规则的几何形状和相同的隙宽。且不考虑在两条裂隙交点上的水头损失。3.1裂隙的基本特征3.1.1裂隙的产状裂隙的产状通常用走向和倾角来定义。走向是对裂隙与水平面的交线的描绘，其方向可以用交线的方位角来确定，也即由北顺时针旋转的角度。如方位角320゜则对应方向N40W。倾角指的是裂隙面与水平面间的夹角。在二维裂隙网络中，裂隙的方向可只用一个角度来表示。3.1.2裂隙尺寸迹线是指岩石露头、边坡或地下洞室的开挖面上观察到的裂隙面与出露的岩面的交线。就目前的技术条件，很难从现场获得裂隙面大小的准确信息，只能取得迹长的量测数据。而所谓裂隙尺寸，指的是通过迹长的测量与概率统计的方法来求得的裂隙迹长或其半径的统计参数及其所服从的分布规律。3.1.3裂隙密度（间距或频率）裂隙间距系指同一组裂隙（其产状围绕其优势产状呈随机变化，但幅度不大）的垂直距离。同一组各裂隙的间距也是不等的，通常用均值间距表示统计意义上的平均间距。同一组裂隙因其尺寸及其在空间的分布也是随机的，对一条裂隙迹线，在其端点作垂直裂隙的延长线，可能并不与距离最近的另一条裂隙相交。因此，裂隙间距不能恰当反映裂隙的密度，文献中常用裂隙频率来描述裂隙的密度。裂隙频率有几种定义：①线频率f1，即与某一方向单位长度测线相交的同一组裂隙的数目；②面频率f2，即某一平面单位面积上与同一组裂隙迹长相交的数目；③体积频率f3：即单位体积内所包含的与同一组裂隙相交的数目。3.1.4裂隙开度本文没有直接使用裂隙开度，而是对裂隙指定了一个常数导水率(m2/s)。因为裂隙开度通常被认为服从此分布(Snow,1970)，所以这里也假定其服从对数正态分布。26 三峡大学硕士学位论文3.2实验室裂隙生成实验试验目的：本次试验的目的是通过室内裂隙开展试验来取得裂隙相关参数的统计分布规律。试验设备：不锈钢铁盘6个，温度计1支，相机1台，秒表1只，支撑架1件（实验室自制）。试验方案：制作31cm长22cm宽1cm厚的土膏，在自然风干作用和人工烘干作用两种条件下的裂隙开展试验。观察裂隙发展规律，并用相机记录其开展过程。最后统计裂隙的发展规律。试验过程：将盛放土膏的铁盘、相机的相对位置调整好，固定在自制支撑架上，在自然风干条件下，或者烘干条件下，每隔一定的时间，拍摄土膏的开裂状态，直至裂隙开展速度趋于0，结束拍摄，整理相片，统计裂隙相关参数的规律。本次裂隙开展实验共计六组，取其中一组风干实验过程进行分析。取裂隙发展全过程的三个时段（如表3.1），统计了四个时刻时裂隙的相关参数的统计分布规律，以及每时段内，裂隙发展的平均速度。表3.1裂隙实验基本资料平均速度编号拍照时间∆t(min)（mm/min）74842009.9.2321：2674922009.9.2411：278410.011575022009.9.2422：166490.044575032009.9.2507：585820.0221如图3.1所示，初始时刻，整块土样上没有裂隙。经过841min的间隔后，土样表面产生了3.1b图所示的少量裂隙，且裂隙开展的平均速度较慢。与上一时刻相隔649min后，土样表面的裂隙如图3.1c所示，这一时段内裂隙发展的平均速度现对较快。由图3.1c到图3.1d经历了582min，裂隙开展的速度明显降低，将时刻d作为终了时刻。27 三峡大学硕士学位论文abcd图3.1裂隙发展过程下面将相关参数：裂隙中点坐标X、Y、角度a以及长度L的统计规律分时段进行分析。初始时刻a，土样还没有产生裂隙。28 三峡大学硕士学位论文时刻b，裂隙的各参数柱状统计图如下：图3.2a裂隙中点坐标X、Y统计分布柱状图（时刻b）图3.2b裂隙长度L和lnL统计分布柱状图（时刻b）图3.2c裂隙角度a统计分布柱状图（时刻b）图3.2d裂隙角度a分组后各组数据统计分布柱状图（时刻b）29 三峡大学硕士学位论文表3.2裂隙中点坐标X、Y及长度L的统计分布规律拟合（时刻b）XYLDistributionADPLRTPADPLRTPADPLRTPNormal8.082<0.0053.928<0.0055.666<0.005Lognormal4.126<0.00517.786<0.0050.9010.0213-ParameterLognormal1.776*03.972*00.547*0.097Exponential2.4770.00314.695<0.00318.608<0.0032-ParameterExponential2.513<0.0100.22214.709<0.0100.18411.32<0.0100Weibull1.083<0.0108.028<0.0101.13<0.0103-ParameterWeibull1.0750.0090.8094.133<0.00500.4820.2390SmallestExtremeValue13.113<0.0104.866<0.01022.134<0.010LargestExtremeValue3.811<0.0104.86<0.0100.9440.018Gamma1.0940.019.962<0.0050.382>0.2503-ParameterGamma1.094*0.9513.998*00.285*0.057Logistic6.896<0.0054.037<0.0053.174<0.005Loglogistic2.851<0.00511.934<0.0051.072<0.0053-ParameterLoglogistic2.124*0.0184.041*00.982*0.626表中，当P值大于0.005时，则表示数据服从该统计分布，且P值越大符合程度越好。而LRTP值则用来说明同种分布下两参数的分布是否比三参数的拟合程度更高。根据上表3.2中的p值来看，裂隙中点横坐标X的分布与指数分布（Exponential）（p=0.003）、3参数韦布分布（3-ParameterWeibull）（p=0.009）以及Gamma分布（p=0.01）接近，且三种分布比较可得，X的分布与Gamma分布更为吻合。纵坐标Y拟合后所得的p值，均小于0.010，所以与表3.2中所列的统计分布都不符合；而根据其柱状分布图，可以得出纵坐标Y基本符合均匀分布。同样的，根据裂隙长度L拟合后得出的p值，可以看出其与对数正态分布（Lognormal）（p=0.021）、3参数韦布分布（3-ParameterWeibull）（p=0.239）、最大极值分布（LargestExtremeValue）（p=0.018）以及Gamma分布（p>0.250）接近，比较得，裂隙长度L与Gamma分布最为吻合。30 三峡大学硕士学位论文表3.3裂隙角度a及其分组后的统计分布规律拟合（时刻b）a(分组前)a(a≤90゜)a(a>90゜)DistributionADPADPADPLRTPNormal5.221<0.0052.937<0.0052.125<0.005Lognormal----2.961<0.0053-ParameterLognormal3.951*3.016*2.173*0.001Exponential----30.926<0.0032-ParameterExponential7.491<0.0106.443<0.0109.72<0.0100Weibull----1.86<0.0103-ParameterWeibull3.954<0.0053.353<0.0051.583<0.0050.16SmallestExtremeValue7.662<0.0103.171<0.0101.577<0.010LargestExtremeValue3.991<0.0103.293<0.0103.379<0.010Gamma----2.669<0.0053-ParameterGamma4.315*3.603*2.425*0.058Logistic5.103<0.0053.026<0.0052.022<0.005Loglogistic----2.597<0.0053-ParameterLoglogistic4.021*3.059*2.039*0.006在对角度a的统计过程中，实行了不分组统计和根据a是否大于90゜的分组统计两种统计方式。根据上表中p值可以得出，裂隙的方向角a及其分组后的统计均不符合表中所列出的统计分布规律，那么根据其柱状分布图来看，角度a是大致符合均匀分布的。31 三峡大学硕士学位论文图3..3裂隙参数X、L拟合分布图（时刻b）图3.3显示了上述变量的拟合中，各自变量与其分布相符合的情况。图中，数据点越靠近中心线，拟合程度越高，中心线两边的线条分表代表了可偏离的最大程度，超出这两边的线条，则表示不符合该种分布。时刻c，裂隙的各参数柱状统计图如下：32 三峡大学硕士学位论文图3.4a裂隙中点坐标X、Y统计分布柱状图（时刻c）图3.4b裂隙长度L和lnL统计分布柱状图（时刻c）图3.4c裂隙角度a统计分布柱状图（时刻c）图3.4d裂隙角度a分组后各组数据统计分布柱状图（时刻c）33 三峡大学硕士学位论文表3.4裂隙中点坐标X、Y及长度L的统计分布规律拟合（时刻c）XYLDistributionADPLRTPADPLRTPADPLRTPNormal12.944<0.00511.071<0.00526.083<0.005Lognormal35.84<0.00553.703<0.0051.805<0.0053-ParameterLognormal10.277*011.308*01.451*0.242Exponential23.697<0.00338.177<0.00326.731<0.0032-ParameterExponential23.632<0.0100.06238.173<0.0100.12910.706<0.0100Weibull13.046<0.01025.624<0.0103.087<0.0103-ParameterWeibull10.47<0.0050.00212.874<0.00500.5830.1360SmallestExtremeValue19.414<0.01012.146<0.01081.686<0.010LargestExtremeValue10.339<0.01014.841<0.0107.291<0.010Gamma15.663<0.00530.133<0.0051.604<0.0053-ParameterGamma11.826*0.00911.667*00.254*0Logistic13.06<0.00511.486<0.00516.764<0.005Loglogistic22.959<0.00537.54<0.0052.506<0.0053-ParameterLoglogistic10.826*011.631*02.747*0.289上表3.4中裂隙中点横坐标X拟合得到的p值、纵坐标Y拟合得到的p值均小于0.010，那么它们均不符合表中所列的统计分布规律。而跟据各自的柱状分布图可以得出，横坐标X、纵坐标Y都大致符合均匀分布。同样，由裂隙长度L拟合得到的p值来看，其规律与3参数韦布分布（3-ParameterWeibull）（p=0.136）符合。拟合分布图如图3.5所示。图3.5裂隙参数L拟合分布图（时刻c）34 三峡大学硕士学位论文表3.5裂隙角度a及其分组后的统计分布规律拟合（时刻c）a(分组前)a(a≤90゜)a(a>90゜)DistributionADPADPADPLRTPNormal6.865<0.00510.147<0.0056.685<0.005Lognormal----7.719<0.0053-ParameterLognormal6.96*10.363*6.751*0.004Exponential----105.834<0.0032-ParameterExponential45.59<0.01021.926<0.01017.57<0.0100Weibull----6.345<0.0103-ParameterWeibull7.72<0.00511.035<0.0057.921<0.0050SmallestExtremeValue11.23<0.0109.476<0.0106.823<0.010LargestExtremeValue10.913<0.01011.68<0.0108.264<0.010Gamma----7.248<0.0053-ParameterGamma7.265*30.869*7.092*0.465Logistic7.34<0.00510.025<0.0056.756<0.005Loglogistic----7.59<0.0053-ParameterLoglogistic7.381*10.177*6.776*0.012根据上表中p值可以得出，裂隙的方向角a及其分组后的统计均不符合表中所列出的统计分布规律，那么根据其柱状分布图来看，角度a是大致符合均匀分布的。终了时刻d，裂隙的各参数柱状统计图如下：35 三峡大学硕士学位论文图3.6a裂隙中点坐标X、Y统计分布柱状图（时刻d）图3.6b裂隙长度L和lnL统计分布柱状图（时刻d）图3.6c裂隙角度a统计分布柱状图（时刻d）图3.6d裂隙角度a分组后各组数据统计分布柱状图（时刻d）36 三峡大学硕士学位论文表3.6裂隙中点坐标X、Y及长度L的统计分布规律拟合（时刻d）XYLDistributionADPLRTPADPLRTPADPLRTPNormal13.889<0.00514.039<0.00533.318<0.005Lognormal45.126<0.00553.69<0.0051.55<0.0053-ParameterLognormal12.458*014.262*01.525*0.848Exponential34.186<0.00332.649<0.00326.057<0.0032-ParameterExponential34.092<0.0100.05432.669<0.0100.0978.639<0.0100Weibull17.248<0.01026.542<0.0104.473<0.0103-ParameterWeibull12.483<0.005019.168<0.00500.9350.020SmallestExtremeValue19.924<0.01016.316<0.01096.595<0.010LargestExtremeValue12.483<0.01016.54<0.01010.78<0.010Gamma21.277<0.00529.143<0.0052.821<0.0053-ParameterGamma12.74*027.704*0.4920.466*0Logistic14.199<0.00514.381<0.00521.843<0.005Loglogistic28.527<0.00538.903<0.0052.739<0.0053-ParameterLoglogistic12.818*014.561*02.976*0.033上表3.6中裂隙中点横坐标X拟合得到的p值、纵坐标Y拟合得到的p值均小于0.010，那么它们均不符合表中所列的统计分布规律。而跟据各自的柱状分布图可以得出，横坐标X、纵坐标Y都大致符合均匀分布。同样，由裂隙长度L拟合得到的p值来看，其规律符合3参数韦布分布（3-ParameterWeibull）（p=0.2）如图3.7所示。图3.7裂隙参数L拟合分布图（时刻d）37 三峡大学硕士学位论文表3.7裂隙角度a及其分组后的统计分布规律拟合（时刻d）a(分组前)a(a≤90゜)a(a>90゜)DistributionADPADPADPLRTPNormal6.786<0.00512.724<0.0057.087<0.005Lognormal----7.678<0.0053-ParameterLognormal7.015*13.019*7.142*0.048Exponential----119.003<0.0032-ParameterExponential54.374<0.01026.06<0.01018.368<0.0100Weibull----6.911<0.0103-ParameterWeibull7.767<0.00513.2<0.0058.052<0.0050SmallestExtremeValue11.344<0.01011.498<0.0108.011<0.010LargestExtremeValue12.778<0.01014.764<0.0108.185<0.010Gamma----7.331<0.0053-ParameterGamma7.701*38.744*8.22*0.386Logistic7.302<0.00512.39<0.0057.221<0.005Loglogistic----7.706<0.0053-ParameterLoglogistic7.454*12.598*7.226*0.115根据上表中p值可以得出，裂隙的方向角a及其分组后的统计均不符合表中所列出的统计分布规律，那么根据其柱状分布图来看，角度a是大致符合均匀分布的。综合上述统计数据，裂隙中点横坐标X、纵坐标Y、长度L及角度a的统计分布规律在三个不同时刻的分布情况如下：表3.8裂隙各参数统计分布规律随时间的变化时刻b时刻c时刻dX伽马分布均匀分布均匀分布Y均匀分布均匀分布均匀分布L伽马分布3参数韦布分布3参数韦布分布a均匀分布均匀分布均匀分布裂隙中点横坐标X的统计分布规律在第一时段内是符合指数分布、3参数韦布分布、伽马分布的，并以伽马分布最为吻合；随着裂隙的开展，第二时段中规律趋于均38 三峡大学硕士学位论文匀分布；第三时段也稳定在均匀分布上。裂隙中点纵坐标Y的规律随时间的推移并没有多大改变，三个时段都呈现出均匀分布的规律。裂隙的长度L第一时段统计分布规律是符合对数正态分布、3参数韦布分布、最大极值分布、伽马分布的，并以伽马分布最为吻合；随着裂隙的进一步开展，第二时段内分布规律趋于3参数韦布分布；第三时段也稳定在该分布上。裂隙的角度a的规律随时间的变化不变，三个时段内的统计规律都符合均匀分布。3.3基于LHS方法的离散裂隙随机生成随机模拟的概念是和离散裂隙网络模型有密切联系的。可以通过岩体中的每组裂隙几何性质的统计描述来表征网络的几何形状。然后，再根据该统计模型，生成多个离散裂隙网络。本文中将基于拉丁超立方抽样法（LHS）来生成二维裂隙网络。3.3.1拉丁超立方抽样法（LHS）从假定的裂隙的各个几何参数的概率模型分布中进行随机抽样，来得到单个裂隙的具体参数，然后，将这些参数进行组合，实现一个具体的裂隙。从已知的概率分布中的随机抽样即可用拉丁超立方抽样方法模拟实现。下面从分层抽样开始，来了解拉丁超立方抽样。A）分层抽样我们考虑一维的单个变量输入问题：y=fx()，x是一个随机变量。分层抽样通过如下的步骤来进行：1)定义参与计算机运行的抽样数目N；2)将x等概率地分成若干个区域——“bin”（图3.8）,图3.8bins示意图x<<<<<