跨尺度網(wǎng)格識(shí)別方法及在防滲墻分析中的應(yīng)用

余翔，賴遠(yuǎn)平，王鈺軻，屈永倩，鄭浩然

(1.中國(guó)地球物理學(xué)會(huì)工程物探檢測(cè)重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430000；2.鄭州大學(xué) 水利與土木工程學(xué)院，河南 鄭州 450001；3.大連理工大學(xué) 水利工程學(xué)院，遼寧 大連 116024)

土石堤壩工程是我國(guó)重要的基礎(chǔ)設(shè)施工程，為了保證設(shè)計(jì)的可靠性以及工程的安全性，常常會(huì)對(duì)工程結(jié)構(gòu)的工作狀態(tài)進(jìn)行分析.網(wǎng)格剖分是有限元數(shù)值模擬的基礎(chǔ)，因此網(wǎng)格尺寸大小對(duì)數(shù)值計(jì)算結(jié)果的可靠性和精確性影響很大[1].防滲墻作為控制堤壩滲漏的關(guān)鍵結(jié)構(gòu)，是堤壩安全的控制部位.為了準(zhǔn)確定位在不同運(yùn)行條件下，防滲墻的薄弱部位會(huì)對(duì)防滲墻進(jìn)行精細(xì)的網(wǎng)格剖分.防滲墻與堤體尺度跨越巨大，若按照均勻尺度的原則對(duì)堤壩工程進(jìn)行網(wǎng)格剖分，防滲墻的小尺度網(wǎng)格尺寸標(biāo)準(zhǔn)將會(huì)造成數(shù)值模型具有大量自由度，計(jì)算量難以承受.相比于防滲墻，堤體的應(yīng)力梯度一般較為平緩，常規(guī)尺度網(wǎng)格能夠滿足其精度要求.數(shù)值模擬對(duì)防滲墻的核心結(jié)構(gòu)進(jìn)行精細(xì)網(wǎng)格剖分，其它部分采用常規(guī)尺度網(wǎng)格剖分，獲得了跨尺度網(wǎng)格，并引入數(shù)值計(jì)算方法，數(shù)值分析的計(jì)算量和精度均能夠得到可靠保證.數(shù)值分析的計(jì)算量和精度同時(shí)得到保證是目前數(shù)值分析進(jìn)一步推廣應(yīng)用的難題.

根據(jù)比例邊界有限元法（scaled boundary finite element method, SBFEM）獲得的網(wǎng)格單元對(duì)復(fù)雜幾何邊界適應(yīng)能力強(qiáng)，具有實(shí)現(xiàn)跨尺度、精細(xì)化離散的特點(diǎn).鄒德高等[2-10]對(duì)比例邊界有限元法進(jìn)行發(fā)展和改進(jìn)，取得了一定的成果.結(jié)合比例邊界有限元法對(duì)構(gòu)成單元節(jié)點(diǎn)沒有限制的特點(diǎn)以及有限元法常規(guī)單元技術(shù)成熟易于求解的特點(diǎn)，SBFEM-FEM 耦合求解分析逐漸發(fā)展起來.鄒德高等[11-13]將該耦合法應(yīng)用于結(jié)構(gòu)應(yīng)力分析，驗(yàn)證了該耦合法的可行性.殷德勝等[14-17]將耦合法應(yīng)用于裂縫擴(kuò)展和損傷破壞.Ye 等[18]將該耦合法應(yīng)用于彈性板結(jié)構(gòu)與多層無界彈性地基相互作用問題的求解.目前SBFEM-FEM 跨尺度分析所采用的數(shù)值網(wǎng)格通常由四分樹法獲得.四分樹法獲得的網(wǎng)格單元形狀固定、可控性差，無法適應(yīng)土石堤壩中分層填筑和材料分區(qū)復(fù)雜的邊界條件，難以根據(jù)實(shí)際情況相應(yīng)做出合理的網(wǎng)格剖分方式，且可能造成需要精細(xì)模擬的部位會(huì)被忽略[19].

