組合網(wǎng)格法在無壓穩(wěn)定滲流中的應用

張順福，丁留謙

（1.中國水利水電科學研究院 防汛抗旱減災研究所，北京 100038；2.水利部防洪抗旱減災中心，北京 100038）

組合網(wǎng)格法在無壓穩(wěn)定滲流中的應用

張順福1，2，丁留謙1，2

（1.中國水利水電科學研究院 防汛抗旱減災研究所，北京 100038；2.水利部防洪抗旱減災中心，北京 100038）

針對無壓穩(wěn)定滲流分析過程中常遇到的線形結構問題，提出將組合網(wǎng)格法應用于無壓穩(wěn)定滲流計算。采用兩套獨立的網(wǎng)格，對整體區(qū)域采用尺寸較大的粗網(wǎng)格進行模擬，不考慮線形結構的影響，在線形結構附近則采用尺寸比較小的細網(wǎng)格進行模擬，考慮線形結構的影響。計算在粗細兩套網(wǎng)格之間迭代調(diào)整，直到達到收斂精度。組合網(wǎng)格法大大簡化離散工作量，且適應于非規(guī)則網(wǎng)格。數(shù)值計算結果表明，組合網(wǎng)格法合理可行，計算程序可靠，這為線形結構的無壓穩(wěn)定滲流分析提供了一種思路和方法。

無壓穩(wěn)定滲流；組合網(wǎng)格；線形結構；改進初流量法；有限單元法

1 研究背景

在飽和無壓穩(wěn)定滲流有限元數(shù)值模擬中，常常會遇到這樣一些結構，它們在某一方向的尺寸與整體計算區(qū)域相比非常小，比如巖體中大量存在的斷層、蝕變帶、軟弱夾層等結構面；土石壩中常用的土工膜、防滲面板、防滲心墻、防滲帷幕等局部構造等等（為方便描述，下文中統(tǒng)稱這些結構為線形結構或奇異區(qū)域）。當這些結構中的水流運動可以忽略不計的時候，可以采用結點雙編號［1］或內(nèi)部隔水邊界［2］對這些結構進行處理，否則對這些結構需要進行精細離散，但這些結構在某一方向的尺寸可能達到厘米甚至毫米級別，與整體模擬區(qū)域的尺寸相差太遠，根據(jù)有限元誤差理論，為保持計算精度其它區(qū)域必然也需要精細的離散，從而導致離散所需結點數(shù)目和單元數(shù)目將會大大增加，且考慮到結構的空間特性，整體計算區(qū)域的離散難度將會大幅度提升，計算所需時間也大大增加，甚至到達一般計算機難以承受的程度。

對于這類局部特性強烈的線形結構，一般是采用非規(guī)則網(wǎng)格局部加密的方法對這些結構及其附近區(qū)域進行局部的離散加密，但局部加密網(wǎng)格法的要求比較高，不僅要求整個計算區(qū)域的網(wǎng)格質(zhì)量很好，而且還要求網(wǎng)格的尺寸過渡非常光滑，這些條件在二維情況下尚且能夠很好的滿足，但在三維情況下，對復雜區(qū)域的非規(guī)則網(wǎng)格劃分常常遇到下述問題：直接劃分失敗，劃分成功但單元質(zhì)量差，或者劃分成功但單元尺寸過渡不光滑等。

為避免局部加密所遇到的問題，一種可行的方案是采用非擬一致網(wǎng)格（Non-conform ing Grid）技術，非擬一致網(wǎng)格技術的主要思想是將整體區(qū)域分解為相互之間只在區(qū)域表面有交集的不重疊小區(qū)域，區(qū)域采用獨立的網(wǎng)格離散，區(qū)域之間不需要滿足網(wǎng)格協(xié)調(diào)的要求，區(qū)域之間的連續(xù)性條件通過多點約束技術（Multiple-Point Constraints）實現(xiàn)，通過運用這種局部區(qū)域協(xié)調(diào)、整體區(qū)域非協(xié)調(diào)的網(wǎng)格離散方式，達到不同的局部區(qū)域上采用不同的網(wǎng)格離散或者不同的數(shù)值計算方法的目的，也就是說可以在某個區(qū)域采用有限元方法進行計算，而在其它區(qū)域采用有限體積或者有限差分進行計算［3］。

