1、高面板堆石壩三維真非線性地震反應分析方法及模型試驗驗證基金項目：“九五”國家科技攻關項目(96-221-02-03-02)和國家電力青年基金項目(97Q02和SPQKJ06)作者簡介：趙劍明（1970），男，山東肥城人，博士/博士后，高級工程師，巖土及土工抗震研究*趙劍明 汪聞韶 常亞屏 (中國水利水電科學研究院巖土工程研究所，北京，100044)摘要基于所建立的土石料三維粘彈塑性動力本構模型，并采用新型三維各向異性有厚度薄單元來模擬面板和堆石的接觸面特性，結合考慮孔隙水壓力的消散和擴散的有效應力分析方法，開發了高面板堆石壩地震反應分析的三維真非線性有效應力分析方法。通過典型壩例分析，比較了真

2、非線性分析與等效線性分析結果的差異，表明真非線性方法較真實地反應了結構的地震反應，而且能夠直接計算出壩體的殘余變形，在理論上更為合理。采用模型壩土石料低應力狀態下實測的靜力和動力特性參數，對三維模型壩進行了地震反應分析，并與相應的大型振動臺模型試驗結果進行了對比分析，兩者結果相當一致，驗證了本研究所提出的計算方法的可靠性。關鍵詞：高面板堆石壩；動力本構模型；地震反應；真非線性；模型試驗1 前言混凝土面板堆石壩以其安全性、經濟性、適應性等優點受到國內外壩工界的普遍重視。目前相當多的面板堆石壩位于強震區，使得面板堆石壩抗震問題研究的重要性和迫切性越來越突出。面板堆石壩地震反應的二維分析已有了較多的

3、計算方法和相應程序，但由于面板堆石壩多修建在狹谷之中，具有明顯的三維效應，僅做二維分析是不夠的，因此面板堆石壩地震反應的三維動力分析成為必要；隨著計算技術和計算機的快速發展，也使得三維分析成為可能，所以近期面板堆石壩的三維動力分析引起了國內外學者的重視。地震反應分析方法從基于的本構模型來分可分為兩大類，一類是基于等價粘彈性模型的等效線性分析方法，另一類是基于（粘）彈塑性模型的真非線性分析方法。等價粘彈性模型概念明確，使用方便，在參數的確定和應用方面積累了較豐富的試驗資料和工程經驗，目前應用較多。然而，等價粘彈性模型的缺陷是明顯的：等價粘彈性模型在加荷與卸荷時模量相同，因而不能直接計算土體在周期

4、荷載連續作用下發生的殘余變形；不能考慮應力路徑的影響；不能考慮土的各向異性；較大應變時誤差大，偏于不安全。由此可見，基于等價粘彈性模型的等效線性分析方法得到的地震反應并不是土體真實的地震反應，要想得到土體真實的地震反應，應采用基于（粘）彈塑性模型的真非線性分析方法。2 高面板堆石壩三維真非線性地震反應分析方法2.1 土石料的三維真非線性動力本構模型“八五”期間，汪聞韶院士和李萬紅1建立了一個二維的適用于土石壩的粘彈塑性動力本構模型，初步應用于工程分析，取得了一些較好的結果。本研究根據面板壩的特點，結合試驗資料，對該模型進行推廣和完善，建立了適用于面板壩的三維粘彈塑性動力本構模型，并開發了相應三

5、維真非線性分析方法。本研究采用的粘彈塑性模型將土視為粘彈塑性變形材料，模型由初始加荷曲線、移動的骨干曲線和開放的滯回圈組成。這種真非線性模型的特點是：（1）與等效線性粘彈性模型相比，能夠較好地模擬殘余應變，用于動力分析可以直接計算殘余變形；在動力分析中可以隨時計算切線模量并進行非線性計算，這樣得到的動力響應過程能夠更好地接近實際情況。（2）與基于Masing準則的非線性模型相比，增加了初始加荷曲線，對剪應力比超過屈服剪應力比時的剪應力應變關系的描述較為合理；滯回圈是開放的；考慮了振動次數和初始剪應力比等對變形規律的影響。模型的數學表達式如下：初始加荷曲線：(1)骨干曲線：(2)滯回圈：(3)式

