跨尺度網格識別方法及在防滲墻分析中的應用

余翔，賴遠平，王鈺軻，屈永倩，鄭浩然

(1.中國地球物理學會工程物探檢測重點實驗室，湖北 武漢 430000；2.鄭州大學 水利與土木工程學院，河南 鄭州 450001；3.大連理工大學 水利工程學院，遼寧 大連 116024)

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

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

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

本研究結合CAD 人為可控的特點，對圖形線段進行預處理，生成有效節點與線段，根據所生成的節點和線段之間的拓撲關系及封閉域構造，構建可以直接用于數值分析的網格模型，以及二維SBFEM-FEM 耦合網格生成方法.本研究將跨尺度網格識別方法應用于某實際堤壩工程的數值分析研究中，驗證了網格生成方法獲得的SBFEM-FEM 耦合網格的有效性和精度.

1 SBFEM-FEM 耦合分析方法

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

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

對于任意一個扇形區，通過邊界積分點和特征值分解技術可求得單元形函數和應變位移轉換矩陣：

式中： Ψu為位移模式對應的轉換矩陣，可以通過特征值分解求得；Sn為特征值分解獲得的特征值矩陣；n為單元邊數；B(·) 為應變-位移的轉換矩陣；I為單位矩陣.通過采用有限元的原理以及域內積分點，即可求出單元剛度矩陣：

式中：Kep為彈塑性剛度矩陣，Dep為彈塑性本構矩陣，F為每個高斯點所代表的面積，Rint為內力向量，σ 為單元應力場，ε 為應變.

通過網格跨尺度剖分，實現對核心結構網格加密，通常需要用多級過渡，如圖2（a）所示.單元尺寸過渡中必然會引起5 節點及以上節點數目的單元，而這類單元需要單獨構造形函數、積分方法，計算較為麻煩.不需要根據單元的結點數量開發相應單元類型是比例邊界有限元法的一個特點，網格跨尺度通過一級比例邊界有限單元就可以過渡完成，如圖2（b）所示.過渡單元兩側采用有限單元模擬，且采用本研究方法進行網格任意加密，有效地解決了復雜幾何區域的網格剖分.

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

2 SBFEM-FEM 耦合網格生成方法

2.1 CAD 數據處理

本研究提出的SBFEM-FEM 耦合網格生成方法主要分3 步：處理CAD 數據；對節點和線段信息預處理；生成數值網格信息.以下結合圖3 所示的結構，說明耦合網格的生成方法.

圖3 懸臂梁結構的示意圖Fig.3 Schematic diagram of cantilever beam structure

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

圖4 懸臂梁結構的簡化示意圖Fig.4 Simplified schematic diagram of cantilever beam structure

2.2 節點與線段信息預處理

1）創建節點對稱矩陣A，如表1 所示.通過節點關系矩陣A快速確定節點之間的連接關系，相應的節點連接關系，如圖6 所示.若是A(i,j)=1 ，則i節點和j節點可組成ij線段.矩陣A的上三角矩陣中值為1 的個數是整個網格中的線段數量，應與上一步所識別出的有效線段數量相等.

2）創建單個節點連接向量VNL，如表2 所示.VNL(i)用于存儲矩陣A的上三角矩陣第i行中值為1 的個數.

表2 節點連接的數目Tab.2 Number of node connections vector

3）處理共用某一節點的所有線段.創建3 個向量NLN、NLC 與NLA，分別存儲共用節點線段的編號、在相應共用線段中的節點排序（值為1 或2）以及以節點為基準的線段與向量(1,0)的夾角，并根據NLA 由小到大將這些線段進行逆時針排序.

2.3 數值網格信息生成

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

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

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

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

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

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

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

圖8 最終識別出的懸臂梁結構的跨尺度網格Fig.8 Cross-scale mesh of cantilever beam structure of final identification

3 應用案例及分析

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

3.1 計算模型

圖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，壩基、基巖總水平長度290 m.壩體分層填筑和蓄水模擬共分為29 個荷載步：

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

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

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

3.2 材料參數

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

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

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

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

采用Goodman 接觸單元模擬墻與土之間的接觸面，用Clough-Duncan 雙曲線模型[33]表達其應力-應變關系，該模型認為單元切向應力與切向相對位移之間的關系呈雙曲線關系.則切向剛度為

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

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

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

3.3 計算網格模型及驗證

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

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

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

圖12 耦合SBFEM-FEM 網格識別后的單元分布Fig.12 Coupling SBFEM-FEM mesh after mesh recognition

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

表6 節點、單元的信息Tab.6 Information of node and unit

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

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

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

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

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

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

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

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

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

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

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

4 結 論

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

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