非擬一致網(wǎng)格技術很好的將局部特性控制在小范圍內(nèi)，使局部區(qū)域的特性不至于導致整體計算的急劇增長。但有效的適應局部特性目標往往與計算求解過程相沖突，包括：（1）由于離散尺寸的差異變化而導致方程求解器會降低效率甚至求解失?。唬?）由于考慮了非結構網(wǎng)格，數(shù)據(jù)結構變得繁瑣；（3）計算機體系結構可能會不能夠高效的處理非規(guī)則網(wǎng)格的影響（比如，向量化被抑制）；（4）即使是在離散過程中，這個目標也會帶來重重困難，對于有限差分格式而言，很難為非規(guī)則網(wǎng)格構造出相應的高精度差分格式；對于有限元計算而言，一組質(zhì)量好的非結構化網(wǎng)格的自動生成意味著高額的花銷［4-6］。

為避免非擬一致網(wǎng)格的困境，Mc Cormick［6］提出了一種有疊加的區(qū)域分解法——快速自適應組合網(wǎng)格（FAC）方法。FAC方法首先形成一套對應于整體區(qū)域的規(guī)則粗網(wǎng)格，然后對規(guī)則粗網(wǎng)格上相應于關注的局部區(qū)域形成一套或者多套與整體粗網(wǎng)格嵌套的細網(wǎng)格，計算的時候先對粗網(wǎng)格求解，然后將粗網(wǎng)格的計算結果通過延拓算子延伸到細網(wǎng)格上，進而獲得細網(wǎng)格上的誤差，將這些誤差再通過限定算子返回粗網(wǎng)格上，粗細網(wǎng)格之間循環(huán)迭代，反復對誤差進行校正直至收斂［7-8］。

由于規(guī)則網(wǎng)的處理（離散與求解）相對不規(guī)則網(wǎng)而言要容易得多，從而使得快速自適應組合網(wǎng)格（FAC）方法成為非常有效的方法，但是快速自適應組合網(wǎng)格（FAC）方法高效性有2個基本要求：（1）“嵌套”和“結構化或規(guī)則的網(wǎng)格”又往往成為其應用的限制條件，由于工程實際問題中區(qū)域的復雜性使得規(guī)則網(wǎng)格和“嵌套”的規(guī)則網(wǎng)格往往很難甚至無法得到；（2）FAC方法相關的數(shù)值計算格式以“嵌套”為基礎，如果不滿足“嵌套”條件，那么幾乎所有相關的計算都無從下手。這些均使得FAC的應用收到了限制。

為此，本文引入一種改進的快速自適應組合網(wǎng)格法——組合網(wǎng)格法（Composite Grid Method，CGM）。CGM繼承了FAC反復在粗細兩套網(wǎng)格進行誤差校正的思想，并且可以使用非規(guī)則網(wǎng)格，細網(wǎng)格和粗網(wǎng)格相互獨立，不要求網(wǎng)格嵌套組合網(wǎng)絡法，該方法已廣泛應用于水工結構［9］、油田開采數(shù)值模擬［10-11］、電磁場［12］、電磁-機械耦合［13］、焊接數(shù)值模擬［14-16］等領域。

2 無壓穩(wěn)定滲流基本理論

根據(jù)微分體質(zhì)量守恒定律，在滲流區(qū)域D，三維地下水穩(wěn)定滲流微分方程為［17］：

式中：xi為笛卡爾直角坐標系的坐標軸；vi為沿坐標軸xi方向的滲流速度，m/s；W為源或匯項，1/s，以獲得水流為正。

同時滿足Darcy定律：

式中：kij為滲透系數(shù)；H為總水頭，m，；p為孔隙水壓強，Pa；γ為水的容重，N/m3。定解條件如下（如圖1所示）。