為了增強(qiáng)網(wǎng)格剖分的可控性，在CAD 軟件中進(jìn)行人為可控的網(wǎng)格劃分，再識(shí)別并生成數(shù)值網(wǎng)格.李華等[20-22]提出二維有限元網(wǎng)格的識(shí)別和生成方法.田林等[23-27]在CAD 環(huán)境下，對(duì)二維網(wǎng)格生成方法進(jìn)行改善優(yōu)化，提高了網(wǎng)格的質(zhì)量和生成效率.籍冉冉等[28]對(duì)莫爾斯（Morse）算法進(jìn)行優(yōu)化，提高了所生成的四邊形網(wǎng)格質(zhì)量.目前均是針對(duì)三角形和四邊形單元的有限元數(shù)值網(wǎng)格，而有關(guān)非常規(guī)的多邊形單元或SBFEM-FEM 耦合數(shù)值網(wǎng)格識(shí)別與生成的研究十分少見.

本研究結(jié)合CAD 人為可控的特點(diǎn)，對(duì)圖形線段進(jìn)行預(yù)處理，生成有效節(jié)點(diǎn)與線段，根據(jù)所生成的節(jié)點(diǎn)和線段之間的拓?fù)潢P(guān)系及封閉域構(gòu)造，構(gòu)建可以直接用于數(shù)值分析的網(wǎng)格模型，以及二維SBFEM-FEM 耦合網(wǎng)格生成方法.本研究將跨尺度網(wǎng)格識(shí)別方法應(yīng)用于某實(shí)際堤壩工程的數(shù)值分析研究中，驗(yàn)證了網(wǎng)格生成方法獲得的SBFEM-FEM 耦合網(wǎng)格的有效性和精度.

1 SBFEM-FEM 耦合分析方法

比例邊界有限元法是由Wolf 等[29-30]提出，該方法在邊界處利用Galerkin 方法進(jìn)行數(shù)值離散，降低了一個(gè)空間的計(jì)算維度，可以在邊界內(nèi)部不需要基本了解的情況下利用解析的方法求解.比例邊界有限元法不僅具有較高的精度，而且對(duì)單元形狀的容許度更高，不需要根據(jù)單元的結(jié)點(diǎn)數(shù)量開發(fā)相應(yīng)單元類型，突破了常規(guī)有限元法的限制.圖1 為典型的多邊形比例邊界單元，圖中ξ 為單元徑向坐標(biāo)，s為環(huán)向坐標(biāo).

圖1 比例邊界有限元法中的邊界離散單元Fig.1 Boundary discrete element in scaled boundary finite element method

對(duì)于任意一個(gè)扇形區(qū)，通過邊界積分點(diǎn)和特征值分解技術(shù)可求得單元形函數(shù)和應(yīng)變位移轉(zhuǎn)換矩陣：

式中： Ψu為位移模式對(duì)應(yīng)的轉(zhuǎn)換矩陣，可以通過特征值分解求得；Sn為特征值分解獲得的特征值矩陣；n為單元邊數(shù)；B(·) 為應(yīng)變-位移的轉(zhuǎn)換矩陣；I為單位矩陣.通過采用有限元的原理以及域內(nèi)積分點(diǎn)，即可求出單元?jiǎng)偠染仃嚕?/p>

式中：Kep為彈塑性剛度矩陣，Dep為彈塑性本構(gòu)矩陣，F(xiàn)為每個(gè)高斯點(diǎn)所代表的面積，Rint為內(nèi)力向量，σ 為單元應(yīng)力場(chǎng)，ε 為應(yīng)變.

