使用動態(tài)規(guī)劃結合有限元應力場對濕陷性土壩進行瞬態(tài)穩(wěn)定性分析

鄧倫宇

(南昌市防洪排澇工程事務中心，南昌，330009)

0 引言

小型濕陷性土壩常見于我國半干旱地區(qū)[1]。這些水壩的建造是為了提供一個水庫，可以在短暫的雨季積累降水，在漫長的旱季為農業(yè)和公共用水提供用水[2]。由于缺乏合適的壓實設備和水，這些大壩中存在一些在干燥的最佳亞穩(wěn)態(tài)條件下壓實很差的情況[3]。因此，在第一次蓄水過程中，壓實土在濕潤狀態(tài)下會發(fā)生較大的坍塌和抗剪強度損失。破壞機制通常以上游邊坡的楔狀破壞為特征，并結合下游邊坡的管道破壞。

本文對某典型小型濕陷性土壩進行了穩(wěn)定性分析。為了進行分析，將動態(tài)規(guī)劃法(DPM)代碼SAFE-DP與全耦合流體力學有限元代碼COUPSO相結合。假設使用的本構參數(shù)從狀態(tài)表面和非線性完美塑性彈性行為獲得，對土壩的上游和下游的斜坡進行了研究。本文介紹了對滑移面形狀和位置的演變以及施工結束、快速蓄水池充填和充填后階段安全因素演變的分析。此外，在文中討論了DPM的性能。

1 動態(tài)規(guī)劃法

動態(tài)規(guī)劃法(DPM)是線性泛函最大化和最小化的通用方法[4-5]。Baker(1980)介紹了DPM在邊坡穩(wěn)定性分析中的應用。Baker(1980)提出的程序使用了DPM，同時保留了關于層間力的假設。后來，幾位研究人員提出了Baker程序的擴展，使用了有限元應力場。

1.1 優(yōu)化過程

DPM是一個通用的優(yōu)化程序，可用于確定故障機制和相應的最小安全系數(shù)[6]。如圖1所示，滑移面由一系列由階段和狀態(tài)點定義的線性段形成。安全系數(shù)Fs可以以離散形式定義，如下所示：

(1)

圖1 用動態(tài)規(guī)劃進行穩(wěn)定性分析的解析方案

其中，n為段數(shù)；τfi為沿i段土體的抗剪強度；△Li為i段的長度；τi為沿段i作用剪應力。由于該方法不適用于非加性泛函，所以不能用DPM最小化泛函Fs。Fs的最小化可以通過最小化以下輔助泛函來實現(xiàn)：

(2)

式中，G的最小化是利用最優(yōu)原則來實現(xiàn)的。根據(jù)最優(yōu)性原則，后一階段最優(yōu)函數(shù)Hi+1(j)是前一階段最優(yōu)函數(shù)Hi(k)的函數(shù)，如下所示：

Hi+1(j)=min[Hi(k)+DGi+1(j,k)]

(3)

式中，DGi+1(j，k)是連續(xù)階段在兩個狀態(tài)點之間通過的“成本”。初始階段的最優(yōu)函數(shù)值(H1(k)=1)在任何狀態(tài)點都等于零。在最后階段，最優(yōu)函數(shù)對應于G的最小值，如下所示：

(4)

其中，Ri為沿線段i的阻力；Si為沿線段i的剪切力。一旦確定了最后階段的最優(yōu)函數(shù)，通過連接從最后階段追溯到最后階段的最佳狀態(tài)點，找到最佳路徑。臨界滑動面對應于最佳路徑。

由式(1)可以看出，優(yōu)化過程是非線性的。對于給定的容錯性，必須假定F的初始值，并且必須重復計算，直到F的值收斂。

圖2 連接兩個狀態(tài)點的線段上的最佳變量

1.2 計算在分析中使用的應力

可根據(jù)抗剪強度包絡線并使用應力狀態(tài)變量定義沿任何段“i”的抗剪強度?？梢允褂糜邢拊椒ǐ@得兩個正交平面在x和y方向(σx、σy和τxy)上的應力狀態(tài)場。根據(jù)σx、σy和τxy定義的應力狀態(tài)場，可使用以下方程式計算作用于以θ角傾斜的任何管段上的法向應力σθ和剪切應力τθ(見圖2)：

σθ=σxsin2θ+σycos2θ-τxysin2θ

(5)

(6)

2 控制水力學土壤行為的耦合微分方程

