非線性比例邊界有限元在面板壩分析中的應用

鄒德高，陳楷，劉鎖，周揚

(1.大連理工大學 水利工程學院；海岸和近海工程國家重點實驗室，遼寧 大連 116024；2.中國電建集團 華東勘測設計研究院有限公司，杭州 311122)

面板堆石壩由于可就地取材、經濟、抗震性強等優點，一度成為壩工界青睞的首選壩型[1]，但其結構尺寸相差懸殊，面板最小厚度僅為0.3 m，與壩體整體尺寸差別數百倍至千倍之多。面板作為大壩最重要的結構，應力梯度高，應布置較密的網格[2]，且實際研究[3-4]表明，精細的網格密度將獲得更合理的模擬結果。

有限單元法使用簡便、通用性強、工程應用廣[5]，是面板壩等巖土工程結構安全評價的有效技術手段，但該方法單元形式單一，平面網格嚴格限制為三角形和四邊形單元，網格離散限制多，對復雜邊界適應性較差[6]。尤其在尺度跨越大的面板壩工程精細化分析中，需消耗大量時間進行較為繁瑣的人機交互前處理，一定程度上降低了自動化程度和分析效率[7-8]。近年來，Wolf等[9]提出了彈性的比例邊界多邊形有限元(scaled boundary polygon finite element method, SBPFEM)，該方法更加靈活自由，融合了邊界元(boundary element method, BEM)和有限元法(finite element method, FEM)的優點，并規避了兩者的不足，且已被證明精度和收斂性都高于多邊形有限元的數值算法[10]，被學者們應用到多個領域：如電磁問題分析[11]、土與結構相互作用分析[12]、裂紋擴展分析[13-14]、面板動水壓力研究[15-16]、三維復雜多面體應用[17-18]、波的傳播問題分析[19]等。此外，為探索SBPFEM的非線性應用，Chen等[20]利用半解析的彈性解構造單元形函數，采用兩套高斯點方案，實現了比例邊界多邊形有限元的非線性分析，并將其拓展于多孔介質動力液化分析[21]及三維彈塑性巖土工程應用[22-24]。

采用筆者自主開發的非線性比例邊界多邊形單元，對典型面板堆石壩進行了系統的靜力、動力和永久變形分析計算，并與傳統有限單元法對比，以驗證提出的方法在面板壩工程分析中的可靠性和合理性。隨后，與高效的四分樹離散技術無縫耦合，進行面板壩結構快速跨尺度精細化分析應用，以揭示SBPFEM能在面板壩工程精細化分析中脫穎而出的顯著優勢。

1 工程應用探究

計算程序采用作者完全自主研發的大型巖土工程非線性分析軟件平臺GEODYNA[25]，非線性比例邊界多邊形單元被集成到該平臺。

1.1 理論概要

比例邊界有限元是一種彈性的半解析數值方法，其環向離散、徑向解析的特點使其可降低一個計算維度，減少一定的計算量。對任意問題域，可用任意n邊形對求解域離散(n> 2)，如圖1所示。通過對每個子域進行求解可獲得整個求解域的數值結果。理論詳盡推導，可參見文獻[9-10,20]。

圖1 比例邊界多邊形有限單元Fig.1 Scaled boundary polygon finite

彈性框架限制了該方法的應用領域，本著拓展推廣SBFEM的原則，采用非線性化的比例邊界多邊形單元。采用兩套高斯點策略：保留原始的邊界線高斯點，用于計算相關系數矩陣和半解析的彈性形函數；同時，引入內部高斯積分點，可積分獲得考慮了材料非線性的單元剛度矩陣，進而組裝得到計算域的總體剛度矩陣，并采用Newton-Raphson迭代算法用于求解非線性平衡方程，從而實現該方法的非線性分析。開發的單元可以求解常規單元及傳統常規有限元無法計算的多邊形單元(常規凸多邊形和四分樹單元)，見圖2。詳細理論推導和計算可參考文獻[9,20,26-29]。

圖2 比例邊界有限元支持單元類型Fig.2 Supported element types in