通過網(wǎng)格跨尺度剖分，實(shí)現(xiàn)對(duì)核心結(jié)構(gòu)網(wǎng)格加密，通常需要用多級(jí)過渡，如圖2（a）所示.單元尺寸過渡中必然會(huì)引起5 節(jié)點(diǎn)及以上節(jié)點(diǎn)數(shù)目的單元，而這類單元需要單獨(dú)構(gòu)造形函數(shù)、積分方法，計(jì)算較為麻煩.不需要根據(jù)單元的結(jié)點(diǎn)數(shù)量開發(fā)相應(yīng)單元類型是比例邊界有限元法的一個(gè)特點(diǎn)，網(wǎng)格跨尺度通過一級(jí)比例邊界有限單元就可以過渡完成，如圖2（b）所示.過渡單元兩側(cè)采用有限單元模擬，且采用本研究方法進(jìn)行網(wǎng)格任意加密，有效地解決了復(fù)雜幾何區(qū)域的網(wǎng)格剖分.

圖2 跨尺度網(wǎng)格過渡方法的示意圖Fig.2 Schematic diagram of cross-scale mesh transition method

2 SBFEM-FEM 耦合網(wǎng)格生成方法

2.1 CAD 數(shù)據(jù)處理

本研究提出的SBFEM-FEM 耦合網(wǎng)格生成方法主要分3 步：處理CAD 數(shù)據(jù)；對(duì)節(jié)點(diǎn)和線段信息預(yù)處理；生成數(shù)值網(wǎng)格信息.以下結(jié)合圖3 所示的結(jié)構(gòu)，說明耦合網(wǎng)格的生成方法.

圖3 懸臂梁結(jié)構(gòu)的示意圖Fig.3 Schematic diagram of cantilever beam structure

在CAD 軟件中，根據(jù)模擬對(duì)象創(chuàng)建節(jié)點(diǎn)與線段.在CAD 中繪出結(jié)構(gòu)的二維簡(jiǎn)化圖形，如圖4所示.提取圖中6 條線段、12 個(gè)節(jié)點(diǎn)的信息，作進(jìn)一步處理.首先求解線段的交點(diǎn)，并將線段在交點(diǎn)處裁剪為多個(gè)子線段；判斷節(jié)點(diǎn)的公用次數(shù)，若公用次數(shù)不大于1，則該節(jié)點(diǎn)無效，且該節(jié)點(diǎn)所屬線段也無效；最后對(duì)所有有效節(jié)點(diǎn)和線段重新進(jìn)行排序生成基本網(wǎng)格，并存儲(chǔ)節(jié)點(diǎn)和線段信息如圖5 所示.

圖4 懸臂梁結(jié)構(gòu)的簡(jiǎn)化示意圖Fig.4 Simplified schematic diagram of cantilever beam structure

2.2 節(jié)點(diǎn)與線段信息預(yù)處理

1）創(chuàng)建節(jié)點(diǎn)對(duì)稱矩陣A，如表1 所示.通過節(jié)點(diǎn)關(guān)系矩陣A快速確定節(jié)點(diǎn)之間的連接關(guān)系，相應(yīng)的節(jié)點(diǎn)連接關(guān)系，如圖6 所示.若是A(i,j)=1 ，則i節(jié)點(diǎn)和j節(jié)點(diǎn)可組成ij線段.矩陣A的上三角矩陣中值為1 的個(gè)數(shù)是整個(gè)網(wǎng)格中的線段數(shù)量，應(yīng)與上一步所識(shí)別出的有效線段數(shù)量相等.

2）創(chuàng)建單個(gè)節(jié)點(diǎn)連接向量VNL，如表2 所示.VNL(i)用于存儲(chǔ)矩陣A的上三角矩陣第i行中值為1 的個(gè)數(shù).

表2 節(jié)點(diǎn)連接的數(shù)目Tab.2 Number of node connections vector

3）處理共用某一節(jié)點(diǎn)的所有線段.創(chuàng)建3 個(gè)向量NLN、NLC 與NLA，分別存儲(chǔ)共用節(jié)點(diǎn)線段的編號(hào)、在相應(yīng)共用線段中的節(jié)點(diǎn)排序（值為1 或2）以及以節(jié)點(diǎn)為基準(zhǔn)的線段與向量(1,0)的夾角，并根據(jù)NLA 由小到大將這些線段進(jìn)行逆時(shí)針排序.

2.3 數(shù)值網(wǎng)格信息生成