用于濕陷性壩穩(wěn)定性分析的應力狀態(tài)場已使用有限元程序COUPSO計算。該程序COUPSO以完全耦合的方式求解二維平面應變條件的平衡方程和水守恒方程。分析非飽和土的行為，飽和條件是一種特殊條件。本構模型和參數(shù)基于狀態(tài)表面概念和非線性彈性、理想塑性行為。使用應力狀態(tài)變量定義孔隙比、飽和度、水力傳導率和泊松比的狀態(tài)表面以及剪切強度包絡線；即凈總應力(σ-ua)和基質吸力(ua-uw)。基于空氣相是連續(xù)的并與大氣接觸的假設，空氣流動被忽略。

2.1 靜力平衡

通過考慮作用在代表性土壤元素體積(REV)上的力的平衡，可以獲得控制靜態(tài)平衡的偏微分方程。微分方程中的應力可以使用應力-應變關系(廣義胡克定律)替換為應變。分別針對x和y方向獲得以下等式：

(7)

(8)

2.2 水質量守恒

考慮水在土壤中的進出流動，利用水流的達西定律，可以得到水的二維流動和守恒的偏微分方程。得到的方程如下：

(9)

其中，βw1=(m1w/m1s)；βw2=m2w-(m1wm2s/m1s)；m1w=[2(1+μ)/Ew]；m2w=[1/Hw-(E/H)/Ew]；m1s=2(1+μ)(1-2μ)/E；m2s=[1/Hx+1/Hy+2μ/Hz]；Ew為相對于凈平均應力變化的水相體積模量，σ-ua；Hw為相對于基質吸力變化的水相體積模量；εv為總體積應變；為發(fā)散算子；k為導水率；γw為水的容重；y為高度。

3 小型亞穩(wěn)態(tài)大壩穩(wěn)定性分析

所分析的假設土壩代表了我國半干旱地區(qū)遇到的典型小型亞穩(wěn)壩。以下章節(jié)將介紹分析中使用的土壤特性、問題幾何、初始和邊界條件以及穩(wěn)定性分析的結果。本文給出的結果將集中于使用SAFE-DP程序和獲得的應力狀態(tài)有限元場獲得的穩(wěn)定性分析結果。

3.1 土壤性質

本研究使用三軸儀、測力儀和直剪儀，在有無基質吸力控制的情況下，對土壤性質進行了實驗測定。

此處總結的土壤特性基于在單調潤濕路徑下進行的測試，基質吸力從370到0kPa分階段變化。已有學者提出了適合實驗結果的數(shù)學函數(shù)，這些函數(shù)基于應力狀態(tài)變量(σ-ua)和(ua-uw)。所采用的數(shù)學函數(shù)能夠擬合獲得的實驗數(shù)據(jù)，是連續(xù)的，并且具有連續(xù)的導數(shù)。用于定義空隙率e和飽和度S狀態(tài)曲面的函數(shù)如下：

(10)

(11)

其中，eu=0.7697-0.0073ln(σm-ua)，為崩前階段的孔隙比；ef=0.752+{[-0.142/(1+(σm-ua)/75)-3.5]}，為塌陷后階段的孔隙比；b=39.01(σm-ua)-0.6103，是控制坍塌相形狀的參數(shù)；c=0.00094(σm-ua)2+0.07465(σm-ua)+11，是塌陷中期基質吸力的對應參數(shù)；(σm-ua)為凈正常平均應力；S0=0.375，為初始飽和度；d=0.9769，為參數(shù)，定義了塌陷階段潤濕曲線的斜率；cm=20，是塌陷階段中期基質吸力對應的參數(shù)。

水導率的實驗結果采用以下函數(shù)擬合：

(12)

其中kp=-1.4×10-7+6.26×10-8ln(σm-ua)，表示總平均應力變化時導率的變化；ψcr=3.0，為進氣量；λ=2.1，為經驗常數(shù)。

剪切強度包絡線利用具有基質吸力控制的直接剪切儀獲得。由于較高的吸力值的貢獻被發(fā)現(xiàn)是可以忽略的，所以基質吸力值在100kPa到0kPa之間變化。用以下函數(shù)表示剪切強度包絡線：

τff=a1+b1(σ-ua)+c1(ua-uw)
+d1(σ-uu)(uu-uw)p

(13)

其中，a1=5.0，是有效內聚截距；b1=0.1944，為有效摩擦角的正切；并且c1=0.3238、d1=0.09319和p=4.307×10-7是擬合參數(shù)。

3.2 幾何問題、邊界和初始條件

