張宏美

(浙江湖州宏盈建設(shè)工程有限公司，浙江 湖州 313100 )

1 概 述

水資源是人類生存發(fā)展的基礎(chǔ)，作為水資源的重要調(diào)控設(shè)施，水庫(kù)的穩(wěn)定性是保障水資源高效利用的重要基礎(chǔ)[1-3]。水庫(kù)的重要組成部分即為大壩，水庫(kù)的安全運(yùn)營(yíng)與大壩安全穩(wěn)定性息息相關(guān)。已有很多學(xué)者針對(duì)水庫(kù)大壩的靜力穩(wěn)定性分析作過(guò)相當(dāng)多的研究，極大提高了水庫(kù)大壩的建設(shè)與設(shè)計(jì)水平[4-6]。大壩作為重要水利樞紐工程，抗震能力是不可避免需要考慮的重要方面。不同于靜力分析，大壩動(dòng)力響應(yīng)是在多相條件下，包括壩基與流水狀態(tài)，既有固體場(chǎng)亦有流體場(chǎng)。國(guó)內(nèi)外已有一些學(xué)者通過(guò)振動(dòng)臺(tái)模型試驗(yàn)等分析大壩動(dòng)力反應(yīng)特征[7-10]；當(dāng)然在試驗(yàn)條件不足情況下，上世紀(jì)已有一些學(xué)者通過(guò)解析方法，求解壩體地震反應(yīng)特征[11-12]。上述兩種研究手段，一定程度上均需要消耗大量的人力物力成本，而數(shù)值分析作為高效求解手段，針對(duì)于動(dòng)力反應(yīng)分析具有較高的匹配性。因而，基于多相場(chǎng)的邊界有限元分析方法[13-15]，求解大壩動(dòng)力反應(yīng)特征，為準(zhǔn)確評(píng)估大壩動(dòng)力響應(yīng)穩(wěn)定性提供重要的理論參考。

2 二維邊界有限元理論

2.1 基本方程

在動(dòng)力分析理論中，常常會(huì)引入邊界元分析理論，并結(jié)合有限元離散計(jì)算方法，提升求解速率與精確度，為解決工程解析問(wèn)題提供重要手段。在二維邊界有限元理論中，重定義邊界比例，中心點(diǎn)已不是傳統(tǒng)的零點(diǎn)，其幾何示意圖見圖1。

圖1 比例邊界坐標(biāo)體系

在比例邊界坐標(biāo)體系中，傳統(tǒng)直角坐標(biāo)體系中坐標(biāo)可表述為[16]：

(1)

根據(jù)斜率求導(dǎo)原則，可獲得偏微分方程為：

(2)

而邊界有限元坐標(biāo)體系中，邊界點(diǎn)的坐標(biāo)可表示為：

(3)

引入雅可比矩陣J，將上式以矩陣形式表述為：

(4)

上式即為邊界有限元中運(yùn)動(dòng)矩陣方程，根據(jù)上式可演變出靜力問(wèn)題下的應(yīng)力與變形解，彈性變形下的應(yīng)力與變形方程如下：

(5)

(6)

而變形區(qū)域內(nèi)幾何方程為：

(7)

式中：[L]為微分算子式，其代表式如下：

(8)

根據(jù)有限元插值理論可知，比例邊界坐標(biāo)系中邊界與中心點(diǎn)距離可用徑向插值位移表述，其插值服從下式：

u(ξ,η)=u(ξ)N(η)

(9)

聯(lián)立上述應(yīng)力應(yīng)變以及幾何方程，可得坐標(biāo)系中幾何統(tǒng)一方程為：

(10)

式中:B1、B2代表系數(shù)矩陣，其具體代表式如下：

(11)

物理運(yùn)動(dòng)方程可統(tǒng)一表述為：

(12)

2.2 動(dòng)力控制方程

引入能量做功理論來(lái)表述動(dòng)力運(yùn)動(dòng)特征，能量參數(shù)的泛函式如下：

(13)

根據(jù)材料虛功原理可知，能量泛函式左邊為0，即可得到：

(14)

各分項(xiàng)做功分別為：

(15)

經(jīng)簡(jiǎn)化后可得到：

(16)

式中：系數(shù)矩陣均為

(17)

針對(duì)二維狀態(tài)下的動(dòng)力平衡方程中，引入慣性力來(lái)表述頻域，其控制方程為：

E0ξ2u(ξ),ξξ+(E0+(E1)T-E1)ξu(ξ),ξ-E2u(ξ)+ξ2ω2M0u(ξ)+F(ξ)=0