2.3.1 遍歷有效節(jié)點(diǎn) 1）根據(jù)節(jié)點(diǎn)編號(hào)順序，依次判別.當(dāng)節(jié)點(diǎn)編號(hào)為a時(shí)，判別 VNL(a) 大小（VNL(a)為矩陣A的上三角矩陣第a行中值為1 的個(gè)數(shù)），若 VNL(a)≥2 ，則該節(jié)點(diǎn)可生成新單元，否則不可生成，且進(jìn)行下一節(jié)點(diǎn)的判別.以圖7 為例，判別出節(jié)點(diǎn)a可以生成新單元，則進(jìn)入第2.3.2 節(jié)遍歷共用當(dāng)前節(jié)點(diǎn)的線段.2）當(dāng)所有節(jié)點(diǎn)遍歷結(jié)束，則網(wǎng)格識(shí)別完成.

圖7 單元生成的判別方法及流程Fig.7 Identification method and process of unit generation

2.3.2 遍歷共用當(dāng)前節(jié)點(diǎn)的線段 1）識(shí)別共用節(jié)點(diǎn)的線段，并確定過渡點(diǎn).每相鄰的2 條線段(NLN(i+1)與 N LN(i) )為一組，其中必有 N LA(i+1)>NLA(i).確定過渡點(diǎn)后進(jìn)入第2.3.3 步.2）所有相鄰共用線段遍歷結(jié)束，則回至2.3.1 步中第1）步重新判別下一個(gè)節(jié)點(diǎn)的VNL 值大小.

以共用節(jié)點(diǎn)a的線段為例，若單元有n個(gè)節(jié)點(diǎn)，則節(jié)點(diǎn)a為單元的起始點(diǎn)，即第1 個(gè)節(jié)點(diǎn).線段NLN(i+1)為單元的終結(jié)邊，該線段另一節(jié)點(diǎn)b 為單元的終結(jié)點(diǎn)，即第n個(gè)節(jié)點(diǎn)；線段 NLN(i) 為單元的起始邊，該線段另一節(jié)點(diǎn)m為單元的過渡點(diǎn).

2.3.3 判別過渡點(diǎn)與終結(jié)點(diǎn)的連接關(guān)系 1）若過渡點(diǎn)與終結(jié)點(diǎn)可連成線段，即A(i+j)=1 ，則新單元生成，回至第2.3.2 節(jié)判斷共用節(jié)點(diǎn)的下一組線段.2）若A(i+j)≠1 ，則該過渡點(diǎn)為單元中間點(diǎn)，繼續(xù)尋找下一個(gè)過渡點(diǎn).3）順序循環(huán)執(zhí)行本步第1）、2）步.

以上述共用節(jié)點(diǎn)a為例，若是A(b,m)=1 ，則組成單元a-m-b；若是A(b,m)≠1 ，則m組成新單元的中間點(diǎn)，記錄該節(jié)點(diǎn)并繼續(xù)查找組成單元的第n-1 個(gè)節(jié)點(diǎn).查詢共用節(jié)點(diǎn)m的所有線段，確定線段NLN(i)在共用節(jié)點(diǎn)m的所有線段中的排序j，基于組成單元的所有節(jié)點(diǎn)須逆時(shí)針排序的原則，定位線段 NLN(j-1) ，獲得與節(jié)點(diǎn)m共同組成線段NLN(j-1)的另一節(jié)點(diǎn)，該節(jié)點(diǎn)為此時(shí)的過渡節(jié)點(diǎn)m，然后繼續(xù)判別A(b,m) 是否等于1.該例最終生成單元a-m1-m2-m3-b.

用上述方法將圖3 所示結(jié)構(gòu)劃分為用于數(shù)值分析的網(wǎng)格，如圖8 所示.單元尺寸為在CAD 軟件中所畫線段的實(shí)際長(zhǎng)度，通過所提方法的識(shí)別節(jié)點(diǎn)和線段直接確定單元尺寸，材料區(qū)域的數(shù)量對(duì)網(wǎng)格的識(shí)別無任何影響.在網(wǎng)格識(shí)別成功后，根據(jù)材料分區(qū)線、填筑分層線以及界面位置等控制信息自動(dòng)完成材料區(qū)域的劃分，接觸單元生成計(jì)算分析前的關(guān)鍵步驟.