1.2 程序驗證

采用彈性懸臂梁結構端部受彎分析驗證實現程序的正確性和計算精度，結構尺寸及網格如圖3所示，其中l=1.2 m、h=0.4 m、彈性模量E=3×1010、泊松比v=0.2、荷載F=2.4×107N。

圖3 懸臂梁結構及網格Fig.3 Cantilever structure and

采用SBPFEM和FEM分別計算相同結構，并與理論近似解(采用文獻[20]提供的公式)進行對比，分析結果列于表1。從表1可以看出，采用的SBPFEM具有更高的精度。

表1 計算結果與理論解比較Table 1 Comparison between calculation results and theoretical approximate solutions

1.3 工程算例驗證

以典型面板堆石壩為例，壩高為155 m，壩頂寬度為17 m，面板坡比1∶1.8，主體結構材料分區包括面板、墊層區以及堆石料區。為更合理地模擬面板變形和應力，面板與墊層之間添加無厚度的Goodman接觸面單元，面板與趾板之間添加縫單元。常規有限元網格離散單元個數為1 702，節點個數1 816，分析模型見圖4。

圖4 面板壩有限元網格Fig.4 Mesh of concrete face rock-fill

1.4 材料參數

靜、動力均采用堆石料廣義塑性模型[30]模擬筑壩材料變形特性，參數列于表2。面板與墊層間設置的接觸面采用廣義塑性接觸面本構模型[31]，詳細參數見表3。

表2 筑壩材料廣義塑性模型參數Table 2 Generalized plastic model parameters of dam material in static analysis

表3 廣義塑性接觸面參數Table 3 Parameters of the generalized plastic interface model

面板、趾板及基巖采用線彈性模型，參數見表4。縫單元參數采用文獻[32]建議值，其法向壓縮剛度為25 GPa/m，法向拉伸剛度為5 MPa/m，切向剛度為1 MPa/m。靜力計算考慮了壩體的施工填筑和蓄水過程，壩體填筑分22個荷載步完成，隨后分30個荷載步蓄水至152 m高程。

表4 線彈性材料參數Table 4 Parameters of linear model for concrete

1.5 地震動輸入和計算方案

1.5.1 地震動輸入 地震動輸入采用規范譜人工波，順河向峰值加速度為0.2g，豎向峰值加速度為順河向的2/3。地震波加速度時程見圖5。地震波持續時長為25.00 s，計算時間步間隔Δt取0.02 s。

圖5 地震波加速度時程曲線Fig.5 Time history curve of seismic wave

1.6 計算結果及分析

1.6.1 靜力結果 圖6為滿蓄期的壩體位移等值線分布圖。從圖6可以看出，FEM計算結果與NSBPFEM結果吻合較好。圖7為滿蓄期面板應力變形，FEM與NSBPFEM結果基本一致。

表5給出了滿蓄時壩體和面板應力變形的極值。從表5可見，兩種方法計算結果相差較小，證明NSBPEM在模擬大壩填筑蓄水過程中的精度可以得到保證。

圖6 滿蓄期壩體位移(單位：m)Fig.6 Displacement of the dam after impoundment

圖7 滿蓄期面板應力變形Fig.7 Deflection and stress along slope direction of face

方法壩體位移/m上游下游豎向面板應力變形撓度/m應力/MPaFEM0.0650.1250.8120.2189.24NSBPEM0.0650.1260.8160.2189.15相差/%0.0000.7900.4900.0000.98

1.6.2 動力結果 圖8為地震過程中面板順坡向應力隨高程分布規律。FEM計算應力與NSBPFEM計算結果非常接近。

表6列出了兩種方法動力計算中壩頂A點的加速度的極值情況，結果顯示兩者相差不大，說明了NSBPFEM在動力計算中的可靠性。

圖8 地震作用下面板應力(壓為負)Fig.8 Stress along slope direction of concrete face during

方法順河向豎向FEM1.971.99NSBPEM1.982.01相差%0.510.99

