基于近場動力學理論的混凝土重力壩地震響應數值模擬

趙凱麗，邱流潮，李敬軍

（中國農業大學水利與土木工程學院，北京100083）

0 引 言

隨著國際能源發展戰略方向的轉移，大力發展水電等可再生能源成為了許多國家戰略發展的新方向。我國如三峽水電站、糯扎渡水電站等大、中型水電站多數位于西部或西南部高烈度地震區［1］，地震作用下的安全問題是影響大壩可持續發展的重要因素，強烈地震作用下水工結構可能會產生結構破壞甚至垮塌的工程事故，從而造成人員傷亡、資金損失等一系列的災難性后果。此外，我國大多數大壩都是在地震現場數據有限的情況下設計的［2］，大壩的除險加固研究也顯得尤為重要。

隨著計算機技術的高速發展，相較物理模型試驗而言，數值模擬方法具有成本低、效率高、簡單易行、可研究材料微觀特性等一系列優點［3］，為大壩的地震動力響應研究提供了一種可行可靠的研究手段。在各種數值計算方法中，有限元方法的理論框架經過幾十年的發展已相對較為成熟，在大壩的地震動力響應研究中被廣泛應用［4-7］。然而，有限元方法是一種基于連續介質力學的方法，在求解非連續問題時仍存在局限性［5］，求解斷裂破壞等非連續問題時會產生奇異性。針對上述不足，相關學者提出了一些改進方法，例如，Zhang 等［8］基于擴展有限元方法研究重力壩在地震作用下的裂紋擴展情況，Fang 等［9］應用擴展有限元方法模擬Koyna 大壩地震裂縫，該方法相對于有限元法有所改進，可以解決弱不連續問題，但是難以模擬多裂紋擴展問題。Das等［10］基于無網格法預測地震作用下大壩潰壩前破壞機制，該方法在處理網格畸變問題時具有很大優勢，克服了對網格的依賴性，然而其計算量較大且在模擬多裂紋擴展時受到限制［11］。Pekau 等［12］基于非連續介質力學的離散元法對大壩進行裂縫分析，成功模擬了開裂問題，但是需要明確指出裂紋萌生點，在復雜的地震荷載作用下難以實現，并且計算量較大，因此需要一種新的方法來克服上述問題。

本文采用近場動力學（peridynamics，PD）方法進行混凝土重力壩地震響應數值模擬。該方法由美國Sandia 實驗室的Silling［13］于2000年提出，黃丹等［14］2010年首次在國內應用。近場動力學法［15，16］是一種聯系經典連續介質力學和分子動力學的數學力學理論，是分子動力學的連續版本，同時兼有分子動力學和無網格法的優點，在解決斷裂等非連續大變形問題方面具有明顯優勢。近場動力學理論包括鍵基（bond-based）近場動力學理論和態基（state-based）近場動力學理論［17］，該理論已被應用到了許多研究領域，包括邊坡穩定性分析［18］、非均質材料的水力壓裂破壞模擬［19］、冰塊撞擊破壞［20］、巖土類材料［21］的大變形研究以及混凝土結構破壞過程模擬［22-24］等。Gu 等［25］運用近場動力學鍵基理論研究了混凝大壩沖擊破壞；張鈺彬等［26］模擬了高水頭作用下混凝土重力壩的水力劈裂過程，為研究大壩破壞機理提供了新方法，具有重要的工程意義。

本文應用鍵基近場動力學理論對混凝土重力壩的地震動力響應進行了數值研究。論文的余下部分安排如下：在第一節中簡要介紹了鍵基近場動力學理論及數值計算方法；在第二節中首先通過懸臂梁自由端施加集中力和底部輸入周期運動的正弦波時的動力響應兩個算例來驗證模型和程序的精度和有效性，隨后基于驗證過的模型對混凝土重力壩的地震動力響應進行數值模擬；在第三節中列出了本文研究得出的一些結論和進一步研究的方向。

1 鍵基近場動力學方法

1.1 運動方程

近場動力學理論不再基于連續性假設以及微分方程求解力學問題［27］，而是采用非局部作用的積分形式的運動方程，避免了非連續奇異性問題的發生。非局部思想反映在物質的相互作用方面，如圖1 所示，在某一時刻，空間域R內物體的任一物質點xi與以其為中心，近場范圍δ（horizon）為半徑的圓形平面或球體范圍內的其他物質點xj之間產生相互作用，即當|xi-xj|≤δ時，兩物質點之間產生相互作用力，而在近場范圍內的物質點xj構成xi的家族（familiy），用H表示。

圖1 物質點xi和xj變形過程中的相互作用Fig.1 Interaction between material points xi and xj during deformation

根據牛頓第二定律可以得到物質點xi的基于鍵基近場動力學理論的運動方程［28］：

式中：ρ表示密度；表示加速度；dV是物質點xj處的體積微元；b表示單位體積物質上施加的外部荷載，即體力密度為力的密度矢量，表示在t時刻物質點之間的相互作用。