圖8 最終識(shí)別出的懸臂梁結(jié)構(gòu)的跨尺度網(wǎng)格Fig.8 Cross-scale mesh of cantilever beam structure of final identification

3 應(yīng)用案例及分析

某堤壩位于河南省信陽市羅山縣，是以防洪、灌溉為主的?。?）型水庫.該水庫工程等級(jí)為第Ⅳ級(jí)，主要建筑物為4 級(jí).該水庫主壩為均質(zhì)土質(zhì)堤壩，壩基由重粉質(zhì)壤土、細(xì)砂及云母石英片巖組成.壩基土體滲透系數(shù)為1.6×10-6～8.9×10-5cm/s ，平均滲透系數(shù)為 2.6×10-5cm/s ，具有微弱的透水性，在堤壩軸線剖位采用混凝土構(gòu)筑懸掛式防滲墻進(jìn)行防滲.

3.1 計(jì)算模型

圖9 為工程的二維模型，X為順河向，Y為高程方向.壩高15.00 m，壩頂寬5.60 m，壩底寬89.66 m.上游坡比1∶2.73 和1∶2.58，混凝土防滲墻位于壩體的中部，高17 m，嵌入壩基 2 m ，厚度50 cm.壩基高20 m，基巖高2 m，壩基、基巖總水平長(zhǎng)度290 m.壩體分層填筑和蓄水模擬共分為29 個(gè)荷載步：

圖9 堤壩的材料分布及荷載步的示意圖Fig.9 Sketch of material distribution and load steps of embankment

1）分16 步（地基作為初始應(yīng)力）對(duì)壩體從壩基面（高程22 m）填筑至壩頂（高程37 m），墻體材料按土體進(jìn)行計(jì)算；2）第17 步將墻體位置的土體換成混凝土材料進(jìn)行計(jì)算；3）再分12 步從壩基面蓄水至壩頂.其中，蓄水引起的靜水壓力以面力的形式施加于防滲墻的上游面與下游面（嵌入壩基部分），如圖10 所示.

圖10 堤壩防滲系統(tǒng)受力示意圖Fig.10 Schematic diagram of force on anti-seepage system

3.2 材料參數(shù)

堤壩體與堤壩基的靜力分析了采用鄧肯-張E-B 非線性模型[31]，壩體、壩基材料參數(shù)取自實(shí)際試驗(yàn)結(jié)果[32].參數(shù)如表3 所示.其中 γd、 γf分別為干重度和浮重度；c、φ和 ?φ分別為凝聚力、摩擦角和摩擦角隨圍壓增大10 倍的減小量，三者均為強(qiáng)度參數(shù)；k和 β 分別為彈性模量系數(shù)和指數(shù)；Rf為破壞比；kb和q分別為體積模量系數(shù)和指數(shù)參數(shù)；kur和nur分別為卸載模量系數(shù)和指數(shù).

表3 壩體和壩基的靜力參數(shù)Tab.3 Static parameters of dam body and foundation

混凝土防滲墻和基巖材料均采用線彈性本構(gòu)模型計(jì)算，材料參數(shù)如表4 所示，其中 ρ 為密度、E為彈性模量、 ν 為泊松比.

表4 墻、基巖的靜力參數(shù)Tab.4 Static parameters of wall and bedrock

采用Goodman 接觸單元模擬墻與土之間的接觸面，用Clough-Duncan 雙曲線模型[33]表達(dá)其應(yīng)力-應(yīng)變關(guān)系，該模型認(rèn)為單元切向應(yīng)力與切向相對(duì)位移之間的關(guān)系呈雙曲線關(guān)系.則切向剛度為

式中：τ 為切向剛度，ω 為應(yīng)變力，u為實(shí)驗(yàn)參數(shù)，K、Rf和j為非線性參數(shù)，δ 為接觸面的界面摩擦角，γw為水的重度，P為大氣壓，σn為法向應(yīng)力.

