劉 杰，荊帥召

(1.重慶市水利電力建筑勘測設(shè)計研究院有限公司，重慶 400020；2.河海大學(xué)水利水電學(xué)院，江蘇 南京 210098)

近些年來，拱壩由于其安全和經(jīng)濟方面的優(yōu)越性，已成為大壩設(shè)計中的三大優(yōu)選壩型之一[1- 3]。隨著我國水利事業(yè)的不斷發(fā)展，在西部地區(qū)已修建了許多超過200m的拱壩，如小灣、溪洛渡、錦屏、白鶴灘等[4- 5]。而壩基巖體的不均勻性作為一種常見的地質(zhì)缺陷，往往會引起拱壩的失穩(wěn)破壞，造成十分嚴(yán)重的后果[6- 7]。因此，研究不均勻壩基條件下壩體的破壞與失穩(wěn)機制對于保障拱壩的安全運行具有重要意義。

有限元法由于在處理復(fù)雜邊界及荷載條件、復(fù)雜結(jié)構(gòu)和非線性問題方面具有獨特的優(yōu)勢，廣泛應(yīng)用于水工結(jié)構(gòu)破壞分析中[8- 11]。而對于水工混凝土這種復(fù)合材料，由于其抗拉強度較低，在模擬混凝土損傷破壞時常會導(dǎo)致結(jié)構(gòu)出現(xiàn)剛度矩陣的不對稱現(xiàn)象，這使得基于隱式算法的有限元分析易出現(xiàn)不收斂的問題。同時，多數(shù)隱式分析需通過反復(fù)調(diào)整模型參數(shù)來避免計算出現(xiàn)不收斂，易引起較大的計算誤差，且試算過程需消耗大量時間。而顯式有限元方法基于動態(tài)算法則無需迭代計算，對于高度非線性問題不存在收斂困難，同時能夠較為真實地模擬結(jié)構(gòu)的實際加載過程[12]。另一方面，通用有限元軟件ABAQUS以其強大的非線性求解能力而被廣泛應(yīng)用于水工結(jié)構(gòu)仿真計算中[13]，其基于顯式算法開發(fā)的顯式分析模塊ABAQUS/Explicit對于求解各類非線性結(jié)構(gòu)力學(xué)問題非常有效[14]。

鑒于此，本文基于三維非線性顯式有限元方法，利用ABAQUS對位于不均勻壩基上的某混凝土拱壩進行數(shù)值模擬分析，總結(jié)出壩基巖體的非均勻性對混凝土拱壩損傷破壞過程和失效模式的影響，可為拱壩結(jié)構(gòu)的整體優(yōu)化設(shè)計和確定壩基巖體工程加固方案提供技術(shù)參考。

1 基于ABAQUS的非線性有限元顯式分析方法

ABAQUS/Explicit中的顯式分析方法采用時間差分法進行積分[15]，通過上一個增量步的動力計算條件進行計算并獲取后一個增量步的動力計算條件，不必再進行平衡迭代，因此計算速度快，一般不會出現(xiàn)計算不收斂問題。具體求解步驟如下：

增量步開始，程序首先求解動力學(xué)平衡方程，節(jié)點合力的求解方程為：

(1)

增量步開始時(t時刻)，計算加速度為：

(2)

采用中心差分法計算加速度對時間的積分，計算速度變化的過程中假設(shè)加速度為常數(shù)，則當(dāng)前增量步中心點速度的計算公式為：

(3)

當(dāng)前增量步結(jié)束時的位移計算公式為：

(4)

不同于隱式分析的無條件穩(wěn)定，顯式分析作為一種條件穩(wěn)定算法，要求時間增量步長Δt不能大于穩(wěn)定時間步長限制值Δtstable，即Δt≤Δtstable，否則會導(dǎo)致結(jié)構(gòu)響應(yīng)出現(xiàn)波動，計算結(jié)果出現(xiàn)無邊界的振蕩發(fā)散。無阻尼時，Δtstable可由下式估計：