（1）定水頭邊界Γ1：

（2）流量邊界Γ2（以流入為正）：

（3）自由面邊界Γ3：

圖1 滲流計算邊界條件

（4）滲流逸出邊界Γ4：

式中：H0為已知水頭函數(shù)，m；q為已知流量，m3/s；ni為邊界外法線方向余弦。

無壓穩(wěn)定滲流控制方程屬于橢圓型偏微分方程，本文采用基于子域積分改進的初流量法［18］進行求解。

3 組合網(wǎng)格法原理

3.1 組合網(wǎng)格法基本思想 組合網(wǎng)格法是一種有疊加的區(qū)域分解算法，其基本思想是：采用粗網(wǎng)格和細網(wǎng)格兩套網(wǎng)格進行分析計算，粗網(wǎng)格對應整個區(qū)域，不考慮線形結構的影響，細網(wǎng)格則為包含線形結構在內(nèi)的局部區(qū)域。兩套網(wǎng)格之間反復迭代直至收斂，網(wǎng)格之間的信息交互通過網(wǎng)格空間位置插值實現(xiàn)，網(wǎng)格之間無需嵌套，因此適應性強，規(guī)則網(wǎng)格和非規(guī)則網(wǎng)格均可使用［19］。組合網(wǎng)格法是文獻［20-21］所提出的方法的改進從數(shù)學上證明了用兩套網(wǎng)格所求得的數(shù)值解可以在粗網(wǎng)格區(qū)域達到粗網(wǎng)格的精度，在細網(wǎng)格區(qū)域達到細網(wǎng)格的精度，從理論上證明了組合網(wǎng)格法的有效性。

圖2 計算區(qū)域

式中：u、v分別為細網(wǎng)格區(qū)域未知量及細網(wǎng)格以外區(qū)域的水頭值；為相應區(qū)域的偏微分算子；為初始流速［22］，Wp為相應區(qū)域的源或匯項；Γ=Ωc?Ωf為兩個區(qū)域的交界的邊界。

采用尺寸比較大的粗網(wǎng)格TH對整體區(qū)域Ω進行網(wǎng)格剖分，得到有限元離散空間SH；采用尺寸比較小的細網(wǎng)格Th對區(qū)域Ωf也進行網(wǎng)格剖分，得到有限元離散空間Sh。

根據(jù)虛位移原理，式（7）對應的虛功方程為：

推導得到邊界Γ上滿足［19］:

3.3 程序編制 組合網(wǎng)格算法在三維穩(wěn)定無壓滲流計算中實現(xiàn)的具體流程如圖3［25］所示，其中a場代表粗網(wǎng)格區(qū)域，b場代表采用細網(wǎng)格控制參數(shù)的細網(wǎng)格區(qū)域，c場代表采用粗網(wǎng)格控制參數(shù)的細網(wǎng)格區(qū)域，Uf、Uc分別為b場和c場的水頭值。據(jù)此編制了相應的三維有限元程序IIFM3DS-CGM［25］。

某無限長土壩寬100.5m、高30.0m，中部有0.5 m的垂直心墻，心墻距離上游50.0 m，如圖4所示。心墻和壩殼均為各向同性材料，上游水位為30.0m，下游水位為6.0m。壩長方向取單位長度，采用組合網(wǎng)格法，用細網(wǎng)格模擬心墻及其上、下游17.0 m區(qū)域范圍，整體區(qū)域則用粗網(wǎng)格進行離散，其中細網(wǎng)格心墻處的六面體單元尺寸為0.125m×1m×0.25m，心墻上、下游2m范圍內(nèi)單元尺寸為0.25m×1m×0.25m，其它區(qū)域單元尺寸為1m×1m×0.25m，如圖5。粗網(wǎng)格不考慮心墻的影響，單元尺寸大小0.995m×1m×1m，如圖6。為驗證組合網(wǎng)格法的精度，文中同時采用較細的網(wǎng)格對整體區(qū)域進行離散，整體細網(wǎng)格在心墻及其上、下游17 m的范圍內(nèi)采用與圖6一樣的網(wǎng)格剖分，其它區(qū)域用尺寸為1m×1m×0.25m的單元進行離散。

