基于邊界有限元計算方法的多相場重力式大壩動力特性分析研究

張宏美

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

1 概 述

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

2 二維邊界有限元理論

2.1 基本方程

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

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

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

(1)

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

(2)

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

(3)

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

(4)

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

(5)

(6)

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

(7)

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

(8)

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

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

(9)

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

(10)

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

(11)

物理運動方程可統(tǒng)一表述為：

(12)

2.2 動力控制方程

引入能量做功理論來表述動力運動特征，能量參數(shù)的泛函式如下：

(13)

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

(14)

各分項做功分別為：

(15)

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

(16)

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

(17)

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

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

(18)

式中有：

(19)

動力控制方程的求解需要引入對偶變量q：

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

(20)

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

(21)

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

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

(22)

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

(23)

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

(24)

3 工程概況

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

根據(jù)地址鉆孔資料得知，地區(qū)內(nèi)主要為河漫灘平原，地勢較為平坦，最高點為兩側(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)，析水性較強，即滲透系數(shù)較低，分布并不均勻，密實性相差較大，部分地段密實度較高；砂土與卵石相交織，其中取樣表明卵石含量約占該深度取樣總重的51.0%～88.4%，砂土屬全新統(tǒng)，分析表明卵石可能是由于水流搬運作用，從上游補充水源河流沖擊搬運至該處，并持續(xù)沉積。根據(jù)進一步調(diào)查發(fā)現(xiàn)，由于區(qū)域內(nèi)所處板塊間，因而考慮地震動響應(yīng)特征很有必要。

4 多相作用下動力分析

4.1 模型建立

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

4.2 計算結(jié)果分析

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

圖2 數(shù)值模型圖

4.2.1 壩體-流水兩相場

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

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

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

圖4為一組激勵頻率下壩體動水壓力分布曲線。

圖4 一組激勵頻率下壩體動水壓力分布曲線

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

圖5 不同激勵頻率下壩踵處動水壓力變化曲線

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

4.2.2 壩體-水壓力-壩基三相場

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

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

圖6 壩頂位移時程曲線

圖7 下游變坡點位移時程曲線

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

5 結(jié) 論

針對重力式大壩動力響應(yīng)特性問題，引入邊界有限元動力分析方法，計算獲得兩相場與三相場條件下大壩動力特征，結(jié)論如下：

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

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

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

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