基于無網(wǎng)格界面模擬方法的面板壩防滲體跨尺度分析

鄒德高，龔 瑾，孔憲京，劉京茂，屈永倩

（1.大連理工大學(xué) 海岸和近海工程國家重點(diǎn)試驗(yàn)室，遼寧 大連 116024；2.大連理工大學(xué) 水利工程學(xué)院，遼寧 大連 116024）

1 研究背景

面板堆石壩在近些年得到迅猛的發(fā)展[1]，我國也正處在從200 m級到300 m級面板堆石壩的關(guān)鍵研發(fā)階段。隨著工程規(guī)模的逐步增加，給以有限元為主的數(shù)值模擬方法帶來巨大挑戰(zhàn)[2]?；炷撩姘遄鳛榇髩侮P(guān)鍵的防滲結(jié)構(gòu)，其安全性態(tài)至關(guān)重要。因此混凝土面板需采用較密的網(wǎng)格模擬其應(yīng)力分布[3]，然而面板與堆石體尺寸相差懸殊，若在堆石體區(qū)采用一致尺寸的網(wǎng)格，無疑大大的增加了不必要的計(jì)算成本，因此跨尺度建模對面板壩精細(xì)化分析有著重要意義。

在面板壩跨尺度模型中，傳統(tǒng)的點(diǎn)對點(diǎn)界面單元（Goodman單元，薄層單元）[4-5]無法連接不同尺寸的面板網(wǎng)格與墊層區(qū)網(wǎng)格。針對上述問題，學(xué)者們做出相應(yīng)的研究。周墨臻等[6-7]通過接觸力學(xué)算法代替界面單元實(shí)現(xiàn)荷載的在疏密網(wǎng)格間的傳遞，但該方法較難描述界面處在地震荷載下復(fù)雜的變形[8]。鐘紅等[9]采取逐漸過渡的方法實(shí)現(xiàn)跨尺度建模，但是這種方法需要大量的人工干預(yù)且不適合三維模型。Chen 等[10-12]結(jié)合比例邊界有限元和四分樹技術(shù)實(shí)現(xiàn)了在堆石體區(qū)域的疏密網(wǎng)格的自動拋分。Qu等[13-14]開發(fā)的非對稱界面單元實(shí)現(xiàn)了兩側(cè)不同尺寸網(wǎng)格的連接，但該界面單元部分節(jié)點(diǎn)仍需對應(yīng)且不適合局部網(wǎng)格優(yōu)化等問題。

若能實(shí)現(xiàn)界面兩側(cè)節(jié)點(diǎn)自由分布，擺脫點(diǎn)對點(diǎn)限制，則可彌補(bǔ)上述方法的不足，解決面板壩中跨尺度的難題。而與有限元相比，只需要節(jié)點(diǎn)信息的無網(wǎng)格法更能滿足要求。最早出現(xiàn)的無網(wǎng)格法是由Lucy提出的SPH法[15]，應(yīng)用于天體問題。1994年Belystschko等[16]提出了基于背景網(wǎng)格的EFG法，隨后無網(wǎng)格法快速發(fā)展，并在裂紋擴(kuò)展、流固耦合和液化大變形等方面[17-20]得到了廣泛的應(yīng)用。

本文將無網(wǎng)格法[21-22]拓展到土與結(jié)構(gòu)界面的模擬中，以實(shí)現(xiàn)界面兩側(cè)節(jié)點(diǎn)自由分布且無需對應(yīng)，該方法可直接連接尺寸不同的網(wǎng)格，從而簡化了跨尺度建模及局部網(wǎng)格優(yōu)化的過程。同時(shí)基于該方法對一200 m級面板堆石壩進(jìn)行跨尺度分析，并與傳統(tǒng)Goodman單元對比，證明該方法在面板壩靜力分析和動力分析中具備足夠的精度。同時(shí)建立防滲面板與墊層區(qū)跨尺度分析模型，通過對比計(jì)算結(jié)果及計(jì)算時(shí)間，證明無網(wǎng)格界面在保證面板模擬精度的前提下，可大幅度減少堆石體的網(wǎng)格，從而提高計(jì)算效率。該方法很容易擴(kuò)展到三維并為面板壩面板精細(xì)化損傷演化分析提供有效的技術(shù)手段。

