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

三维边坡渗流场计算的变单元法 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)步中自由面出现回弯的问题。称此方法为“修正的变单元渗透系数法”。运用实际的非饱和区的渗透系数来计算渗流场,比假设该区为“干区”具有更严格的理论基础。 经过上述计算,得出的渗流场,即为稳定状态下饱和-非饱