趙凱麗，邱流潮，李敬軍

（中國農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院，北京100083）

0 引 言

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

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

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

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

1 鍵基近場動力學(xué)方法

1.1 運動方程

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

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

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

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

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

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

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

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

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

其中伸長率計算可以看做經(jīng)典連續(xù)介質(zhì)力學(xué)中的應(yīng)變值計算如下［29］：

常數(shù)c1的計算公式如下［29］：

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

1.2 數(shù)值離散與時間積分

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

式中：Vj表示物質(zhì)點xj代表的體積。

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

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

2 數(shù)值算例與結(jié)果

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

2.1 集中力作用下的懸臂梁動力響應(yīng)

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

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

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

數(shù)值計算中，懸臂梁一共離散為337 50 個物質(zhì)點，如圖4所示。時間積分步長大小設(shè)置為Δt=1×10-8s，仿真總時間為0.02 s。本算例中使用的材料參數(shù)如下：密度為ρ=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 之后集中力大小不變時，自由端部數(shù)值模擬結(jié)果在1.5 mm 與2.5 mm 之間波動，即在解析解2 mm 上下波動，與實際相符，并將本文鍵基近場動力學(xué)模擬結(jié)果與有限元軟件ABAQUS 的計算結(jié)果進行了比較，圖6表明二者結(jié)果非常吻合，驗證了本文鍵基近場動力學(xué)模型和程序的正確性。

圖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 正弦加速度激勵下垂直方柱的動力響應(yīng)

在結(jié)構(gòu)的地震動力響應(yīng)分析中，需要輸入地震波，本文計算中假設(shè)地基是剛性的，因而在結(jié)構(gòu)基底輸入地震加速度。本文近場動力學(xué)模擬中將已知加速度施加在三層虛擬邊界層的物質(zhì)點上實現(xiàn)地震波輸入，虛擬邊界層的定義如圖7所示。

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

為了驗證虛擬邊界層物質(zhì)點施加加速度模擬地震波的適用性和正確性，對垂直方柱在其底部正弦波激勵作用下的動力響應(yīng)進行數(shù)值模擬計算。方柱高度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

數(shù)值計算中，方柱一共離散為10 800 個物質(zhì)點，如圖10 所示。時間積分步長大小設(shè)置為Δt= 1× 10-5s，仿真總時間為4 s。本算例中使用的材料參數(shù)為：ρ=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 同一時刻相對應(yīng)。圖12 給出了方柱頂部和底部y向的相對位移時程圖，更直觀的顯示出了柱體的動力響應(yīng)，并將本文鍵基近場動力學(xué)模擬結(jié)果與有限元軟件ABAQUS 的計算結(jié)果進行了比較，二者結(jié)果非常吻合，驗證了本文基于鍵基近場動力學(xué)模型的虛擬邊界層物質(zhì)點施加加速度的方法實現(xiàn)地震波輸入是可行的。

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

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

2.3 混凝土重力壩地震動力響應(yīng)

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

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

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

圖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 給出了混凝土重力壩頂水平向和豎向相對位移時程圖的鍵基近場動力學(xué)模擬結(jié)果，數(shù)值計算中采用鍵基近場動力學(xué)模型在虛擬邊界層物質(zhì)點施加加速度的方法實現(xiàn)地震波輸入，從圖像可以看出，壩體的各個部位在不同方向上的動力響應(yīng)情況隨時間是變化的，由于本模擬厚度方向上取一個單元厚度，故暫不考慮z方向上的動力響應(yīng)情況。從圖18可以看出大約在3 s左右壩體橫向和縱向位移開始有明顯波動，與施加的地震波波動大小相對應(yīng)，很好說明了本文鍵基近場動力學(xué)模型模擬混凝土壩地震動力響應(yīng)是可行的。

圖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 結(jié) 論

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

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