三维边坡渗流场计算的变单元法

三维边坡渗流场计算的变单元法

1有限渗流的求解方法在研究边缘、基础和水库的变形和稳定方面，渗流计算发挥着重要作用。自由面是渗流问题研究的重点和难点之一,渗流自由面的确定,可以明确工程中岩土体的赋存特性。有限单元法和有限差分法是用来求解渗流问题的两种常见方法。变网格法由于要不断修改网格形状,操作相当不便,基本已经被淘汰;固定网格法中常用的方法包括剩余流量法、初流量法、Signorini型变分不等式法和渗透系数调整法等。上述方法在应用到有限差分法来求解渗流自由面的研究相对于有限单元法较少,但有限差分法可以对不同复杂几何边界下的物理问题实现统一的数值求解,程序算法设计简单。在渗流分析,尤其自由面渗流问题的求解中可以很好地避免传统有限元法的缺点。FLAC3D是国际通用的三维岩土力学有限差分计算机程序,它提供了基于fish语言的二次开发平台,尤其擅长流-固耦合问题的求解。传统的渗流分析中,认为自由面以上的区域是干燥的,在自由面上没有渗流作用,然而随着饱和-非饱和理论的发展和完善,非饱和区负孔隙水压力对边坡、大坝等工程的有利影响得到充分认识。然而,在FLAC3D中不能直接求解出渗流的自由面和溢出点,并且其渗透系数只能设置为一固定常数。为了克服上述困难,本文利用fish语言,对FLAC3D的渗流模块进行二次开发,运用修改渗透系数法求解出渗流的自由面和溢出位置。在此基础上,利用非饱和土的特性,不断迭代修改非饱和区的渗透系数,进而算出较为精确的饱和-非饱和渗流场。2计算原则和步骤2.1面水渗流数值模拟非饱和土中的渗流支配方程,在各向同性均质情况下为式中:渗透系数kw为体积含水率wn的函数;Sw为饱和度;Ss为单位贮存量。式(1)是进行饱和-非饱和渗流计算得一般微分方程式。渗流自由水面以上的含水率状态受制于非饱和流动,降雨或灌溉水入渗地下,通常是非饱和渗流。对于土坝和土堤来说实际上也受非饱和渗流的影响。因此,Neuman提出,同时考虑饱和与非饱和进行数值计算得数学模型。此方程可以用压力水头p代替水头h,以体积含水率wn代替饱和度Sw;并定义单位容水量,则式(1)可写为式中:压力水头p是体积含水率wn的函数,在已知的初始和边界条件下即可求解上式,连接计算所得p(28)0各点就是渗流自由面。2.2基于傅氏一元模型的三维地震分析求解饱和非饱和渗流场,最重要的是要求出渗流场饱和区与非饱和区的分界面,即自由面。其次,饱和区的体积含水率不发生变化,因此,渗透系数为常数;而非饱和区的渗透系数,对于不同的土体,会随土体类型、颗粒大小、孔隙率、体积含水率的变化而变化。在实际中,非饱和区的渗透系数,是通过土-水特征曲线试验得出,这里采用Garden建议的关于渗透系数作为基质吸力的函数关系式以作方法上介绍,后人也提出了许多关于基质吸力与渗透系数相互关系的理论,并经试验证明其有效,对于具体的工程,应根据试验来得出各土、岩层的曲线。计算的具体步骤如下:(1)可直接用fish语言建模或从ANSYS等软件中先建立模型然后导入到FLAC3D中,方法简单,建模方法见陈育民等,限于篇幅,不予介绍。由于渗流溢出的位置开始无法确定,首先只定义上、下游水头边界和全区域材料的渗透系数都定义为饱和状态下的渗透系数进行计算。根据经验,由此得出的孔隙水压力为0的面(初始自由面)是要高于真实自由面的,所以自由面的位置处于上、下游水头连线和初始自由面之间,溢出位置位于初始溢出位置以下和下游水头以上。(2)由第(1)步计算得出的渗流场,通过fish语言,寻找孔隙水压力为负值的区域,将该区域渗透系数修改成为远小于饱和状态下渗透系数的值,比如:1×10-18,以保证饱和区与非饱和区的流量交换几乎可以忽略。重复上述迭代(程序默认的收敛精度为1×10-5),直到下一次迭代的步数小于10步时,停止迭代,这样前后两次计算得各点的水头差甚微,完全能满足精度要求,由此计算得出自由面和溢出点的位置。但由于此时从饱和过渡到非饱和的区域的渗透系数变化过大,所以浸润线的形状扭曲,而且有回弯,根据陈益峰等的研究,在连续均质介质中这种情况是不会出现的。(3)在第(2)步的基础上,寻找下游面溢出边界上的节点,重新定义其位置水头为其总水头或定义其孔压为0,根据渗透系数和基质吸力的关系,重新修改非饱和区域的渗透系数(对于非均质、不同土层,可循环修改各位置、各方向的渗透系数,并将在三维渗流分析中体现,方法简便、耗时少),进行迭代计算,收敛准则与第(2)步中相同。计算后所得浸润线较为光滑,解决了第(2)步中自由面出现回弯的问题。称此方法为“修正的变单元渗透系数法”。运用实际的非饱和区的渗透系数来计算渗流场,比假设该区为“干区”具有更严格的理论基础。经过上述计算,得出的渗流场,即为稳定状态下饱和-非饱和渗流场。3计算结果对比为了验证上述方法的正确性,本文采用基于有限单元法的GEO-SLOPE的计算结果与FLAC3D的计算结果进行对比。GEO-SLOPE是一个二维的岩土工程仿真分析软件,其渗流模块SEEP/W最大的优势就是能够定义潜在溢出位置并自动搜索溢出点,应用十分广泛;但它只局限于二维计算,不能进行三维分析,而往往工程中的问题,不能简化为二维问题。3.1孔隙水压力下降的过程本文采用一个宽12m,高16m,厚1m的模型,上游水头为12m,下游水头为4m,假定饱和状态下的渗透系数为1×10-7,将模型剖分成192个网格(12×16×1),如图1所示。定义上、下游水头边界条件,全部区域都定义饱和状态的渗透系数进行计算,找出孔隙水压力为0的等值线,即为求出的初始浸润线,如图2所示,该线以上的区域为非饱和区,以下的区域为饱和区。然后,利用fish语言,找出孔隙水压力小于0的区域,即非饱和区,保持饱和区的渗透系数不变,重新赋值非饱和区的渗透系数kuf0a2(28)1×10-18,以保证饱和区和非饱和区的流量交换很小,可以忽略不计,在循环迭代过程中,不断重复上述修改渗透系数的步骤,如果下一次迭代的步数低于10步,则计算停止这样前后两次计算得各点的水头差甚微。如图3所示,孔隙水压力为0的等值线(浸润线)以上的区域(非饱和区),由于渗透系数很小,基本没有渗流,通过此步骤的计算,得出了弯曲的浸润线,这是由于两个区域所定义的渗透系数突变,而发生了局部水流折射所致,但根据自由面在任一均匀介质内部必连续光滑且不产生回弯的原理,该浸润线必须进行修正。在上述计算的基础上,记录下游溢出边界上的节点,并重新赋值其孔隙水压力为0。从图中粗略地可以看出,溢出边界为z从4m到7m的区域。然后利用渗透系数和基质吸力的关系,重新修改非饱和区域的渗透系数,进行迭代计算,在循环迭代过程中,不断重复上述修改渗透系数的步骤,如果下一次迭代的步数低于10步,则计算停止。最终计算结果如图4所示,其浸润线相当光滑,没有出现回弯,这是因为从饱和区到非饱和区域的渗透系数是一个逐渐变化的过程,没有出现跳跃,这个方法成功地实现了这个平缓的过程,与实际情况相符。对于非饱和渗流区域,从下至上流速矢量由于孔压的降低而不断减小,在最上部终趋近于0。3.2两种算法的对比利用基于有限单元法的GEO-SLOPE程序,计算该模型的饱和-非饱和渗流场,如图5所示,定义同样的上下游水头边界,在下游面上定义潜在溢出边界,范围设置为上、下游水头之间。然而,该程序在定义其渗透系数时,需要取不同孔压值时的渗透系数做代表点,再自动生成曲线,当取的代表点不够时,可能引起非饱和区的渗流场有较大的误差。因此,该程序计算非饱和区渗流时存在不完善性。为了精确地进行对比,在运用GEO-SLOPE进行计算时,定义了30个代表点,以保证其准确性,这个步骤相对比较麻烦。而利用fish语言,可以简便、准确地描述渗透系数和基质吸力之间的关系。经过计算,得出的渗流场如图6所示。将浸润线,流速矢量,压力水头等值线绘制于图中,与图5对比,各渗流要素的形式基本一致,尤其是浸润线的位置,极其相近。通过fish调出FLAC3D计算结果中每个节点的孔隙水压力,与GEO-SLOPE中计算得每个节点的孔隙水压力进行对比分析,取下游面上所以的节点为代表点,将其结果列于表1中。从表1和图8的结果来看,两种算法计算出的各点的孔压值相差甚小。取从8m高程到16m高程的9个节点的孔压差,求其算术平方和为369Pa,而上下游水头差为8m,相当于78.4kPa,这个孔压差只占它的0.471%,这说明,本文提出的方法与GEO-SLOPE计算结果相近,具有有效性。4复杂flac3的接触面单元为了与二维的计算软件GEO-SLOPE进行对比分析,上述计算采用的是单位厚度,模型简单,所以在第(3)步中,获得溢出位置以下、下游水头以上的节点比较容易,手动也可以操作。但对于复杂的三维模型,尤其是边坡等表面弯曲的模型,由于其分水岭处的上游边界水头不是水平的,溢出边界的上边缘常常是曲线,手动寻找这些点是非常困难的,工作量极大。为了解决此问题,可以在下游水头以上的下游面附着一层接触面。FLAC3D中的接触面采用的是无厚度接触面单元,该单元由一系列的3节点的三角形单元构成,接触面单元将三角形面积分配到各个节点中,每个接触面节点都有一个相关的表示面积。每个四边形区域面用两个三角形接触面单元来定义,然后在每个接触面单元顶点上自动生成节点,当另外一个网格面与接触面单元相连,接触面节点就会产生。接触面单元可以通过接触面节点和实体单元表面(称为目标面)之间建立联系。分析复杂三维状态下FLAC3D的饱和-非饱和渗流时,在寻找自由面以下、下游水头以上的这些节点时,采用附着接触面的方法。在潜在溢出区域附着一层接触面,然后借助fish语言,找到连接接触面的节点中,寻找溢出边界上的节点,然后将这些节点的孔隙水压力重新赋值为0。因为前面提到的所有fish程序都是通用的,因此,只需注意这一步就能按2.2中的步骤进行求解。下面以某水利枢纽工程右岸边坡渗流场为例,验证该方法的可行性。该处地形复杂,大坝右肩边坡位置原为一山脊,河道上方位置弯曲不规则,且边坡有多层较厚的风化层,各层渗透系数皆不相同,边坡渗流场分布情况相当复杂。边坡背面水头约为185m,河道水头为96m,假设两侧为不透水边界。河道上方位置相当部分(范围可覆盖整个表面)假定为潜在溢出面,弯曲不规则,计算网格及接触面示意见图9。通过计算,边坡饱和-非饱和渗流场,其压力水头分布情况见图10。地质勘探测到的3个位置的地下水自由面位置,与计算结果十分相近,相差均不大于0.45m水头,误差由于地质十分复杂等而不可避免。上述结果表明,此方法无论是对简单的模型,还是对复杂的