(5)

式中，ωmax—模型最高固有頻率；Le—最小尺寸單元的長度；cd—材料波速，可由材料的彈性模量、泊松比和密度決定。

2 不均勻壩基拱壩損傷破壞分析

2.1 工程概況

為研究不均勻壩基條件下混凝土拱壩失效模式，本文以我國西南部某水平拱圈呈拋物線型的混凝土雙曲拱壩作為研究對象進行顯式有限元分析。

該拱壩最大壩高78m，最低建基面高程441m。壩肩地質(zhì)斷面及開挖線如圖1所示，壩基由多個不同厚度的水平巖層組成。根據(jù)地質(zhì)勘查資料，441～480m高程的下部壩基巖層的力學(xué)性能明顯強于480m高程以上的上部壩基巖層，壩基巖體整體處于不均勻狀態(tài)。

圖1 壩肩地質(zhì)斷面圖

2.2 有限元計算模型及荷載工況

本文采用ABAQUS/Explicit計算混凝土拱壩壩體-壩基的三維有限元計算模型。拱壩與壩基系統(tǒng)整體有限元網(wǎng)格主要由八節(jié)點六面體單元組成，結(jié)點和單元數(shù)分別為154813和141877，如圖2(a)所示。為更真實地模擬拱壩壩體損傷破壞的發(fā)生位置和擴展過程，對壩體模型進行了較為精細的網(wǎng)格離散。壩體部分的最小單元尺寸為1.5m，結(jié)點和單元數(shù)分別為24660和21064，圖2(b)所示。

此外，本文采用水密度超載法[16]開展超載工況分析，可間接模擬由于材料強度降低導(dǎo)致大壩失效的過程，具體介紹如下：假定六組超載系數(shù)，超載系數(shù)為施加的水荷載與正常荷載的比值，即將上游水容重增加到正常工況下的1、2、3、4、5、6倍，分別試算拱壩超載1～6倍水壓荷載，分析得到超載系數(shù)下相應(yīng)的應(yīng)力位移分布情況及塑性破壞區(qū)域。

圖2 拱壩三維有限元模型

2.3 材料本構(gòu)模型及力學(xué)參數(shù)

壩體混凝土材料在拉壓過程中因塑性積累和剛度的退化，性能變化極其復(fù)雜。而混凝土連續(xù)損傷塑性模型(CDP)[17]作為一種典型的非線性損傷模型，考慮了混凝土材料拉壓性能的差異，可較好地模擬混凝土材料在外荷載作用下由于損傷引起的剛度退化。因此，為了能較好地描述混凝土材料的力學(xué)特性，本文采用ABABUS內(nèi)置的CDP模型作為壩體混凝土材料的本構(gòu)模型，其他相關(guān)混凝土力學(xué)性能參數(shù)見表1。

表1 混凝土材料力學(xué)性能參數(shù)

對于壩基巖體，本文采用ABAQUS軟件提供的線性Drucker-Prager模型[18- 19]作為本構(gòu)模型，其考慮了中間主應(yīng)力和靜水壓力的影響，可較為準(zhǔn)確地描述壩基巖體的力學(xué)特性。同時，為體現(xiàn)壩基巖體的不均勻性，力學(xué)性能相對較好的壩基巖體采用實際的力學(xué)參數(shù)，而對于力學(xué)性能較差的壩基巖體則采用參數(shù)折減法描述其力學(xué)性能。上部軟弱壩基巖體力學(xué)參數(shù)為折減系數(shù)與下部壩基巖體力學(xué)參數(shù)的乘積，其中軟弱壩基巖體密度與泊松比的值保持不變。折減系數(shù)值R分別取為1.0、0.7、0.5、0.3，不同R值對應(yīng)的軟弱壩基巖體的力學(xué)參數(shù)見表2。

表2 壩基巖體力學(xué)性能參數(shù)

2.4 結(jié)構(gòu)準(zhǔn)靜態(tài)模擬分析