董景剛[34]對(duì)接觸面單元的法向剛度應(yīng)取值為相鄰材料彈性模量的50～100 倍，接觸面模型的法向剛度與切向剛度的關(guān)系大致為.法向剛度 = ， 是50 倍的混凝土彈性模量.初始切向剛度= ，是法向剛度的0.4 倍.接觸面材料參數(shù)取自實(shí)際試驗(yàn)結(jié)果[35]參數(shù)如表5 所示.

表5 接觸面的靜力參數(shù)Tab.5 Static parameters of contact surface

3.3 計(jì)算網(wǎng)格模型及驗(yàn)證

為了驗(yàn)證本研究數(shù)值網(wǎng)格生成方法的有效性以及SBFEM-FEM 耦合數(shù)值網(wǎng)格的精度，分別根據(jù)以下3 種情況剖分堤壩模型，并采用本研究建立的方法進(jìn)行數(shù)值網(wǎng)格的生成.堤壩主體耦合網(wǎng)格剖分圖和網(wǎng)格剖分細(xì)圖如圖11 所示.模型1：大尺度常規(guī)單元網(wǎng)格（墻體單元豎向長(zhǎng)度與土體單元保持一致）.模型2：耦合的SBFEM-FEM 單元網(wǎng)格（只在墻體進(jìn)行局部加密，土體單元與模型1 一致）.模型3：精細(xì)尺寸常規(guī)單元網(wǎng)格（以模型2 中墻體單元尺寸為基準(zhǔn)，均勻剖分堤壩）.

圖11 墻體與土體局部單元示意圖Fig.11 Diagram of local elements of wall and soil

圖12 為模型2 耦合SBFEM-FEM 網(wǎng)格通過本研究方法識(shí)別生成的有限元模型局部網(wǎng)格示意圖.其中，常規(guī)單元單元類型號(hào)為1、多邊形單元單元類型號(hào)為44、接觸單元單元類型號(hào)為4 表示.可以看出，圖12 與圖11（b）完全一致，表明了本研究建立的數(shù)值網(wǎng)格生成方法的有效性.

圖12 耦合SBFEM-FEM 網(wǎng)格識(shí)別后的單元分布Fig.12 Coupling SBFEM-FEM mesh after mesh recognition

表6 和圖13 為3 種網(wǎng)格下堤壩節(jié)點(diǎn)、單元信息和墻體單元信息，其中N*為節(jié)點(diǎn)總數(shù)量，N為模型總單元數(shù)量，n*為防滲墻單元數(shù)量，M為模型編號(hào)，H1和H2分別為墻體豎向和橫向長(zhǎng)度.模型2 耦合SBFEM-FEM 網(wǎng)格單元數(shù)量與模型1 相差很小，而墻體單元數(shù)量是模型1 的3.8 倍.模型2墻體單元數(shù)量與模型3 墻體單元數(shù)量相等，而總單元數(shù)量是模型3 的約1/3.

表6 節(jié)點(diǎn)、單元的信息Tab.6 Information of node and unit

圖13 3 種模型下堤壩單元和墻體的單元數(shù)量Fig.13 Element number of dam and wall for three conditions

本軟件生成的計(jì)算模型采用大連理工大學(xué)抗震研究所自主研發(fā)的大型巖土工程靜、動(dòng)力非線性分析軟件GEODYNA 進(jìn)行數(shù)值計(jì)算.圖14 為堤壩蓄水完成后，3 種模型下堤體的水平、豎向位移.可以看出網(wǎng)格對(duì)堤壩的整體變形結(jié)果影響很小，3 種模型下的變形等值線基本吻合.當(dāng)研究堤壩整體變形分布規(guī)律時(shí)，常規(guī)有限元模擬及分析方法也能獲得較高的精度.

圖14 蓄水后堤體的位移分布圖Fig.14 Distribution of dam displacement after impoundment