6、(2)和式(3)中，在加荷時取（）、（），在卸荷時取（）、（）。在此非線性動力模型中，骨干曲線和滯回圈的原點不斷移動產生殘余變形，即有：(4)以上諸式中：和為剪應力和剪應變，max為極限剪應力，max=f /Rf ；Rf為破壞比；f為破壞剪應力；' 為有效內摩擦角；' 為有效正應力；0為骨干曲線和滯回圈原點相應的剪應變或稱塑性剪應變；h是以0為零點的剪應變；A、B為模型參數；DRSd為動剪應力比幅值；DRS動剪應力比，DRS=RS-RS0，RS=/'，RS0為初始剪應力比。模型參數A和B以及0可以用剪應力比控制的循環三軸試驗來測定1，主要受振次、動剪應力比幅值和初始剪應

7、力比的影響比較大。模型參數A和B也可由等效線性粘彈性模型參數換算近似得到，換算原則是使兩變形模型的骨干曲線重合和滯回圈包圍的面積相等。2.2 單元型式堆石壩體及地基主要采用三維八結點六面體等參單元來模擬，在邊界不規則處采用六結點五面體三棱柱單元來填充。混凝土面板可用薄板單元來模擬，不過根據已有的研究結果，在二維問題中，面板采用梁單元或四邊形單元來模擬，計算結果相差不大，表明面板隨壩體一起運動時彎曲應力不大；而用板單元模擬面板時存在著板單元與六面體等參元位移模式不協調的誤差。另一方面，盡管面板厚度小，與其它兩個方向的尺寸相差大，但面板單元形態規則，不至于在有限元計算中造成大的誤差，面板用六面體單

8、元離散能夠滿足計算精度要求。基于以上原因，并考慮到計算方便，本研究選用了同類型的三維八結點六面體等參元來模擬。為模擬混凝土面板三維受力狀態下的應力應變關系的非線性，并考慮計算方便，對面板采用分段線彈性模型。2.3 接觸面的模擬本研究中，在文獻2的基礎上，建立了一種三維各向異性有厚度薄單元3來模擬面板和堆石的接觸面特性，如圖1所示。這種單元剛陣與一般的三維六面體等參單元在形式上是相同的。yzx圖1 三維有厚度接觸面薄單元接觸面上的變形可以分為基本變形和破壞變形兩部分。基本變形與其它土體的變形一樣，不管破壞與否都是存在的；破壞變形包括滑動破壞和拉裂破壞，只有當剪應力達到抗剪強度產生了沿接觸面的滑動

9、破壞，或接觸面受拉產生了拉裂破壞時才存在。在三維接觸面中，有三個方向的可能破壞位移：接觸面法向的張裂和沿接觸面的兩個方向的滑移。對于圖1所示的三維有厚度接觸面薄單元，則：(5)2.4 接縫的模擬為了模擬面板周邊縫和垂直縫的特性，在三維有限元計算中設置了無厚度的六面體縫間連接單元4。2.5 動水壓力關于動水壓力，本文采用了附加質量法，具體是用一種廣義邊界元法來處理庫水無界性的問題，給出單位加速度時，不同壩坡條件下的壩面動水壓力分布系數，在計算時將其轉化為相應的附加質量進行分析。2.6 基于真非線性模型的動力反應分析算法鑒于本文采用的真非線性粘彈塑性模型的特點，為了更有效地進行動力反應分析，本文采

10、用了增量法和全量法交替進行的算法以控制增量法的誤差積累。根據非線性粘彈塑性模型及有限元原理，推導出結構的增量和全量方程分別為：（6）（7）式中，和分別為結點位移、速度和加速度，為彈性位移，代表增量，M為質量矩陣，Ct和Cs分別為切線和割線阻尼矩陣，Kt和Ks分別為切線和割線剛度矩陣，Fa為地震力，Fe為應力超過強度時加以修正的等價結點力（超越力）。具體求解按增量步進行。對每一增量步，先求解增量方程（6），然后如果為奇數增量步，則在假定不變的條件下，由全量方程（7）計算彈性位移；如果為偶數增量步，則在假定不變的條件下計算加速度，并用此加速度校正方程（6）中的，以減少用增量法解方程產生的誤差積累。