表1給出了壩殼滲透系數(shù)k1與心墻滲透系數(shù)k2不同比值下，組合網(wǎng)格法計算得到的心墻下游出水

圖3 計算流程

4 算例

圖4 計算模型

圖5 不考慮心墻的整體區(qū)域及相應粗網(wǎng)格剖分（組合網(wǎng)格法）

圖6 考慮心墻的局部區(qū)域及相應細網(wǎng)格剖分（組合網(wǎng)格法）

表1 心墻下游出水點高程比較

點高程與整體區(qū)域細網(wǎng)格的計算結果比較，從表1可以看到，組合網(wǎng)格法計算的結果與整體區(qū)域細網(wǎng)格的計算結果基本一致，最大誤差為2.3%。圖7給出了不同滲透系數(shù)比值條件下，組合網(wǎng)格法計算得到的心墻附近自由面與整體區(qū)域細網(wǎng)格計算結果的比較。從圖7可以看到，組合網(wǎng)格法的計算結果與整體區(qū)域細網(wǎng)格的結果非常接近。

5 結語

無壓穩(wěn)定滲流有限元分析中，線形結構增加了計算區(qū)域的離散難度。本文介紹了組合網(wǎng)格法的基本原理，并將其引入到穩(wěn)定無壓滲流有限元分析中，為線形結構的模擬提供了一種新思路。

圖7 組合網(wǎng)格與整體細網(wǎng)格的自由面比較

由于組合網(wǎng)格的粗、細網(wǎng)格可以獨立生成，因此組合網(wǎng)格法對線形結構的網(wǎng)格剖分具有較強的適應性，能極大的簡化無壓滲流計算的網(wǎng)格剖分工作。組合網(wǎng)格法在粗網(wǎng)格區(qū)域達到粗網(wǎng)格的精度，在細網(wǎng)格區(qū)域達到細網(wǎng)格的精度。

結合改進的初流量法給出了穩(wěn)定無壓滲流組合網(wǎng)格算法格式和計算流程，并編制了相應的三維有限元程序IIFM3DS-CGM。采用垂直心墻算例對算法和程序進行了驗證，結果證明算法是合理的，程序是有效的、可信的。

本文的研究成果還需結合實際工程進行更深入的研究。

［1］ 王浩然，朱國榮 .淄博市孝婦河源區(qū)地下水資源的開發(fā)利用研究［J］.高校地質(zhì)學報，2000，4（2）：198-204.

［2］ 薛禹群，謝春紅.水文地質(zhì)學的數(shù)值法［M］.北京：煤炭工業(yè)出版社，1980.

［3］ Bernardi C，Maday Y，Patera A.Domain decomposition by the mortar finite element method［J］.Asymptotic and Numerical Methods for PDEs with Critical Parameters，1993，384：269-286.

［4］ Mc Cormick S，Thomas J.The fast adaptive composite grid（FAC）method for elliptic equations［J］.Mathematics of Computation，1986，46：439-456.

［5］ 呂濤，石濟民，林振寶.區(qū)域分解算法［M］.北京：科學出版社，1992.

［6］ Mc Cormick S.Fast adaptive composite grid（FAC）methods：theory for the variational case［J］.Defect Correction Methods Computing Supplementum，1984，5：115-121.

［7］ 江思珉，朱國榮，王浩然，等 .FAC方法在地下水數(shù)值模擬中的應用［J］.水利學報，2006，37（11）：1389-1392.

［8］ 陶文銓.數(shù)值傳熱學［M］.第2版.西安：西安交通大學出版社，2001.

［9］ 郭勝山，李德玉，唐菊珍.水工結構動接觸問題的組合網(wǎng)格算法［J］.水力發(fā)電，2009（5）：24-26.

