基于COMSOL Multiphysics的高面板堆石壩應(yīng)力變形分析

鄧七生

(鄱陽縣水利局 江西上饒 333100)

混凝土面板堆石壩因其對環(huán)境適應(yīng)性強(qiáng)、施工方便、造價(jià)相對較低、抗震性能好等優(yōu)點(diǎn)，在世界范圍內(nèi)得到廣泛應(yīng)用。大壩變形是面板堆石壩的一個(gè)嚴(yán)重問題，因此對混凝土面板堆石壩的應(yīng)力變形進(jìn)行分析具有重要意義[1]。目前，對混凝土面板堆石壩的應(yīng)力和變形分析進(jìn)行了一些研究。黃向志采用三維非線性靜力有限元法分析壩體和防滲體的應(yīng)力和變形。李福忠等采用鄧肯E-B非線性本構(gòu)模型計(jì)算和分析了面板堆石壩的有限元應(yīng)力和變形。張雙梅等人使用ANSYS軟件和有限元法對大壩進(jìn)行了模擬和分析[2]。劉紅等人使用ABAQUS 建立了有限元模型，以分析和研究大壩位移、大壩應(yīng)力和板應(yīng)力。李斌通過數(shù)值模擬分析了復(fù)雜覆蓋層上面板堆石壩的應(yīng)力變形。部分研究對象仍然是100 m 以下的面板堆石壩，隨著施工技術(shù)的不斷進(jìn)步和完善，許多面板堆石大壩的高度已經(jīng)超過100 m。因此，分析高面板堆石壩的應(yīng)力和變形具有重要意義。該文采用Duncan-Zhang E-B 非線性彈性雙曲線模型，利用有限元軟件COMSOL Multiphysics 建立模型，利用MATLAB函數(shù)進(jìn)行二次開發(fā)，對某水庫混凝土面板堆石壩的應(yīng)力和變形進(jìn)行了分析[3]。該水庫流域面積1 915.00 km2，總庫容6.08 億m3，水庫死水位208.00 m，正常水位254.00 m。大壩為混凝土面板堆石壩，壩頂高程275.70 m，最大壩高102.20 m。

1 數(shù)值模擬模型

1.1 堆石體本構(gòu)模型

在有限元數(shù)值模擬中，混凝土面板堆石壩壩體采用鄧肯-張E-B 非線性彈性雙曲線模型。非線性彈性雙曲線模型的計(jì)算公式如下。

式（1）和（2）中，Pa為實(shí)際標(biāo)準(zhǔn)大氣壓力；K為彈性模量；n為彈性模量；σ1為第一主應(yīng)力；σ3為第三主應(yīng)力（圍壓）。

根據(jù)莫爾-庫侖強(qiáng)度準(zhǔn)則：

式（3）中，C和φ均為抗剪強(qiáng)度的指標(biāo)。體積變形模量的計(jì)算公式為

式（4）中，Kb為體積模量系數(shù)；m為體積模量指數(shù)。摩擦φ隨圍壓σ3變化的公式如下：

式（5）中，φ0為σ3=Pa條件下的內(nèi)摩擦角。

1.2 接觸模型

采用無厚度的Goodman單元模擬面板與襯墊之間的接觸狀態(tài)。Goodman單元以兩側(cè)對應(yīng)節(jié)點(diǎn)的相對位移為變量，可以很好地模擬接觸面的位錯(cuò)、滑移或張開，并可以考慮接觸變形的非線性特征[4]。

2 初始應(yīng)力場計(jì)算

根據(jù)以下方法確定新填充層的初始應(yīng)力狀態(tài)：

式（6）和（7）中，σ1和σ3分別為新填充層接合處的主應(yīng)力和次主應(yīng)力；γ為新填土層的重力；h為元素質(zhì)心土壤表面以下的深度。式（8）中，K0為土壤的靜態(tài)側(cè)壓力系數(shù)；φ為材料的內(nèi)摩擦角。

3 基本方程的解

采用分步加載法模擬施工過程中的實(shí)際應(yīng)力應(yīng)變狀態(tài)?；痉匠滩捎煤奢d分級中點(diǎn)增量法求解。將載荷分為若干載荷增量，并對每個(gè)增量進(jìn)行兩次有限元計(jì)算。首先取荷載增量的一半，以求平均應(yīng)力和相應(yīng)的切線模量。然后，將其作為該類荷載增量的平均值，并求解在全荷載增量作用下的位移、應(yīng)力和應(yīng)變增量，該增量疊加在總位移、應(yīng)力和應(yīng)變上[5]。

4 計(jì)算軟件及其二次開發(fā)