(18)

式中有：

(19)

動(dòng)力控制方程的求解需要引入對(duì)偶變量q：

q(ξ)=E0ξu(ξ),ξ+(E1)Tu(ξ)

(20)

通過(guò)矩陣轉(zhuǎn)換以及特征值換算，獲得控制方程解為：

(21)

式中：c為常數(shù)矩陣。

由于本文分析的是多相場(chǎng)作用下動(dòng)力反應(yīng)特征，因而劃分單元網(wǎng)格時(shí)考慮多場(chǎng)單元體模型，其形函數(shù)可表述為：

(22)

根據(jù)形函數(shù)表述形式，獲得動(dòng)力特征中應(yīng)力應(yīng)變轉(zhuǎn)換矩陣表達(dá)式為：

(23)

因而，可在上述基礎(chǔ)上，獲得比例邊界坐標(biāo)體系中剛度矩陣K與質(zhì)量矩陣M分別為：

(24)

3 工程概況

某水庫(kù)建設(shè)有混凝土重力壩，壩頂最高處達(dá)80 m，壩頂寬度為3.5 m，已修建有通行公路，上下游坡度均為0.75，水庫(kù)最大庫(kù)容可達(dá)1 200×104m3，承擔(dān)著需水灌溉及工農(nóng)業(yè)生活用水，在枯水期作為空間上水資源調(diào)配的重要利器。大壩包括有液壓控制閘門，啟閉度以及高度均由電腦系統(tǒng)精確控制，輸送水資源量可精確到0.01 m3/s，閘門上下游分別設(shè)置有陡坎與消力池，減弱大水流對(duì)壩體沖蝕作用。壩體整體防滲性較高，通過(guò)在上游壩基覆蓋層以上部分設(shè)置有厚60 cm的防滲墻，壩身上安置有止水面板以及土工格柵，降低水流對(duì)壩體滲漏破壞性。水庫(kù)正常蓄水位設(shè)計(jì)為72.8 m，所有水位設(shè)計(jì)均按照大壩允許承受荷載計(jì)算。由于混凝土施工關(guān)乎著重力式大壩安全可靠性，因而采用BIM技術(shù)對(duì)大壩混凝土施工均進(jìn)行過(guò)多次模擬，獲得混凝土最佳振搗次數(shù)為2次，減少由于混凝土裂縫帶來(lái)的壩體失穩(wěn)與滲漏威脅。

根據(jù)地址鉆孔資料得知，地區(qū)內(nèi)主要為河漫灘平原，地勢(shì)較為平坦，最高點(diǎn)為兩側(cè)的丘陵山坡，壩體兩翼巖體處于較穩(wěn)定狀態(tài)，滑坡及坍塌等地質(zhì)災(zāi)害并不顯著，基巖層為白云質(zhì)灰?guī)r，在局部壩體區(qū)段處還含有夾礫巖，粒徑為15～60 mm，完整性較好，地表出露主要在右側(cè)坡體上。表面覆蓋土層包括人工填土、粉質(zhì)壤土與砂土，其中粉質(zhì)壤土夾有粉、黏粒，含水量達(dá)37%，部分鉆孔取樣顯示有結(jié)核，具有韌性表現(xiàn)，分布范圍基本涵蓋壩體所有區(qū)段，厚度約為2.2～6.9 m；人工填土主要為種植土，較薄，最厚處僅為2.5 m，夾有細(xì)砂結(jié)構(gòu)，析水性較強(qiáng)，即滲透系數(shù)較低，分布并不均勻，密實(shí)性相差較大，部分地段密實(shí)度較高；砂土與卵石相交織，其中取樣表明卵石含量約占該深度取樣總重的51.0%～88.4%，砂土屬全新統(tǒng)，分析表明卵石可能是由于水流搬運(yùn)作用，從上游補(bǔ)充水源河流沖擊搬運(yùn)至該處，并持續(xù)沉積。根據(jù)進(jìn)一步調(diào)查發(fā)現(xiàn)，由于區(qū)域內(nèi)所處板塊間，因而考慮地震動(dòng)響應(yīng)特征很有必要。

4 多相作用下動(dòng)力分析

4.1 模型建立

以壩體平面狀態(tài)作為建?；A(chǔ)，壩基作為剛性變形體，以SOLID65作為其基本單元體模型，壩體采用PLANE42模型，利用ANSYS建立數(shù)值模型。圖2為大壩樁號(hào)1+150～1+160區(qū)段的數(shù)值模型圖，該區(qū)段劃分單元體網(wǎng)格，共獲得18 368個(gè)，節(jié)點(diǎn)數(shù)16 538個(gè)，動(dòng)荷載通過(guò)引入地震荷載開展分析。