2 基于無網(wǎng)格的非點(diǎn)對點(diǎn)界面模擬方法

2.1 理論概要及相關(guān)公式推導(dǎo)基于無網(wǎng)格的非點(diǎn)對點(diǎn)界面，兩側(cè)節(jié)點(diǎn)自由分布且無需對應(yīng)。高斯點(diǎn)的生成和積分過程均在引入的背景網(wǎng)格線上進(jìn)行。圖1、圖2為一個(gè)典型無網(wǎng)格界面及高斯點(diǎn)分布。根據(jù)相對距離，確定高斯點(diǎn)支持域內(nèi)節(jié)點(diǎn)，并通過徑向插值函數(shù)（RPIM）[22]計(jì)算相應(yīng)插值點(diǎn)位移。圖3所示為一高斯點(diǎn)分別在界面兩側(cè)支持域內(nèi)的節(jié)點(diǎn)。插值點(diǎn)在界面上、下兩側(cè)位移可表達(dá)為（這里僅以上側(cè)位移為例）：

圖1 基于無網(wǎng)格的非點(diǎn)對點(diǎn)界面

圖2 背景網(wǎng)格線及高斯點(diǎn)分布

圖3 高斯點(diǎn)支持域內(nèi)節(jié)點(diǎn)

式中：ai、bj為待求常數(shù)；l為各節(jié)點(diǎn)局部坐標(biāo)；rk為各節(jié)點(diǎn)與高斯點(diǎn)距離；Bi（rk）為徑向基函數(shù)（RBF），常見的徑向基函數(shù)有許多種，本文采用MQ基函數(shù)，MQ基函數(shù)中包含c、q兩個(gè)形狀參數(shù)，根據(jù)Liu等[23]的數(shù)值試驗(yàn)，本文中c=3，q=1.03；Pj（l）為在徑向基函數(shù)中增加的線性多項(xiàng)式，增加的線性多項(xiàng)式為[1l]；n為徑向基函數(shù)個(gè)數(shù)，即表示高斯點(diǎn)支持域內(nèi)節(jié)點(diǎn)數(shù)量；m為多項(xiàng)式基函數(shù)個(gè)數(shù)，本文中m=2。

式中：dav為支持域內(nèi)節(jié)點(diǎn)平均間距。

聯(lián)立式（1）和式（2）可得到如下矩陣方程：

上式中各矩陣表達(dá)式為：

將式（4）—式（6）帶入式（3）中，求解可得：

將式（7）和式（8）帶入式（1），則插值點(diǎn)在界面上側(cè)位移可表示為：

式中：φi為插值點(diǎn)在界面上側(cè)支持域內(nèi)各節(jié)點(diǎn)處形函數(shù)；ue為各節(jié)點(diǎn)處位移插值點(diǎn)處相對位移，可表示為：

由此可得到應(yīng)變轉(zhuǎn)化矩陣Bi：

式中：φi、φi'分別代表插值點(diǎn)在界面上、下兩側(cè)支持域內(nèi)各節(jié)點(diǎn)處形函數(shù)；n、n′分別代表插值點(diǎn)在界面兩側(cè)支持域內(nèi)節(jié)點(diǎn)個(gè)數(shù)。

求得應(yīng)變矩陣后，引入界面本構(gòu)[24]矩陣Di，并在背景網(wǎng)格線內(nèi)進(jìn)行高斯積分得到剛度矩陣：

以上公式均在局部坐標(biāo)下求得，通過坐標(biāo)轉(zhuǎn)化矩陣可將上述剛度陣、內(nèi)外力矩陣轉(zhuǎn)化到整體坐標(biāo)系下。最后根據(jù)節(jié)點(diǎn)編號，將上述矩陣組裝整體剛度陣及不平衡力中，由平衡方程求解各節(jié)點(diǎn)處位移：

從上述的推導(dǎo)中可以看出，無網(wǎng)格界面與Goodman 界面單元主要的區(qū)別在于高斯點(diǎn)位移的插值模式分別為徑向插值函數(shù)和線性插值及界面兩側(cè)節(jié)點(diǎn)分布上，兩種方法均可通過界面本構(gòu)矩陣來反映土體與結(jié)構(gòu)之間的切向、法向接觸關(guān)系。