對于鍵基近場動力學理論，物質點之間的相互作用力大小相等，方向相反，如圖1所示。力密度矢量形式［29］如下：

式中：C是一個與物質點xi和xj之間相對伸長量以及近場動力學范圍有關的輔助參數和表示物質點xi和xj變形后的位置。則有物質點xi的運動方程為：

式中：f表示xi和xj之間的相互作用，用“鍵”來表示，由其構成的函數稱作本構力函數。

兩物質點之間的相對位置用ξ=xj-xi來表示，相對位移表示為：η=uj-ui，式（4）可進一步寫為：

對于各向同性材料，可以假設力密度矢量f與物質點之間的拉伸成線性相關［29，30］：

其中伸長率計算可以看做經典連續介質力學中的應變值計算如下［29］：

常數c1的計算公式如下［29］：

式中：c為鍵常數微觀模量；體積模量，E是彈性模量，υ是泊松比。

1.2 數值離散與時間積分

近場動力學方法求解積分形式的運動方程，通常將所研究物體離散為一系列物質點，于是得到控制方程（4）在空間域的離散表達式為：

式中：Vj表示物質點xj代表的體積。

本文對（9）式采用顯式的向前差分時間積分求解，即按下列步驟計算：

式中：n表示時間步數；Δt是時間積分步長，tn+1=tn+Δt，對于顯式時間積分法，時間步長的選取需要滿足穩定條件［30］。

2 數值算例與結果

本文通過使用FORTRAN 90 語言編程實現了近場動力學的數值計算。首先通過懸臂梁自由端受集中力響應以及在底部輸入周期運動時動力響應兩個算例驗證來模型和程序的精度和有效性，接著對混凝土重力壩的地震動力響應進行了數值模擬。

2.1 集中力作用下的懸臂梁動力響應

這里首先對懸臂梁自由端在集中力作用下的動力響應進行模擬。懸臂梁的長度L=0.02 m，寬度W=0.002 m 和高度H=0.002 m，如圖2 所示。懸臂梁自由端施加動荷載P隨時間變化的關系如圖3所示。

圖2 懸臂梁尺寸及示意圖（單位：mm）Fig.2 Dimensions and schematic diagram of cantilever beam

圖3 荷載與時間之間的關系圖Fig.3 The relation diagram between load and time

數值計算中，懸臂梁一共離散為337 50 個物質點，如圖4所示。時間積分步長大小設置為Δt=1×10-8s，仿真總時間為0.02 s。本算例中使用的材料參數如下：密度為ρ=1 000 kg/m3，彈性模量E=1 GPa，泊松比為υ= 0.25。本算例將集中力作為體力密度施加在懸臂梁最外層單元上，并且暫時不考慮阻尼影響。

圖4 懸臂梁初始時刻的離散效果圖Fig.4 Discrete rendering of the cantilever at the initial moment

圖5 給出了t=5×10-4s、t=2×10-3s、t=5.5×10-3s 和t=2×10-2s四個時刻懸臂梁的豎向位移云圖。圖6給出了懸臂梁自由端的豎向位移時程圖，由于不考慮阻尼影響，1×10-3s 之后集中力大小不變時，自由端部數值模擬結果在1.5 mm 與2.5 mm 之間波動，即在解析解2 mm 上下波動，與實際相符，并將本文鍵基近場動力學模擬結果與有限元軟件ABAQUS 的計算結果進行了比較，圖6表明二者結果非常吻合，驗證了本文鍵基近場動力學模型和程序的正確性。

圖5 懸臂梁在不同時刻z方向上的位移Fig.5 The displacement of the cantilever beam in the z direction at differenct moments

圖6 自由端位移曲線比較Fig.6 The comparison of free end displacement curves

2.2 正弦加速度激勵下垂直方柱的動力響應

在結構的地震動力響應分析中，需要輸入地震波，本文計算中假設地基是剛性的，因而在結構基底輸入地震加速度。本文近場動力學模擬中將已知加速度施加在三層虛擬邊界層的物質點上實現地震波輸入，虛擬邊界層的定義如圖7所示。

圖7 虛擬邊界層的定義Fig.7 The definition of virtual boundary layer

為了驗證虛擬邊界層物質點施加加速度模擬地震波的適用性和正確性，對垂直方柱在其底部正弦波激勵作用下的動力響應進行數值模擬計算。方柱高度L=1.0 m，寬度W=0.02 m，厚度H=0.02 mm，如圖8 所示。在本算例中，在方柱的底部虛擬邊界三層粒子的y方向上施加總時長為4 s 的正弦加速度如圖9所示。

圖8 垂直懸臂柱體尺寸圖（單位：m）Fig.8 The dimension of vertical cantilever column

圖9 正弦加速度Fig.9 Sinusoidal acceleration