11、盡管混凝土面板壩中壩體大部分是非飽和的，但對于建造在深厚覆蓋層上的面板壩以及下游水位較高時，還是必須考慮振動孔隙水壓力的影響。因此，本文采用了考慮孔壓消散和擴散的三維非線性有效應力地震反應分析方法5。2.7 真非線性分析與等效線性分析的若干比較為了驗證本文提出的真非線性分析方法的特點，針對一修筑于梯形河谷中壩高100m的典型面板壩，分別采用等效線性分析方法和真非線性分析方法進行了地震反應分析6，在順河向輸入El Centro波，最大加速度取為0.2g。Disp.(cm)等效線性分析真非線性分析圖2 典型結點動位移時程曲線圖3 典型單元動剪應變時程曲線Strain(10-3)等效線性分析真非線性

12、分析圖2是兩種分析方法得到的同一典型結點(698號結點)的動位移時程曲線。圖3是兩種分析方法得到的同一典型單元（715號單元）的動剪應變時程曲線。由圖2和圖3可見，等效線性方法和真非線性方法得到的應變和位移地震反應有著明顯的區別：等效線性分析得出的動應變和位移圍繞零點振動，沒有偏移，無殘余變形產生；真非線性分析得出的動應變和位移在振動過程中偏離零點，產生殘余變形，并且地震過程中殘余變形不斷積累和增長。可見，本文提出的面板壩真非線性動力分析方法和等效非線性方法在概念上有著本質的區別，在計算結果上存在差異，真非線性方法較真實地反應了結構的地震反應，而且能夠直接計算出壩體的殘余變形，在理論上更為合理

13、。3 三維真非線性地震反應分析方法的模型試驗驗證3.1 模型試驗及模型壩計算網格目前，面板堆石壩實際震害和記錄資料很少，特別是缺乏強震區高面板堆石壩的資料。為了研究高面板堆石壩的動力特性和在地震作用下的動力反應性狀，同時為驗證和改進計算方法及程序，中國水科院利用具有世界先進水平的5×5m大型三向模擬地震振動臺，對實際面板壩進行了三維大型振動臺模型試驗6。模型試驗中采用的三維模型壩的幾何比尺為123.5，壩高100cm，壩軸線長315.44cm，上游壩坡坡比為1：1.55，下游壩坡設三級馬道，分別寬2.43cm，壩坡呈上緩下陡，壩坡比分別為1：1.50、1：1.40、1：1.40和1：

14、1.40。在模型試驗中，對模型壩分別輸入了人工合成原波（簡稱人工原波）、壓縮人工合成原波（簡稱壓縮人工波）和壓縮松潘波，地震波特性及波形曲線見文獻6。研究內容包括在不同地震波作用下模型壩的反應加速度、面板應力和應變、壩體變形等地震反應性狀。為了進一步驗證本研究提出的三維真非線性分析方法的可靠性，采用模型土石料低應力狀態下實測的靜力和動力特性參數，具體參數見文獻6，對模型試驗中采用的三維整體模型進行了三維真非線性地震反應分析，并與相應的試驗結果作了比較驗證。在對模型壩進行三維真非線性地震反應分析中，整個模型壩共劃分了1996個結點和1583個單元，包括面板單元、壩體單元、接觸面單元和接縫單元。模

15、型壩橫剖面和縱向單元劃分情況分別如圖4a和圖4b所示。x、y及z分別表示順河向、壩軸向及豎直向。輸入地震波采用了試驗輸入的地震波，即分別為人工原波、壓縮人工波和壓縮松潘波。計算中的輸入形式和方式與試驗中的保持一致，動力分析中的時間步長取為0.002s。100.00155.01152.98面板xz0+69.040+384.48yz圖4 三維模型壩單元剖分圖 (cm)圖4b 縱向（上游立視）(cm)圖4a 橫剖面（0+234.67）3.2 動力分析結果 本研究計算了各種地震波作用下，三維模型壩的動力反應，并與試驗結果進行了比較，取得了較為一致的結論，下面給出其中的典型結果。加速度反應 圖5是三維模