為解決顯式算法分析的條件穩(wěn)定性問題，需對拱壩壩體-壩基系統(tǒng)模型進行準(zhǔn)靜態(tài)模擬分析，將結(jié)構(gòu)響應(yīng)引起的波動性控制在工程可接收的范圍內(nèi)。而滿足顯式計算準(zhǔn)靜態(tài)分析要求的主要因素和要點包括加載幅值曲線、加載時長以及模型的網(wǎng)格劃分3個方面。

由于顯式分析是基于自然時間的求解過程，且由式(5)可知時間增量Δt數(shù)值上非常小，因此模型加載時長若取值過小會引起動態(tài)效應(yīng)造成計算誤差，過大則會提高計算的時間成本。而另一方面，在結(jié)構(gòu)準(zhǔn)靜態(tài)加載過程中，最小自振周期對結(jié)構(gòu)的響應(yīng)具有控制性作用。因此，可將模型加載時長增加到系統(tǒng)最小自振周期的10倍左右[20]，以減小荷載施加過快引起的動態(tài)效應(yīng)，使計算結(jié)果的準(zhǔn)確度和計算效率滿足要求。經(jīng)過計算，拱壩壩體-壩基系統(tǒng)結(jié)構(gòu)的最小自振周期為0.796s，綜合考慮，取加載時長為20s。

其次，在模型加載過程中，若加載速度出現(xiàn)突變，則會引起結(jié)構(gòu)的震蕩，造成計算結(jié)果出現(xiàn)較大偏差[21]。因此，為使準(zhǔn)靜態(tài)分析結(jié)果的誤差更小、計算效率更高，本文ABAQUS中選用如圖3所示加載過程平滑、波動性較小的光滑函數(shù)曲線。該曲線將兩個自定義的幅值間以五階多項式進行過渡，在自定義的幅值點處速度和加速度為0。

此外，由于最小網(wǎng)格尺寸與數(shù)值模擬的精度成反比，且由式(5)可知，穩(wěn)定極限大致與最短的單元尺寸成比例，因此最小單元尺寸的合理選取對顯式計算結(jié)果影響重大，且結(jié)構(gòu)網(wǎng)格的劃分應(yīng)盡量均勻。本文限于篇幅和計算量，只選取一種較為合理的網(wǎng)格劃分方式。

由于實現(xiàn)準(zhǔn)靜態(tài)分析的核心要求是選擇合理的加載速度，為判斷上述方法是否滿足準(zhǔn)靜態(tài)分析要求，本文以加載過程中結(jié)構(gòu)的動能與內(nèi)能的比值不超過10%作為判斷準(zhǔn)則[22]。其中，動能表征結(jié)構(gòu)運動速度，內(nèi)能表征結(jié)構(gòu)變形程度。

圖3 平滑幅值曲線

3 結(jié)果分析

3.1 準(zhǔn)靜態(tài)分析驗證

由2.4節(jié)可知，在ABAQUS后處理中可繪制結(jié)構(gòu)動能與內(nèi)能的比值隨時間變化的曲線用于檢驗結(jié)構(gòu)在運算過程中是否處于準(zhǔn)靜態(tài)。本次分析結(jié)果如圖4所示，結(jié)構(gòu)動能與內(nèi)能的比值在整個加載過程中的最大值僅為0.0113%，滿足準(zhǔn)靜態(tài)分析要求。

圖4 動能/內(nèi)能隨時間變化值

3.2 超載破壞分析結(jié)果

超載分析中，本文以塑性區(qū)貫通作為拱壩結(jié)構(gòu)完全破壞的標(biāo)志，以等效塑性應(yīng)變超過100με作為混凝土材料進入塑性的指標(biāo)[23]。針對不同R值下的壩基條件，圖5—8分別展示了塑性區(qū)貫通時相應(yīng)超載倍數(shù)下拱壩的塑性區(qū)(紅色區(qū)域)分布。