2.2 改進(jìn)RPIM插值函數(shù)徑向插值函數(shù)（RPIM）在求解域內(nèi)部具有很高精度，然而在邊界附近由于一側(cè)支持域內(nèi)節(jié)點(diǎn)減少，RPIM插值精度有所降低。通過引入虛設(shè)節(jié)點(diǎn)，結(jié)合有限元線性插值的優(yōu)勢可改進(jìn)RPIM插值函數(shù)，圖4所示為該方法的主要流程。

（1）首先在求解域內(nèi)邊界增設(shè)虛節(jié)點(diǎn)。

（2）通過RPIM計(jì)算虛節(jié)點(diǎn)處位移，具體過程參考式（1）—式（10），虛節(jié)點(diǎn)處位移可表示為：

（3）判斷插值點(diǎn)位置，如果插值點(diǎn)位于求解域內(nèi)部，依然通過RPIM求該點(diǎn)位移，若插值點(diǎn)位于邊界附近，則通過虛節(jié)點(diǎn)和邊界處節(jié)點(diǎn)線性插值計(jì)算：

將式（15）帶入式（16），最終插值點(diǎn)位移可由各節(jié)點(diǎn)處位移表示：

由改進(jìn)的RPIM插值函數(shù)計(jì)算得到節(jié)點(diǎn)3處的形函數(shù)如圖5所示。修正后的RPIM繼承原插值函數(shù)基本性質(zhì)，即：（1）δ函數(shù)性質(zhì)；（2）單位分解性；（3）再生性；（4）緊致性；（5）連續(xù)性。

通過曲線擬合測試驗(yàn)證改進(jìn)RPIM插值函數(shù)精度。本文選擇函數(shù)y=sinπx，在[-2，2]區(qū)間設(shè)置11個(gè)均勻的計(jì)算點(diǎn)并求其函數(shù)值，通過與理論值對比，圖6繪制了兩種插值函數(shù)誤差曲線。表1總結(jié)了在（-2，-1.9]區(qū)間兩者誤差對比。

改進(jìn)的RPIM 插值函數(shù)結(jié)合RPIM和有限元線性插值的優(yōu)點(diǎn)，在保證原插值函數(shù)優(yōu)勢的前提下，進(jìn)一步改善了邊界附近的插值的精度，并將改進(jìn)的RPIM插值函數(shù)引入非點(diǎn)對點(diǎn)界面中。文中無網(wǎng)格界面插值函數(shù)均采用改進(jìn)的RPIM計(jì)算。

圖4 改進(jìn)的RPIM插值函數(shù)實(shí)現(xiàn)過程

圖5 形函數(shù)對比

圖6 誤差分析對比

表1 兩種插值方法在（-2，-1.9]誤差對比

2.3 數(shù)值程序?qū)崿F(xiàn)本文基于課題組自主開發(fā)的GEODYNA 軟件平臺，采用面向?qū)ο蟪绦蛟O(shè)計(jì)思想及C++語言，將基于無網(wǎng)格的界面模擬方法封裝成一個(gè)超單元并整合到GEODYNA中，實(shí)現(xiàn)與FEM無縫耦合，并可應(yīng)用復(fù)雜彈塑性界面本構(gòu)。非點(diǎn)對點(diǎn)界面程序框架如圖7所示。首先抽象出其與有限元的共同特征（計(jì)算形函數(shù)，求解剛度陣、質(zhì)量陣、阻尼陣等），進(jìn)而設(shè)置相應(yīng)的成員函數(shù)及內(nèi)部繼承關(guān)系，其次根據(jù)非點(diǎn)對點(diǎn)界面自身的性質(zhì)，設(shè)置外部調(diào)用接口及特有的成員函數(shù)。最終將基于無網(wǎng)格的非點(diǎn)對點(diǎn)界面在限元框架下實(shí)現(xiàn)。

3 面板堆石壩中的應(yīng)用