16、型壩在壓縮人工地震波單向（x）作用下，主斷面沿中心線、上游坡和下游坡加速度放大倍數分布情況。圖6是三維模型壩在不同地震波作用下壩頂（中心線）加速度放大倍數與輸入地震加速度最大幅值的關系。圖5 三維模型壩在壓縮人工地震波單向（x）作用下，主斷面沿中心線、上游坡和下游坡加速度放大倍數分布 B1109x1(0.110) 計算結果 B1109x2(0.186) 計算結果 B1109x7(0.554) 計算結果 B1109x1(0.110) 試驗結果 B1109x2(0.186) 試驗結果 B1109x7(0.554) 試驗結果 中心線00.00.51.0Ratioh/H12345上游坡0123450.

17、00.51.0h/HRatio02468100.00.51.0h/HRatio下游坡圖6 三維模型壩中心線上壩頂加速度放大倍數與輸入加速度最大幅值的關系0.00.10.20.30.40.50.60.7Amax(g)0.02.04.06.08.0Ratio 1109x單向輸入 計算結果 1109x單向輸入 試驗結果 SP1109x 單向輸入 計算結果 SP1109x 單向輸入 試驗結果 人工原波x單向輸入 計算結果 人工原波x單向輸入 試驗結果混凝土面板應變和應力反應 圖7是三維模型壩在雙向地震波作用下面板順坡向應變沿壩高分布情況。三維模型壩在單向（x）地震波作用下面板應變與輸入加速度最大值的關

18、系如圖8所示。圖9是三維模型壩在三向地震波作用下面板最大主應力沿壩高的分布情況。圖8 三維模型壩在單向（x）地震波作用下面板應變與輸入加速度最大值的關系Strain (10-4)Amax (g) 0.43H高程處x向應變 計算 0.43H高程處x向應變 試驗 0.70H高程處x向應變 計算 0.70H高程處x向應變 試驗 0.43H高程處y向應變 計算 0.43H高程處y向應變 試驗圖7 三維模型壩在雙向地震波作用下面板順坡向應變沿壩高分布情況Strainxmax (10-4) 1109xy1(0.127, 0.095) 計算 1109xy1(0.127, 0.095) 試驗 1109xy2(

19、0.178, 0.146) 計算 1109xy2(0.178, 0.146) 試驗 1109xz1(0.095, 0.076) 計算 1109xz1(0.095, 0.076) 試驗 1109xz2(0.164, 0.117) 計算 1109xz2(0.164, 0.117) 試驗0.000.050.100.150.200.250.300.00.20.40.60.81.0h/HStress 1max (100kPa) 1109xyz1(0.139, 0.143, 0.146) 計算 1109xyz1(0.139, 0.143, 0.146) 試驗 1109xyz2(0.229, 0.227,

20、0.151) 計算 1109xyz2(0.229, 0.227, 0.151) 試驗 1109xyz2(0.315, 0.298, 0.193) 計算 1109xyz2(0.315, 0.298, 0.193) 試驗0.000.050.100.150.200.250.300.350.400.00.20.40.60.81.0h/H圖9 三維模型壩在三向地震波作用下面板最大主應力沿壩高分布情況壩體殘余變形 本研究計算了三維模型壩在各種地震波作用下的地震殘余變形，其中整理得到的三維模型壩在地震動結束后壩頂中心線沉降變形沿壩軸線的分布情況如圖10所示。圖10 三維模型壩在地震動結束后壩頂中心線沉降變形

21、沿壩軸線的分布0.0040.0080.00120.00160.00200.00240.00280.00320.00B(cm)0.000.200.400.600.801.00S(mm) 1109x2(0.149) 計算 1109x2(0.149) 試驗 1109xz2(0.164, 0.117) 計算 1109xz2(0.164, 0.117) 試驗 1109xyz4(0.315, 0.298, 0.193) 計算 1109xyz4(0.315, 0.298, 0.193) 試驗 b1109x8(0.664) 計算 b1109x8(0.664) 試驗由各圖可見，計算結果與試驗結果相當一致。具體結