1.6.3 永久變形結果 圖9為地震后壩體永久變形計算結果。從圖9可以看出，兩者變形輪廓基本完全一致，說明NSBPFEM可用于震后永久變形計算分析。

圖9 堆石體震后永久變形(放大10倍)Fig.9 Permanent deformation of the dam body after

2 高效的跨尺度精細化分析

除可用于常規單元(三角形和四邊形)分析外，SBPFEM的優勢在于其可與傳統有限單元不能直接處理的四分樹網格離散技術無縫結合，進行高效的跨尺度精細化分析。四分樹可快速實現尺度跨越，從1 m級到1 mm級僅需10層遞歸劃分(210= 1 024)[33]。可自動進行高質量單元離散，且劃分網格以正方形為主，僅邊界和過渡區為多節點單元。

2.1 計算模型與參數

采用四分樹技術，并結合SBPFEM對面板壩進行高效跨尺度精細化分析應用，其中，面板單元尺寸設定為0.5 m，壩體單元大小為8.0 m。如圖10～圖11所示，整個模型剖分時間僅為0.216 s。該技術可高效地實現跨尺度精細化分析，避免了文獻[7-8]的不足。值得注意的是，該方法進行網格二次劃分時，僅需重新定義最大和最小網格尺寸即可，顯著改善了網格離散效率。

圖10 面板壩四分樹網格Fig.10 Quadtree mesh of concrete face

圖11 面板壩四分樹網格局部Fig.11 Partial quadtree mesh of concrete face

靜力計算分25個荷載步模擬實際施工過程，分30步緩慢蓄水至152 m高程，靜力計算所得的應力變形作為初始應力輸入，進行動力時程分析，相應的計算參數列于表2～表4。地震動激勵輸入如圖5所示。

2.2 計算結果分析

圖12 滿蓄期面板應力變形Fig.12 Deflection and stress along slope direction of

地震動作用下，面板動應力分布情況繪制于圖13。從圖13可看出：面板應力有少量差別，其中，SBPFEM計算的最大順坡向應力最大值為13.2 MPa，FEM計算的最大值則為12.2 MPa。最小順坡向應力計算中，SBPFEM最大壓應力為6.91 MPa，最大拉應力為2.21 MPa；FEM計算的最大壓應力為5.92 MPa，最大拉應力為2.96 MPa，詳細結果列于表7。兩種方法計算所得規律基本一致，僅數值有少量差別。

圖13 地震作用下面板應力(壓為負)Fig.13 Stress along slope direction of concrete face

方法最大順坡向最大壓應力/MPa最小順坡向最大壓應力/MPa最大拉應力/MPaFEM12.25.922.96NSBPEM13.26.912.21

聯合NSBPFE和四分樹網格技術，對面板壩進行跨尺度精細化靜動力數值分析，計算結果符合堆石壩正常變形的基本規律，且FEM計算結果與NSBPFE結果很接近，采用跨尺度方案后，對滿蓄期壩體位移及面板應力變形影響不大，但對地震作用下面板應力影響較為明顯，最大應力差別約1.0 MPa。

總體來看，可以認為NSBPFEM用于跨尺度精細化分析是合理和可行的，可為下一步研究損傷演化和局部破壞提供了技術支撐。

3 結論

采用自主開發的非線性比例邊界多邊形單元(NSBPFE)對典型面板堆石壩系統地進行了靜力、動力計算分析，結果表明：

1)NSBPFE可簡單、便捷地處理三角形、四邊形以及傳統有限元不能直接處理的五邊形、六邊形等多邊形單元，對堆石壩結構進行系統的靜力、動力和永久變形分析，表明該方法結果合理，精度較高，可用于面板壩數值分析。

2)NSBPFE具有更強的靈活性、高效性和通用性優勢，主要體現為該方法可與四分樹網格技術(多節點單元)自動地無縫結合，以進行面板壩高效的跨尺度精細化分析。

3)研究工作為進一步分析面板壩防滲面板結構局部損傷演化和漸進破壞提供了有力的技術支撐，且具有較強的通用性，可以拓展到其他土木和水利工程的安全評價。