(1)當(dāng)R=1.0時，壩基巖體可視為均勻巖體，超載倍數(shù)達到6時壩體塑性區(qū)貫通，如圖5所示，拱壩的塑性區(qū)主要集中在拱壩壩體底部、上游面拱端處、下游面上部拱端和靠近壩頂?shù)闹行膮^(qū)域。

圖5 折減系數(shù)R=1.0，超載倍數(shù)6.0

(2)當(dāng)R=0.7時，超載倍數(shù)達到4.5時壩體塑性破壞區(qū)貫通，如圖6所示，且與R=1.0時拱壩上下游面塑性破壞區(qū)分布位置基本一致。

圖6 折減系數(shù)R=0.7，超載倍數(shù)4.5

(3)當(dāng)R=0.5時，拱壩在超載倍數(shù)為3.5時塑性區(qū)貫通，如圖7所示，相比于R≥0.7時，上游面塑性區(qū)在480m高程以上分布有所增加，而其他塑性區(qū)分布位置基本保持不變；下游面塑性破壞區(qū)在480m高程處從拱端向中心區(qū)域發(fā)展，但仍未貫通，而靠近壩頂中心區(qū)域的塑性區(qū)減少。

圖7 折減系數(shù)R=0.5，超載倍數(shù)3.5

(4)當(dāng)R=0.3時，拱壩在超載倍數(shù)達到2.5倍后塑性破壞區(qū)貫通，如圖8所示，上下游塑性破壞區(qū)皆貫通于480m高程拱圈附近區(qū)域。上游面底部由于中心區(qū)域的屈服破壞使得梁的作用失效，并未出現(xiàn)塑性屈服。

圖8 折減系數(shù)R=0.3，超載倍數(shù)2.5

3.3 失穩(wěn)模式分析

(1)失效模式Ⅰ：當(dāng)R≥0.7時，壩基巖體不均勻程度較低，塑性區(qū)分布位置隨R值的改變并無顯著變化，這意味著下部壩基巖體對拱壩的失效模式影響更大。此種損傷破壞模式可定義為失效模式Ⅰ。在失效模式Ⅰ下，拱壩的損傷破壞主要受到壩體與下部壩基之間相互作用的影響。

(2)失效模式Ⅱ：相較于R≥0.7的情況，R=0.5時，在480m高程拱端附近也出現(xiàn)了塑性區(qū)，此時壩基巖體的不均勻程度對拱壩破壞機制產(chǎn)生了一定影響，此種損傷破壞模式可定義為失效模式Ⅱ。對于失效模式Ⅱ，拱壩的失穩(wěn)破壞主要受到壩體與上部和下部壩基共同作用的影響。

(3)失效模式Ⅲ：當(dāng)R=0.3時，壩基巖體不均勻程度較高，拱壩的塑性區(qū)主要集中在480m高程處對應(yīng)的拱圈附近區(qū)域。相較于失效模式Ⅰ，拱壩的破壞位置和形式發(fā)生了明顯變化，此種損傷破壞模式可定義為失效模式Ⅲ。對于失效模式Ⅲ，拱壩的失穩(wěn)破壞主要受到壩體與上部壩基之間相互作用的影響。

4 結(jié)語

為研究不均勻壩基條件下混凝土拱壩的失效模式，基于實際工程對位于不均勻壩基上的混凝土拱壩開展了一系列的數(shù)值計算，主要結(jié)論如下：

(1)基于顯示有限元的準(zhǔn)靜態(tài)分析方法得到混凝土拱壩失效的全過程。

(2)受壩基不均勻程度的影響，混凝土拱壩存在3種失效模式。

(3)對于不同的失效模式，提出了當(dāng)折減系數(shù)值R達到一定數(shù)值后，一般性的工程措施很難滿足大壩安全需要，其余可根據(jù)混凝土的失效過程對壩基采用相對應(yīng)的加固措施，可為后續(xù)類似工程設(shè)計提供技術(shù)參考，具有較為重要的學(xué)術(shù)意義與工程應(yīng)用價值。