4.2 計(jì)算結(jié)果分析

為分析方便，本文給出地震動(dòng)響應(yīng)特征解析解與邊界有限元計(jì)算結(jié)果對(duì)比，為評(píng)價(jià)邊界有限元方法的科學(xué)合理性提供重要依據(jù)。以壩體-水壓力、壩體-壩基的兩相場(chǎng)以及壩體-水壓力-壩基三相場(chǎng)分別進(jìn)行計(jì)算分析，根據(jù)相場(chǎng)特點(diǎn)，與圖2建模類似，分別取其中壩體段開展計(jì)算分析。

圖2 數(shù)值模型圖

4.2.1 壩體-流水兩相場(chǎng)

通過(guò)施加地震動(dòng)荷載后，獲得各階荷載下自振頻率。本文以一階自振頻率下不同反應(yīng)系數(shù)下動(dòng)水壓力的變化曲線作為對(duì)比分析參數(shù)，其計(jì)算結(jié)果見圖3，其中一階自振頻率為32.31 rad/s。

圖3 不同反應(yīng)系數(shù)下動(dòng)水壓力的變化曲線

從圖3可看出，反應(yīng)系數(shù)為0.25時(shí)，曲線相關(guān)系數(shù)達(dá)到0.99，即邊界有限元計(jì)算結(jié)果與解析解基本一致，動(dòng)水壓力整體上呈從壩踵至壩頂，動(dòng)水壓力逐漸降低，在壩踵處，基于邊界有限元方法計(jì)算出的動(dòng)水壓力最大可達(dá)到0.026 5 MPa，與高程0.8H處時(shí)相比來(lái)說(shuō)下降52.8%；當(dāng)反應(yīng)系數(shù)為0.5時(shí)，壩踵處動(dòng)水壓力為0.037 MPa，相比反應(yīng)系數(shù)0.25時(shí)增大39.6%，反應(yīng)系數(shù)0.5時(shí)壩體高程上動(dòng)水壓力分布基本與小反應(yīng)系數(shù)時(shí)基本一致；當(dāng)反應(yīng)系數(shù)增大至0.75，甚至更大的0.925時(shí)，動(dòng)水壓力在高程上分布趨勢(shì)并未發(fā)生較大改變，且壩踵處動(dòng)水壓力分布呈現(xiàn)隨反應(yīng)系數(shù)增大逐漸增大，其中反應(yīng)系數(shù)0.925時(shí)壩踵動(dòng)水壓力相比反應(yīng)系數(shù)0.25時(shí)增大326.4%，達(dá)到0.113 MPa。

圖4為一組激勵(lì)頻率下壩體動(dòng)水壓力分布曲線。

圖4 一組激勵(lì)頻率下壩體動(dòng)水壓力分布曲線

從圖4中可看出，最大動(dòng)水壓力曲線出現(xiàn)在激勵(lì)頻率與自振頻率相等時(shí)，即壩體出現(xiàn)自振現(xiàn)象，動(dòng)水壓力在壩體上分布處于最高，且壩踵處為動(dòng)水壓力中分布最大值。由于壩踵處動(dòng)水壓力分布最大，給出不同激勵(lì)頻率下壩踵處動(dòng)水壓力變化曲線，見圖5。

圖5 不同激勵(lì)頻率下壩踵處動(dòng)水壓力變化曲線

從圖5中可看出，隨著激勵(lì)頻率變化，動(dòng)水壓力整體變化呈先增大后逐漸降低，其中最大動(dòng)水壓力出現(xiàn)在一階自振頻率時(shí)，即頻率5 Hz下動(dòng)水壓力最大，達(dá)0.058 MPa。在該頻率后，壩踵動(dòng)水壓力逐漸降低，但當(dāng)降低至二階自振頻率時(shí)，會(huì)出現(xiàn)陡降現(xiàn)象，即頻率15 Hz時(shí)，下降幅度超過(guò)70%，且在二階自振頻率后動(dòng)水壓力持續(xù)處于較低水平。綜上分析表明，兩相場(chǎng)條件下，壩體動(dòng)水壓力分布受地震動(dòng)荷載頻率與動(dòng)力反應(yīng)系數(shù)影響，壩踵處動(dòng)水壓力乃壩體上最大值。

4.2.2 壩體-水壓力-壩基三相場(chǎng)