為了有效模擬和評估混凝土面板堆石壩施工和運(yùn)行期間的結(jié)構(gòu)安全狀態(tài)，考慮到前處理和后處理的需要以及堆石壩非線性本構(gòu)模型的二次開發(fā)，選擇COMSOL Multiphysics 有限元計(jì)算軟件進(jìn)行研究。在鄧肯拉伸E-B 模型中，堆石材料的切向彈性模量和體積變形模量通過加載和卸載問題計(jì)算，并且應(yīng)力狀態(tài)隨每次加載而變化?；贑OMSOL和MATLAB聯(lián)合開發(fā)函數(shù)，使用MATLAB 函數(shù)計(jì)算每個(gè)加載步驟的切線彈性模量和體積變形模量。通過在COMSOL計(jì)算中自動(dòng)調(diào)用MATLAB 函數(shù)，可以獲得相應(yīng)的切線彈性模量和體積變形模量[6]。對于初始應(yīng)力場的計(jì)算，也基于COMSOL 與MATLAB 聯(lián)合開發(fā)功能，并用MATLAB 進(jìn)行相應(yīng)的二次開發(fā)。

5 有限元模型

考慮盤石頭水庫混凝土面板堆石壩的實(shí)際斷面結(jié)構(gòu)和材料分區(qū)，建立了壩體的二維有限元模型。首先在CAD中建立超元模型（如圖1所示），然后將CAD模型導(dǎo)入COMSOL Multiphysics 有限元軟件，并采用極精細(xì)網(wǎng)格細(xì)分模式自動(dòng)細(xì)分有限元網(wǎng)格（如圖2 所示）。大壩底部固定。細(xì)分單元總數(shù)為20 792個(gè)。大壩填筑分為11個(gè)階段。水負(fù)荷分為5個(gè)階段。

圖1 CAD模型圖

圖2 有限元網(wǎng)格

6 計(jì)算工況和參數(shù)

不同材料的混凝土和堆石體計(jì)算參數(shù)如表1所示。

表1 計(jì)算參數(shù)

7 結(jié)果分析

堆石體的水平和垂直位移云圖以及主要和次要主應(yīng)力云圖如圖3～圖6所示。

圖3 正常蓄水條件下的垂直位移和水平位移云圖

圖4 正常蓄水條件下的主要主應(yīng)力和次要主應(yīng)力云圖

圖5 校核洪水位條件下的垂直位移和水平位移云圖

圖6 校核洪水位條件下的主要主應(yīng)力和次要主應(yīng)力云圖

大壩各種條件下堆石的最大位移和最大應(yīng)力如表2所示。

表2 最大位移和最大應(yīng)力

完工期的板撓度較小，為10.3 cm，出現(xiàn)在覆蓋層頂部的高程處。在正常蓄水條件下，板的最大撓度出現(xiàn)在板的中間，值為34.2 cm，周邊接縫的沉降為1.04 cm，周邊縫的開口位移為1.18 cm。在設(shè)計(jì)洪水位條件下，板的最大撓度出現(xiàn)在板的中間，值為31.4 cm，周邊接縫的沉降為1.22 cm，周邊縫的開口位移為1.30 cm。在校核洪水位條件下，板的最大撓度出現(xiàn)在板的中間，值為53.9 cm，周邊接縫的沉降為1.28 cm，周邊縫的開口位移為1.36 cm。當(dāng)水壓作用在面板上時(shí)，面板法線方向的應(yīng)力基本一致，約等于面板前水壓的大小，其值小于混凝土的抗壓強(qiáng)度。面板主要承受沿壩坡的拉應(yīng)力。在正常蓄水條件下，面板上的最大拉應(yīng)力為2.58 MPa，但周邊接縫附近的位置除外，該位置出現(xiàn)在194.5 m高程的面板表面。從258.37 m 標(biāo)高開始，標(biāo)高以下的拉應(yīng)力超過C25混凝土抗拉強(qiáng)度的設(shè)計(jì)值。在設(shè)計(jì)洪水位條件下，面板上的最大拉應(yīng)力為2.86 MPa，但周邊接縫附近的位置除外，該位置出現(xiàn)在194.5 m高程的面板表面。從261.92 m 標(biāo)高開始，標(biāo)高以下的拉應(yīng)力超過C25 混凝土抗拉強(qiáng)度的設(shè)計(jì)值。在校核洪水位條件下，面板上的最大拉應(yīng)力為3.14 MPa，出現(xiàn)在194.5 m高程的面板表面。從標(biāo)高262.27 m 開始，標(biāo)高以下的拉應(yīng)力超過C25混凝土抗拉強(qiáng)度的設(shè)計(jì)值。

8 結(jié)語

該文基于COMSOL 與MATLAB 的聯(lián)合開發(fā)函數(shù)，利用MATLAB 函數(shù)進(jìn)行二次開發(fā)，模擬分析了不同工況下堆石壩的應(yīng)力和變形特性。計(jì)算結(jié)果與實(shí)際相符，表明模型和計(jì)算方法科學(xué)合理，為相關(guān)工程提供了一種新的計(jì)算方法。有限元計(jì)算和分析表明，由于下游堆石區(qū)的軟巖材料，大壩的沉降明顯從大壩位移分布向下游傾斜。大壩沉降不大，但上下游大壩沉降差較大。周邊接縫在正常范圍內(nèi)變形。大壩中沒有明顯的剪切破壞區(qū)。根據(jù)有限元計(jì)算結(jié)果，在水壓作用下，面板的小主應(yīng)力會(huì)出現(xiàn)超過混凝土抗拉強(qiáng)度的拉應(yīng)力。因此，應(yīng)加強(qiáng)該小組的日常維護(hù)和監(jiān)測。