22、論見下一小節。3.3對比驗證及小結本文應用本研究提出的三維真非線性地震反應分析方法，采用模型土石料低應力狀態下實測的靜力和動力特性參數，對三維模型壩進行了地震反應分析，并將計算結果對試驗結果進行了對比分析，分析結果表明，計算得到的地震作用下壩體的加速度反應、面板的應力和應變及壩體殘余變形與試驗結果相當一致。從而，在計算和試驗相互印證的同時，也進一步驗證了本研究提出的三維真非線性地震反應分析方法的可靠性，并可以得出如下的一些結論：1）對相同高程，下游坡面上的加速度反應比上游面板上及中心線上的加速度反應大，表現出較明顯的壩體表層加速度放大效應，在工程上應注意加強下游坡的抗震保護，尤其是上部表面；而

23、上游面加速度反應較下游為小，體現了面板的約束作用。2）地震波振幅對加速度放大倍數有較大影響，對于整個壩高，加速度放大倍數有隨輸入加速度最大幅值Agxmax的增大而減小的趨勢，而且使加速度分布也有所變化。3）地震波類型等對加速度放大倍數分布有較大影響，其中人工原波得到的放大倍數顯著地小于經過時間比尺壓縮的人工波所得到的放大倍數，這主要是由于人工原波卓越頻率遠小于壩體基本頻率的緣故；同樣時間比尺Ct的壓縮人工波與壓縮松潘波由于它們卓越頻率段不同，在模型上產生的加速度放大倍數和分布也有所不同。4）計算得到面板順坡向應變及主應變和主應力最大值發生位置雖然有所變化，但最大值位置基本發生在0.50.7H的

24、范圍內。5）地震波類型和加速度幅值對面板應變和應力有重要影響，其影響規律與影響加速度的規律類似。6）壩頂地震沉降在壩軸向中部較大，但總體上地震沉陷量并不大，在本研究輸入的各種地震波作用下均不超過壩高的0.1%，表明面板堆石壩具有較好的抗震性能。4 結論本文基于所建立的土石料三維粘彈塑性動力本構模型，并采用新型三維各向異性有厚度薄單元來模擬面板和堆石的接觸面特性，結合考慮孔隙水壓力的消散和擴散的有效應力分析方法，開發了高面板堆石壩地震反應分析的三維真非線性有效應力分析方法。通過典型壩例分析，比較了真非線性分析與等效線性分析結果的差異，表明真非線性方法較真實地反應了結構的地震反應，而且能夠直接計算

25、出壩體的殘余變形，在理論上更為合理。采用模型壩土石料低應力狀態下實測的靜力和動力特性參數，對三維模型壩進行了地震反應分析，并與相應的大型振動臺模型試驗結果進行了對比分析，兩者結果相當一致，驗證了本研究所提出的計算方法的可靠性。根據模型試驗和計算分析成果，得出了一些有關面板堆石壩地震反應性狀、抗震設計方法及抗震措施方面的研究結論，可為工程建設提供參考。參考文獻1 李萬紅、汪聞韶，無粘性土動力剪應變模型J。水利學報，1993年第9期。2 殷宗澤、朱泓、許國華，土與結構材料接觸面的變形及其數學模型J。巖土工程學報，1994，16（2）。3 趙劍明、汪聞韶，混凝土面板堆石壩面板地震反應分析J。巖石力學

26、與工程學報，2001年，20（s2）。4 顧淦臣，土石壩地震工程M。河海大學出版社，1989。5 趙劍明、汪聞韶、張崇文，土石壩振動孔壓影響因素的研究J。水利學報，2000年第5期。6中國水利水電科學研究院，強震區面板壩大型振動臺模型試驗及動力分析R。“九五”國家科技攻關專題研究報告，2000。3-D Authentic Nonlinear Dynamic Analysis Method of High CFRD and Verification by Model Test DataZhao Jianming Wang Wenshao Chang Yaping(Geotechnical Engineering Department, IWHR, Beijing, 100044)Based on the proposed 3-D authentic nonlinear dynamic constitutive model of soil and applying a new kind of 3-D anisotropic thin element to simulate t