圖3顯示了用于分析施工、水庫快速填充和后填充階段的幾何問題和有限元網(wǎng)格。施工階段使用5層進行模擬，并考慮線性彈性行為，μ=0.3和E=5800kPa。施工階段唯一的作用荷載是與土壤單位重量相對應的體荷載(γnat=14.75kN/m3)。下邊界的邊界條件對應于剛性基礎，沒有相對滑動(x和y方向沒有位移)。施工階段的孔隙水壓力變化被忽略。

假設370kPa的基質吸力和施工階段模擬獲得的應力為快速水庫蓄水階段模擬的初始條件。通過零檢驗獲得初始基質吸力值后，采用相同的基質吸力值和從快速水庫蓄水分析中獲得的應力作為蓄水后階段模擬的初始條件。填充后分析假設的水流邊界條件如圖3所示。在快速蓄水過程中，不允許出現(xiàn)水滲透和孔隙水壓力變化。時間離散包括17d的初始時間步，然后是0.2d的較小時間步。為了防止數(shù)值不穩(wěn)定，需要較大的初始時間步長。

圖3 大壩橫截面、有限元網(wǎng)格和邊界條件

3.3 結果與討論

圖4和圖5分別為上、下游邊坡臨界滑面形態(tài)和位置的演化。可以看出，施工階段結束時，上下游邊坡的臨界滑移面和安全系數(shù)是一致的，且與預期一致。

水庫快速填筑階段后上游邊坡的臨界滑動面略深于施工階段結束時的臨界滑動面。這是由于水庫荷載導致上游邊坡附近的圍壓較高。水庫荷載的另一個影響是，由于摩擦產生的較高剪切強度，安全系數(shù)顯著增加?？紤]到水庫荷載引起的圍壓增加隨上游邊坡距離的增加而減小，施工結束和快速填筑階段獲得的相關結果是正確的。水庫快速蓄水階段后，下游邊坡的臨界滑動面和安全系數(shù)沒有發(fā)生任何顯著變化。

圖4和圖5還顯示了充填后階段的臨界滑動面。濕潤鋒提前導致大壩不均勻沉降。因此，應力分布會發(fā)生變化，以適應此類沉降。水庫蓄水后58d內，上游臨界滑動面相對較淺。這一結果與濕潤鋒的位置以及該區(qū)域剪切強度的相應降低一致。在此階段，下游邊坡的臨界滑動面沒有顯著變化。

水庫蓄水后98d的上游臨界滑動面比前一個滑動面(58d)更深。滑動面的形狀與圓形相差很大。臨界滑動面的加深和安全系數(shù)的降低是由于濕潤鋒的進一步推進。下游斜坡再次受到濕潤鋒的顯著影響，濕潤鋒仍然局限于上游斜坡附近的區(qū)域。

圖4 上游邊坡穩(wěn)定性分析結果

圖5 下游邊坡穩(wěn)定性分析結果

水庫蓄水后138d內獲得的上游臨界滑動面對應于大壩理論失效的階段(Fs<1)。隨著濕潤鋒的推進，獲得的臨界滑動面比之前的表面更深?；瑒用嫘螤蠲黠@偏離圓形，由剛性壩基礎的位置控制。未觀察到下游邊坡安全系數(shù)和臨界滑動面的顯著變化。由于不均勻沉降引起的應力重新分布，安全系數(shù)幾乎沒有降低。圖6總結了上述安全演變因素。

圖6 下游和上游邊坡安全系數(shù)隨時間的演化

4 結論

本文采用動態(tài)規(guī)劃法(DPM)結合有限元應力狀態(tài)場分析了濕陷性土壩的穩(wěn)定性。在施工結束、快速水庫蓄水和蓄水后階段進行穩(wěn)定性分析。利用計算機程序COUPSO，通過應力和水流的完全耦合分析，獲得了應力狀態(tài)場。DPM和耦合分析的結合允許放松滑動面形狀限制，并對土壤水力力學行為進行相對嚴格的建模。

盡管瞬態(tài)潤濕過程中亞穩(wěn)大壩的應力場十分復雜，但動態(tài)規(guī)劃程序SAFE-DP能夠識別現(xiàn)場觀測到的破裂機制。安全系數(shù)結果也與期望值定量一致。濕陷性大壩在建壩和快速充庫階段是穩(wěn)定的，濕潤鋒通過壩后變得不穩(wěn)定。結果表明：蓄水后138d，水庫安全系數(shù)達到1以下。