［10］ 韓修廷，梁國平，吳恩成，等.油田套管損壞數(shù)值模擬專用軟件開發(fā)及應用［J］.數(shù)字石油和化工，2006（4）：33-36.

［11］ 王伯軍，張士誠，張勁，等.考慮注采關系的三維地應力場數(shù)值模擬研究［J］.西南石油大學學報：自然科學版，2006（6）：33-35.

［12］ 甘艷，阮江軍，張宇，等.組合網(wǎng)格法及其在電磁問題中的應用［J］.電工技術學報，2008，23（8）：8-14.

［13］ 張宇，阮江軍，劉兵，等.組合網(wǎng)格法在電磁-機械耦合問題中的應用［J］.中國電機工程學報，2007，27（37）：42-47.

［14］ 陳文平.組合網(wǎng)格法及其在焊接數(shù)值模擬中的應用［D］.福州：福建師范大學，2009.

［15］ 陳文平，馬昌鳳，蔣利華.組合網(wǎng)格法在攪拌摩擦焊接數(shù)值模擬中的應用［J］.安徽理工大學學報：自然科學版，2009（2）：57-61.

［16］ 陳文平，馬昌鳳，蔣利華.激光焊接問題的組合網(wǎng)格法［J］.福建師范大學學報：自然科學版，2009（3）：29-32

［17］ 王洪濤.多孔介質(zhì)污染物遷移動力學［M］.北京：高等教育出版社，2008.

［18］ 張順福，丁留謙，劉昌軍，等.基于子域積分的初流量法改進［J］.水電能源科學，2013，31（3）：15-18.

［19］ 王德生.組合網(wǎng)格法和非結構化網(wǎng)格自動生成［D］.北京：中國科學院數(shù)學研究所，2001.

［20］ Xu J.Iterative Methods by Space Decomposition and Subspace Correction［J］.SIAM Review，1992，34（4）：581-613

［21］ Xu J.The auxiliary space method and optimal multigrid preconditioning techniques for unstructured grids［J］. Computing，1996，56（3）：215-236.

［22］ 張有天，陳平，王鐳.有自由面滲流分析的初流量法［J］.水利學報，1988（8）：18-26.

［23］ 李新強.無壓滲流有限元分析的改進初流量法［J］.水利學報，2007，38（8）：961-965.

［24］ 王媛.求解有自由面滲流問題的初流量法的改進［J］.水利學報，1998（3）：68-74.

［25］ 張順福.子域積分法及組合網(wǎng)格法在無壓滲流中的應用［D］.北京：中國水利水電科學研究院，2013.

Application of composite grid method in unconfined seepage analysis

ZHANG Shunfu1，2，DING Liuqian1，2

（1.China Institute of Water Resources and Hydro Power Research，Beijing 100038，China；2.Research Center on Flood&Drought Disaster Reduction of the Ministry of Water Resources，Beijing 100038，China）

To solve the“l(fā)ine-style structure”problem often encountered in unconfined steady seepage analysis，Composite Grid Method is adopted.Two sets of independent grids are used in numerical simulation. One is the coarse grid of larger size simulating the entire region without consideration of“l(fā)ine-style structure”and the other is the fine grid of relatively small size simulating the domain including the“l(fā)ine-style structure”.Adjustments are carried out iteratively between the coarse grid and fine grid until desired convergence precision is achieved.With Composite Grid Method，the amount of discretization is greatly reduced. The Composite Grid Method can be adapted to non-regular grid.Numerical simulation result shows that the theory is reasonable and the program is reliable.Composite Grid Method provides an approach to solving the“l(fā)ine-style structure”problem in unconfined steady seepage analysis.

unconfined steady seepage；Composite Grid Method；line-style structure；improved Initial Flow method；Finite element method

TV139.14

A

10.13244/j.cnki.jiwhr.2016.01.003

1672-3031（2016）01-0016-07

（責任編輯：李 琳）

2015-08-04

中國科技部國際合作項目（2010DFA74520）

張順福（1982-），男，廣西人，博士，工程師，主要從事滲流與控制研究。E-mail：zhangshunfu@live.cn