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

权利要求书

1.一种基于实域总势能计算渗流溢出点的方法，其特征在于，包括如下步骤：

⑴支配方程和边界条件

根据达西定律和地下水运动的连续性条件，不考虑土和水的压缩性，二维均质各向异性土体的稳定渗流满足偏微分方程

∂∂x(kx∂h(x,z)∂x)+∂∂z(kz∂h(x,z)∂z)=0---(1)]]>

式中，h(x,z)是水头函数，k_x、k_z分别为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)]]>

式中，E_e为每延米上单个三角形单元的势能(J/m)，ρ为水的密度(kg/m³)，g为重力加速度(m/s²)，N_i为三角形单元的形函数，h_i为单元节点水头值(m)，n为土的孔隙率；

由三角形单元的形函数的特点可知

∫∫Nidxdz=A3---(5)]]>

式中，A为三角形单元的面积；

每延米上全域总势能为

E=ΣEe=Σnρg(hie+hje+hme)Ae3---(6)]]> 1

⑷

⑷实域总势能及公式推导

将渗流场分为湿区、过渡区和干区，若单元中节点水头与坐标的差值均不小于零，则该单元为湿单元；若单元中节点水头与坐标的差值均不大于零，则该单元为干单元；余下的单元则为过渡单元；实域总势能E_实是指湿区和过渡区中自由面以下区域的势能总和，即实域总势能

E_实＝E₀-αE_过 (7)

式中，E₀为全域总势能，E_实为实域总势能，E_过为过渡单元的总势能，α为过渡单元中虚区面积占单元总面积的比；

⑸计算步骤

①对土石坝进行有限元划分，首先对坝体进行竖向和横向剖分，得到规则的四边形单元和三角形单元，再对各四边形单元进行同方向对角线连接，剖分为三角形单元，按照由上向下，由左向右的顺序依次编号，或直接生成三角形单元；

②将已知上、下游各节点的水头值和假定溢出点的水头值(h＝z)代入式(3)，为编程方便，将已知水头值节点在总渗透矩阵中所对应的主对角线元素变为1，其余为0，将已知水头值代入所对应的水头列向量，无需对渗透矩阵、水头列向量及自由项进行降阶，再依据式(3)计算出各待求节点的水头值；

③计算各节点的水头值与纵坐标的差值，将渗流场划分为湿区、过渡区和干区，通过对过渡区单元的计算，依据式(8)计算出过渡区单元的面积比α；

④将干区单元丢弃，保持赋予溢出点的位置不变，依据式(3)重新求解各节点水头值，依据式(6)求解出此时的全域总势能，将丢单元前计算出的面积比α代入式(7)中，无需调整自由面，计算出在当前溢出点位置时的实域总势能E_实；

⑤将假定溢出点位置升高一个节点，重复步骤②、③、④求解出实域总势能E_p，p为执行此步骤的序号，判断E_p与E_p-1的大小，若E_P＞E_p-1，执行下一步，若E_P＜E_p-1，重复计算步骤⑤；

⑥执行步骤序号p-2、p-1、p所对应的节点坐标分别为(x₁,z₁)、(x₂,z₂)、(x₃,z₃)，令x₂＝x₁+(x₃-x₁)qβ/(z₃-z₁)，z₂＝z₁+qβ(β为溢出点最小误差限值)，并将(x₂,z₂)作为新的溢出点坐标，重复步骤②、③、④求解实域总势能E_q，q为执行此步骤的序号；

⑦判断E_q与E_q-1的大小，若E_q＜E_q-1，令q＝q+1，回到步骤⑥，若E_q＞E_q-1，此时执行此步骤序号q-1所对应的节点坐标即为真实溢出点坐标。