[发明专利]一种基于实域总势能计算渗流溢出点的方法

摘要

本发明涉及一种基于实域总势能计算渗流溢出点的方法，其属于水利工程、岩土工程、边坡工程等问题的渗流分析计算领域。其实域总势能E实＝E0‑αE过。本发明的有益效果是将虚域势能从全域总势能中扣除，提出了实域总势能的概念，减少了因虚域势能造成的误差，提高了求解溢出点的计算精度；解决了求解溢出点时迭代条件复杂且不易收敛的问题，避免了以往求解溢出点时需同自由面一起迭代，提高了求解效率及计算稳定性；通过求解计算模型表明，可更精确地确定溢出点位置，通过采取预防措施改变渗流场的水头值、渗流方向、渗流高度及溢出点坡降等，减轻渗流对边坡稳定造成地危害，防止滑坡、管涌及流土等渗透破坏的发生。

主权项

一种基于实域总势能计算渗流溢出点的方法，其特征在于，包括如下步骤：⑴支配方程和边界条件根据达西定律和地下水运动的连续性条件，不考虑土和水的压缩性，二维均质各向异性土体的稳定渗流满足偏微分方程∂∂x(kx∂h(x,z)∂x)+∂∂z(kz∂h(x,z)∂z)=0---(1)]]>式中，h(x,z)是水头函数，kx、kz分别为x、z方向的渗透系数；在稳定渗流计算中，渗流溢出点的位置是未知的，且必须同时满足第一类边界条件和第二类边界条件，即①水头等于位置高程，即h＝z；②有外渗流量，即⑵应用有限单元法求解渗流场由变分原理可知，渗流场基本微分方程(1)的定解问题等价于求解渗流能量泛函的极值问题，构造如下泛函：I(h)=∫∫Ω[12kx(∂h∂x)2+12kz(∂h∂z)2]dxdz---(2)]]>分别为x、z方向的渗流速度，分别为x、z方向的水力坡降；根据变分原理，对泛函求极值渗流控制方程转化为如下线性方程组：[K]{h}＝{f}                                              (3)[K]为总渗透矩阵，m为节点数，{h}为待求m个节点的水头向量，{f}为自由项，对于稳定渗流场{f}＝{0}；求解线性方程组(3)即可得到渗流场内各节点的水头值；⑶全域总势能确定溢出点原理在求解全域总势能时，以上游坡面底部A为坐标原点，水平方向向右为x轴正方向，铅垂向上方向为z轴正方向，建立直角坐标系，选取三角形单元对全域进行划分，定义每延米上的单个三角形单元的势能为Ee=nρg∫∫eΣi=1mNihidxdz---(4)]]>式中，Ee为每延米上单个三角形单元的势能(J/m)，ρ为水的密度(kg/m3)，g为重力加速度(m/s2)，Ni为三角形单元的形函数，hi为单元节点水头值(m)，n为土的孔隙率；由三角形单元的形函数的特点可知∫∫Nidxdz=A3---(5)]]>式中，A为三角形单元的面积；每延米上全域总势能为E=ΣEe=Σnρg(hie+hje+hme)Ae3---(6)]]>⑷⑷实域总势能及公式推导将渗流场分为湿区、过渡区和干区，若单元中节点水头与坐标的差值均不小于零，则该单元为湿单元；若单元中节点水头与坐标的差值均不大于零，则该单元为干单元；余下的单元则为过渡单元；实域总势能E实是指湿区和过渡区中自由面以下区域的势能总和，即实域总势能E实＝E0‑αE过           (7)式中，E0为全域总势能，E实为实域总势能，E过为过渡单元的总势能，α为过渡单元中虚区面积占单元总面积的比；⑸计算步骤①对土石坝进行有限元划分，首先对坝体进行竖向和横向剖分，得到规则的四边形单元和三角形单元，再对各四边形单元进行同方向对角线连接，剖分为三角形单元，按照由上向下，由左向右的顺序依次编号，或直接生成三角形单元；②将已知上、下游各节点的水头值和假定溢出点的水头值(h＝z)代入式(3)，为编程方便，将已知水头值节点在总渗透矩阵中所对应的主对角线元素变为1，其余为0，将已知水头值代入所对应的水头列向量，无需对渗透矩阵、水头列向量及自由项进行降阶，再依据式(3)计算出各待求节点的水头值；③计算各节点的水头值与纵坐标的差值，将渗流场划分为湿区、过渡区和干区，通过对过渡区单元的计算，依据式(8)计算出过渡区单元的面积比α；④将干区单元丢弃，保持赋予溢出点的位置不变，依据式(3)重新求解各节点水头值，依据式(6)求解出此时的全域总势能，将丢单元前计算出的面积比α代入式(7)中，无需调整自由面，计算出在当前溢出点位置时的实域总势能E实；⑤将假定溢出点位置升高一个节点，重复步骤②、③、④求解出实域总势能Ep，p为执行此步骤的序号，判断Ep与Ep‑1的大小，若EP＞Ep‑1，执行下一步，若EP＜Ep‑1，重复计算步骤⑤；⑥执行步骤序号p‑2、p‑1、p所对应的节点坐标分别为(x1,z1)、(x2,z2)、(x3,z3)，令x2＝x1+(x3‑x1)qβ/(z3‑z1)，z2＝z1+qβ(β为溢出点最小误差限值)，并将(x2,z2)作为新的溢出点坐标，重复步骤②、③、④求解实域总势能Eq，q为执行此步骤的序号；⑦判断Eq与Eq‑1的大小，若Eq＜Eq‑1，令q＝q+1，回到步骤⑥，若Eq＞Eq‑1，此时执行此步骤序号q‑1所对应的节点坐标即为真实溢出点坐标。