[發明專利]一種基于實域總勢能計算滲流溢出點的方法

在審

• 申請號: 201710556728.6
• 申請日: 2017-07-10
• 公開（公告）號: CN107451399A
• 公開（公告）日: 2017-12-08
• 發明（設計）人: 侯興民;孫偉建
• 申請（專利權）人: 煙臺大學
• 主分類號: G06F19/00
• 分類號: G06F19/00
• 代理公司: 青島發思特專利商標代理有限公司37212
• 代理人: 耿霞

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所對應的節點坐標即為真實溢出點坐標。