圖15 為蓄水后墻體的水平位移.圖中 ?S為水平位移，H為墻體豎向高程.從圖15 可以看出，墻體位移分布規(guī)律沒有變化.墻體最大位移出現(xiàn)在墻頂，墻體最大位移為7.3 cm，這是蓄水后上游水壓力及上游堤體浮托力共同作用的結(jié)果.防滲墻位于堤壩中部，其變形完全受堤體控制，墻體變形對(duì)單元尺寸不敏感.當(dāng)研究堤壩防滲墻變形規(guī)律時(shí)，采用單元尺寸較大的大尺度網(wǎng)格模擬.

圖15 蓄水后墻體的水平位移示意圖Fig.15 Horizontal displacement of wall after impoundment

圖16～18 為蓄水后墻體上、下游側(cè)和蓄水后堤壩體（不含防滲墻）最大、最小主應(yīng)力分布曲線.圖中 σ1為 最大主應(yīng)力， σ3為最小主應(yīng)力.從圖16～18可以看出，3 種網(wǎng)格獲得的應(yīng)力分布規(guī)律相同，但大尺度常規(guī)網(wǎng)格的計(jì)算結(jié)果與其它2 種差別較大，本研究耦合SBFEM-FEM 網(wǎng)格的計(jì)算結(jié)果與精細(xì)常規(guī)網(wǎng)格計(jì)算結(jié)果相差甚微.大尺度常規(guī)網(wǎng)格計(jì)算結(jié)果低估了墻體峰值，其計(jì)算結(jié)果不能作為評(píng)價(jià)墻體安全性能的依據(jù).

圖16 蓄水后墻體上游面的最大、最小的主應(yīng)力分布圖Fig.16 Maximum and minimum principal stress distribution of wall in upstream after impoundment

圖17 蓄水后墻體下游面的最大、最小主應(yīng)力分布圖Fig.17 Maximum and minimum principal stress distribution of wall in downstream after impoundment

圖18 蓄水后堤壩的最大、最小主應(yīng)力分布圖Fig.18 Maximum and minimum principal stress distribution of embankment in downstream after impoundment

表7 為蓄水后防滲墻主應(yīng)力的峰值.其中F為峰值，以壓應(yīng)力為正，拉應(yīng)力為負(fù)；D為誤差.與精細(xì)網(wǎng)格相比，大尺度網(wǎng)格獲得的墻體大主應(yīng)力誤差可以達(dá)到29.5%，小主應(yīng)力誤差接近50%，而耦合網(wǎng)格的誤差均小于5%.耦合SBFEMFEM 模擬與分析方法可用少量網(wǎng)格獲得較高精度的結(jié)果.

表7 3 種模型下墻體的主應(yīng)力Tab.7 Principal stress of wall for three conditions

4 結(jié) 論

（1）對(duì)于堤壩工程，可以采用常規(guī)大尺度有限元網(wǎng)格分析堤體整體變形規(guī)律以及堤壩中部防滲墻的變形規(guī)律.當(dāng)研究防滲墻應(yīng)力特性時(shí)，防滲墻須采用精細(xì)網(wǎng)格剖分.與常規(guī)精細(xì)化有限元網(wǎng)格的結(jié)果相比，常規(guī)大尺度有限元網(wǎng)格模擬獲得主應(yīng)力的誤差可以超過48%.

（2）基于SBFEM-FEM 耦合數(shù)值網(wǎng)格的數(shù)值模擬，計(jì)算效率及精度均較高.SBFEM-FEM 耦合數(shù)值模擬分析融合了SBFEM 的靈活通用性以及FEM 求解高效的特點(diǎn)，簡(jiǎn)便地完成了大尺度堤體有限元網(wǎng)格與精細(xì)墻體有限元網(wǎng)格的過渡.相應(yīng)跨尺度數(shù)值網(wǎng)格獲得的墻體應(yīng)力與常規(guī)精細(xì)化有限元網(wǎng)格的結(jié)果相比，誤差不超過5%，而跨尺度數(shù)值網(wǎng)格單元量大幅度降低.