數值計算中，方柱一共離散為10 800 個物質點，如圖10 所示。時間積分步長大小設置為Δt= 1× 10-5s，仿真總時間為4 s。本算例中使用的材料參數為：ρ=7 780 kg/m3，彈性模量E=1 GPa，泊松比為υ= 0.25。計算中不考慮阻尼影響。

圖10 垂直懸臂柱體離散圖Fig.10 The discrete diagram of vertical cantilever column

圖11 給出了t=1 s、t=2 s、t=3 s 和t=4 s 共4 個時刻方柱在y向位移云圖，顯示了柱體在施加周期波后在y方向上的擺動情況，且從圖像可以看出懸臂柱體在以上4 個時刻頂部位移要大于底部位移，與圖12 同一時刻相對應。圖12 給出了方柱頂部和底部y向的相對位移時程圖，更直觀的顯示出了柱體的動力響應，并將本文鍵基近場動力學模擬結果與有限元軟件ABAQUS 的計算結果進行了比較，二者結果非常吻合，驗證了本文基于鍵基近場動力學模型的虛擬邊界層物質點施加加速度的方法實現地震波輸入是可行的。

圖11 懸臂梁在不同時刻y方向上的位移Fig.11 The displacement of the cantilever in the y direction at different moments

圖12 近場動力學模型與ABAQUS軟件計算得到的柱頂部和柱底部相對位移Fig.12 The relative displacement of column top and column bottom was calculated by the peridynamics model and ABAQUS software

2.3 混凝土重力壩地震動力響應

為了驗證本文鍵基近場動力學模型模擬混凝土重力壩地震動力響應的可行性，本小節對Koyna 混凝土重力壩在水平及豎向地震波共同作用下的動力響應進行模擬。研究壩斷面高度103 m，底寬70 m，頂寬14.8 m，原壩型形狀及尺寸如圖13所示。

圖13 大壩基本尺寸示意圖（單位：m）Fig.13 The schematic diagram of basic dam dimensions

數值計算中按照原壩型進行建模，為簡化計算，將厚度方向取0.5 m，為一個單元高度大小，混凝土重力壩一共離散為28 840 個物質點，如圖14 所示。時間積分步長大小設置為Δt=1× 10-5s，仿真總時間為10 s。本算例中使用的材料參數為：ρ=2 643 kg/m3，彈性模量E=31.02 GPa，泊松比為υ= 0.25。在本算例中，在壩的底部虛擬邊界三層粒子的x及y方向上施加圖15所示的地震加速度。計算不考慮庫水及阻尼的影響，同時假設壩基為剛體。

圖14 壩體離散圖Fig.14 The discrete diagram of the dam

圖15 地震加速度Fig.15 seismic acceleration

圖16 和17 分別給出了t=2 s、t=4 s、t=8 s 和t=10 s 共4 個時刻混凝土壩水平向和豎向位移云圖。圖18 給出了混凝土重力壩頂水平向和豎向相對位移時程圖的鍵基近場動力學模擬結果，數值計算中采用鍵基近場動力學模型在虛擬邊界層物質點施加加速度的方法實現地震波輸入，從圖像可以看出，壩體的各個部位在不同方向上的動力響應情況隨時間是變化的，由于本模擬厚度方向上取一個單元厚度，故暫不考慮z方向上的動力響應情況。從圖18可以看出大約在3 s左右壩體橫向和縱向位移開始有明顯波動，與施加的地震波波動大小相對應，很好說明了本文鍵基近場動力學模型模擬混凝土壩地震動力響應是可行的。

圖16 重力壩在不同時刻沿x方向的位移Fig.16 The displacement of the gravity dam in the x direction at different moments

圖17 重力壩在不同時刻沿y方向的位移Fig.17 The displacement of the gravity dam in the y direction at different moments

圖18 壩頂相對位移時程圖Fig.18 The relative displacement time history of dam crest

3 結 論

本文建立了模擬混凝土重力壩地震動力響應的鍵基近場動力學模型，并采用FORTRAN90語言編程實現。其中，地震波的輸入是通過在三層虛擬節點施加地震加速度來實現的。本文通過對比懸臂梁自由端和底部受動力和周期運動正弦波的近場動力學數值模型與ABAQUS 軟件的模擬結果，證明該近場動力學模型適宜模擬動力響應問題。接著利用本文模型對混凝土重力壩的地震動力響應進行數值模擬，結果表明：本文鍵基近場動力學模型模擬混凝土重力壩的地震動力響應是可行的，為混凝土壩地震動力響應數值模擬方法提供了新的思路。

近場動力學需要計算近場范圍內的所有物質點的相互作用，較有限元等傳統的網格方法需要計算的物質點多，因此計算時間較長，今后有待于進一步開展并行算法研究以提高計算效率。本文選用較為簡單、應用較廣的近場動力學鍵基模型，存在泊松比的限制，一些情況下可能產生誤差，后續可以采用態基近場動力學模型進行研究。另外，本文計算中沒有考慮壩基的變形以及輻射阻尼的影響，將在后期研究中進一步考慮。