基于近場(chǎng)動(dòng)力學(xué)理論的混凝土重力壩地震響應(yīng)數(shù)值模擬

趙凱麗，邱流潮，李敬軍

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

0 引 言

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

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

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

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

1 鍵基近場(chǎng)動(dòng)力學(xué)方法

1.1 運(yùn)動(dòng)方程

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖4 懸臂梁初始時(shí)刻的離散效果圖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四個(gè)時(shí)刻懸臂梁的豎向位移云圖。圖6給出了懸臂梁自由端的豎向位移時(shí)程圖，由于不考慮阻尼影響，1×10-3s 之后集中力大小不變時(shí)，自由端部數(shù)值模擬結(jié)果在1.5 mm 與2.5 mm 之間波動(dòng)，即在解析解2 mm 上下波動(dòng)，與實(shí)際相符，并將本文鍵基近場(chǎng)動(dòng)力學(xué)模擬結(jié)果與有限元軟件ABAQUS 的計(jì)算結(jié)果進(jìn)行了比較，圖6表明二者結(jié)果非常吻合，驗(yàn)證了本文鍵基近場(chǎng)動(dòng)力學(xué)模型和程序的正確性。

圖5 懸臂梁在不同時(shí)刻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 正弦加速度激勵(lì)下垂直方柱的動(dòng)力響應(yīng)

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

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

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

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

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

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

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

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

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

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

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

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

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

數(shù)值計(jì)算中按照原壩型進(jìn)行建模，為簡(jiǎn)化計(jì)算，將厚度方向取0.5 m，為一個(gè)單元高度大小，混凝土重力壩一共離散為28 840 個(gè)物質(zhì)點(diǎn)，如圖14 所示。時(shí)間積分步長(zhǎng)大小設(shè)置為Δt=1× 10-5s，仿真總時(shí)間為10 s。本算例中使用的材料參數(shù)為：ρ=2 643 kg/m3，彈性模量E=31.02 GPa，泊松比為υ= 0.25。在本算例中，在壩的底部虛擬邊界三層粒子的x及y方向上施加圖15所示的地震加速度。計(jì)算不考慮庫(kù)水及阻尼的影響，同時(shí)假設(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 個(gè)時(shí)刻混凝土壩水平向和豎向位移云圖。圖18 給出了混凝土重力壩頂水平向和豎向相對(duì)位移時(shí)程圖的鍵基近場(chǎng)動(dòng)力學(xué)模擬結(jié)果，數(shù)值計(jì)算中采用鍵基近場(chǎng)動(dòng)力學(xué)模型在虛擬邊界層物質(zhì)點(diǎn)施加加速度的方法實(shí)現(xiàn)地震波輸入，從圖像可以看出，壩體的各個(gè)部位在不同方向上的動(dòng)力響應(yīng)情況隨時(shí)間是變化的，由于本模擬厚度方向上取一個(gè)單元厚度，故暫不考慮z方向上的動(dòng)力響應(yīng)情況。從圖18可以看出大約在3 s左右壩體橫向和縱向位移開始有明顯波動(dòng)，與施加的地震波波動(dòng)大小相對(duì)應(yīng)，很好說(shuō)明了本文鍵基近場(chǎng)動(dòng)力學(xué)模型模擬混凝土壩地震動(dòng)力響應(yīng)是可行的。

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

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

圖18 壩頂相對(duì)位移時(shí)程圖Fig.18 The relative displacement time history of dam crest

3 結(jié) 論

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

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