3.1 算例模型某面板堆石壩，壩高200 m，上游坡比1∶1.5，下游坡比1∶1.6，壩頂寬度18 m，面板頂部寬度0.3 m，壩體部分為主堆石區(qū)及墊層區(qū)。為保證面板在靜、動力作用下的模擬精度，面板沿厚度方向分為5層單元，沿壩坡方向長度2 m。面板與墊層界面分別用有厚度Goodman單元及基于無網(wǎng)格的非點(diǎn)對點(diǎn)界面來進(jìn)行模擬，兩種方法均不同于傳統(tǒng)無厚度Goodman單元[4]，法向剛度不再是個(gè)很大的常數(shù)[25]，其法向剛度可與周圍土體剛度大小相兼容，可較好地模擬接觸帶的法向剪縮、剪脹變形（或法向剛度）特性。算例局部網(wǎng)格如圖8所示。

3.2 材料參數(shù)在靜、動力計(jì)算中，堆石料及墊層材料均采用堆石料廣義塑性模型[26-27]，材料參數(shù)分別列于表2、表3，表2、表3中G0、K0、ms、mv為與彈性模量相關(guān)系數(shù)，Mg、Mf、af、ag為與塑性流動方向相關(guān)系數(shù)，其余為與塑性模量相關(guān)系數(shù)。面板與墊層間的界面采用界面廣義塑性模型[28]，參數(shù)見表4，表4中Ds0、Dn0為與彈性模量相關(guān)系數(shù)，km、Mf、α、rd為與塑性流動方向相關(guān)系數(shù)，a、b、c為與顆粒破相關(guān)系數(shù)，Mc、er0、λ為與臨界狀態(tài)相關(guān)系數(shù)，k、H0、fh為與塑性模量相關(guān)系數(shù)。在界面廣義塑性本構(gòu)模型中可直接引入厚度t來考慮接觸帶沿厚度方向的特性。面板、趾板及基巖部分采用線彈性模型，參數(shù)見表5。靜力計(jì)算時(shí)，壩體經(jīng)40個(gè)荷載步完成填筑，并分38個(gè)荷載步蓄水至190 m高程。

圖7 基于無網(wǎng)格的非點(diǎn)對點(diǎn)界面程序框架

圖8 CFRD網(wǎng)格對比

表2 堆石料廣義塑性模型參數(shù)

表3 墊層料廣義塑性模型參數(shù)

表4 界面廣義塑性模型參數(shù)

表5 線彈性材料參數(shù)

3.3 地震動輸入根據(jù)現(xiàn)行的《水工建筑物抗震設(shè)計(jì)規(guī)范》規(guī)定的標(biāo)準(zhǔn)反應(yīng)譜，生成水平向及豎向人工模擬地震波。順河向地震動峰值為0.3g，豎向地震動峰值為0.2g，地震波持時(shí)25 s。地震波加速度時(shí)程如圖9所示，地震采用波動法輸入，計(jì)算步長Δt=0.01s。

圖9 地震動加速度時(shí)程

3.4 計(jì)算結(jié)果分析首先給出靜力計(jì)算結(jié)果，兩種模型在滿蓄期主應(yīng)力及面板順坡向應(yīng)力的對比。圖10為堆石壩滿蓄期壩體最大、最小主應(yīng)力。圖11為滿蓄期面板（面板中間層單元）順坡向應(yīng)力（壓為負(fù)）對比。

無網(wǎng)格界面實(shí)現(xiàn)了兩側(cè)節(jié)點(diǎn)自由分布且無需對應(yīng)。圖10表明，無網(wǎng)格界面與Goodman單元兩種計(jì)算方法相差較小。證明了無網(wǎng)格界面在面板堆石壩靜力分析中具備足夠精度。

同時(shí)將無網(wǎng)格界面應(yīng)用于面板壩動力分析中。圖12為在地震荷載作用下下面板最大、最小順坡向應(yīng)力（壓為負(fù)），圖13為面板頂部加速度時(shí)程曲線。

由圖12、圖13可以看出，兩種算法在地震中面板應(yīng)力相差較小，面板頂部加速度時(shí)程一致，證明了基于無網(wǎng)格的非點(diǎn)對點(diǎn)界面在面板壩動力分析的可行性。

圖10 滿蓄期壩體應(yīng)力（單位：MPa）

圖11 滿蓄期面板順坡向應(yīng)力

圖12 地震作用下面板順坡向應(yīng)力

圖13 面板頂部A點(diǎn)加速度時(shí)程曲線