前文計(jì)算結(jié)果為兩相場(chǎng)條件下壩體動(dòng)力反應(yīng)特征，但壩體在一定程度上能否抗震，很大程度上取決于壩基的穩(wěn)定性，故而將壩基納入分析動(dòng)力反應(yīng)特性中很有必要。即構(gòu)建起三相場(chǎng)條件下動(dòng)力反應(yīng)體系，壩體區(qū)段取圖2所示的數(shù)值模型，相關(guān)參數(shù)按照土工試驗(yàn)結(jié)果給定。

為分析方便，本文給出大壩兩個(gè)特征點(diǎn)處X、Y相位移時(shí)程曲線，兩個(gè)點(diǎn)分別是壩頂與下游坡度變化點(diǎn)，見圖6和圖7。

圖6 壩頂位移時(shí)程曲線

圖7 下游變坡點(diǎn)位移時(shí)程曲線

從圖6和圖7中可看出，壩頂處X向位移值整體高于Y向位移，在時(shí)間為第二個(gè)峰值點(diǎn)時(shí)，即時(shí)間1.2 s,X向位移達(dá)到0.066 m，Y向位移僅為前者的30.3%；同理，在谷值點(diǎn)時(shí)亦是如此，在第二個(gè)谷值點(diǎn)處，Y向位移達(dá)到負(fù)方向的0.017 m，而X方向相比增長(zhǎng)2.2倍，達(dá)到0.054 m。從X、Y向位移整體來(lái)看，各振蕩次數(shù)下的峰值點(diǎn)位移、谷值點(diǎn)位移之間的周期基本相近，且振蕩表現(xiàn)逐漸趨于穩(wěn)定，即邊界有限元所計(jì)算出來(lái)的位移時(shí)程結(jié)果整體是穩(wěn)定的。相比壩頂處位移結(jié)果，下游變坡點(diǎn)位移整體表現(xiàn)與之一致，不論是X向亦或是Y向，均處于較為穩(wěn)定狀態(tài)；從位移量值來(lái)看，下游變坡點(diǎn)位移值亦是X向位移高于Y向位移，在第二個(gè)峰值位移點(diǎn)時(shí)，X向位移高出Y向位移4.5倍，且谷值點(diǎn)處亦是相差近5倍，即變坡點(diǎn)處X向與Y向位移間差距相比壩頂處更大。

5 結(jié) 論

針對(duì)重力式大壩動(dòng)力響應(yīng)特性問(wèn)題，引入邊界有限元?jiǎng)恿Ψ治龇椒?，?jì)算獲得兩相場(chǎng)與三相場(chǎng)條件下大壩動(dòng)力特征，結(jié)論如下：

1) 基于邊界有限元?jiǎng)恿碚撚?jì)算的結(jié)果與解析解基本一致，不論動(dòng)力反應(yīng)系數(shù)如何變化，兩者結(jié)果均一致，相關(guān)系數(shù)均超過(guò)0.99，邊界有限元理論對(duì)多相場(chǎng)重力壩求解具有較高的自適應(yīng)性。

2) 一階自振頻率下動(dòng)水壓力從壩踵至壩頂逐漸降低，反應(yīng)系數(shù)為0.25時(shí)壩踵動(dòng)水壓力最大，達(dá)0.026 5 MPa，而高程0.8 H時(shí)降低52.8%；當(dāng)動(dòng)力反應(yīng)系數(shù)增大時(shí)，壩踵處最大動(dòng)水壓力逐漸增大，反應(yīng)系數(shù)0.925時(shí)壩踵動(dòng)水壓力相比反應(yīng)系數(shù)0.25時(shí)增大326.4%。

3) 當(dāng)激勵(lì)頻率與自振頻率一致時(shí)，壩體動(dòng)水壓力最高；隨激勵(lì)頻率增大，壩踵處動(dòng)水壓力呈先增大后逐漸降低，其中最高點(diǎn)位于一階自振頻率時(shí)，達(dá)0.058 MPa，二階自振頻率時(shí)動(dòng)水壓力會(huì)出現(xiàn)陡降，降低幅度超過(guò)70%。

4) 三相場(chǎng)條件下X、Y向位移周期與相位均穩(wěn)定，且壩頂與下游變坡點(diǎn)均呈X向位移高于Y向位移，第二個(gè)峰值點(diǎn)時(shí)，壩頂處X向位移達(dá)到0.066 m，Y向位移僅為前者的30.3%，但下游變坡點(diǎn)兩者差距4.5倍，X、Y向位移間差距比壩頂處更大。