圖14 不同比例的面板和墊層跨尺度模型

圖15 面板順坡向應(yīng)力

表6 平均相對偏差

表7 計(jì)算效率對比

3.5 高效跨尺度建模及分析面板采用精細(xì)化網(wǎng)格模擬后，由于在堆石體區(qū)采用一致尺寸的網(wǎng)格，總單元數(shù)達(dá)31 768 個(gè)，無疑增加了不必要的計(jì)算成本。非點(diǎn)對點(diǎn)界兩側(cè)節(jié)點(diǎn)自由分布且無需對應(yīng)，因此可連接不同尺寸的面板與墊層區(qū)網(wǎng)格，進(jìn)而實(shí)現(xiàn)面板壩的跨尺度分析。在圖8（c）基礎(chǔ)上保持面板網(wǎng)格不變，將堆石體用相對較疏的網(wǎng)格替換，通過替換不同尺寸的網(wǎng)格，實(shí)現(xiàn)不同比例的跨尺度建模。圖14為通過非點(diǎn)對點(diǎn)界面實(shí)現(xiàn)的面板和墊層順坡向尺度跨越（分別為1∶2和1∶3）。

采用上述的材料及地震輸入，分別對不同比例的跨尺度模型進(jìn)行靜動力分析。圖15為各模型面板順坡向應(yīng)力比較（壓為負(fù)）。表6給出不同比例跨尺度模型與有限元一致尺寸模型計(jì)算結(jié)果的平均偏差，表7給出了各模型的單元總數(shù)及靜、動力計(jì)算時(shí)間。

結(jié)果表明，無網(wǎng)格界面很容易實(shí)現(xiàn)了面板壩跨尺度分析，并可大幅提高計(jì)算效率。在跨尺度比例為1∶2時(shí)，結(jié)果與有限元一致尺寸網(wǎng)格基本一致，而1∶3時(shí)差計(jì)算結(jié)果產(chǎn)生一定偏差但控制5%左右。因此，在應(yīng)用非點(diǎn)對點(diǎn)界面對面板壩進(jìn)行分析時(shí)，面板與墊層區(qū)網(wǎng)格順坡向尺寸跨越可控制在1∶2到1∶3 之間，在保證精度的前提下，顯著提高計(jì)算效率。如果繼續(xù)考慮在堆石區(qū)內(nèi)進(jìn)行疏密網(wǎng)格過渡，可進(jìn)一步提高計(jì)算效率[10-12]。

4 結(jié)論

通過引入背景網(wǎng)格線及徑向插值函數(shù)（RPIM），開發(fā)了基于無網(wǎng)格的界面模擬方法，實(shí)現(xiàn)了界面兩側(cè)節(jié)點(diǎn)自由分布且無需對應(yīng)，并可靈活地連接不同尺寸面板與堆石體網(wǎng)格，進(jìn)而建立了面板與墊層跨尺度分析模型。（1）基于無網(wǎng)格的非點(diǎn)對點(diǎn)界面，可通過自由分布且無需對應(yīng)的節(jié)點(diǎn)連接不同尺度的面板與墊層區(qū)網(wǎng)格，從而實(shí)現(xiàn)對面板壩跨尺度分析。（2）通過更新無網(wǎng)格界面兩側(cè)的節(jié)點(diǎn)，替換后的堆石體網(wǎng)格可直接與面板原始網(wǎng)格相連，簡化了局部網(wǎng)格精細(xì)化的過程。（3）通過改進(jìn)的徑向插值函數(shù)（RPIM）計(jì)算得到無網(wǎng)格界面的形函數(shù)，在保持無網(wǎng)格界面優(yōu)勢基礎(chǔ)上，提高了邊界附近的插值精度。（4）計(jì)算結(jié)果表明，采用無網(wǎng)格界面對面板壩進(jìn)行分析時(shí)，面板與墊層區(qū)網(wǎng)格順坡向尺度跨越可控制在1∶2到1∶3之間。在保證計(jì)算精度的前提下，大幅度提高計(jì)算效率。

該方法很容易拓展到三維界面的模擬中，效率提高程度將比二維更加明顯，可以為面板壩面板損傷演化精細(xì)化分析提供有效的技術(shù)手段。