張順福,丁留謙
(1.中國(guó)水利水電科學(xué)研究院 防汛抗旱減災(zāi)研究所,北京 100038;2.水利部防洪抗旱減災(zāi)中心,北京 100038)
組合網(wǎng)格法在無(wú)壓穩(wěn)定滲流中的應(yīng)用
張順福1,2,丁留謙1,2
(1.中國(guó)水利水電科學(xué)研究院 防汛抗旱減災(zāi)研究所,北京 100038;2.水利部防洪抗旱減災(zāi)中心,北京 100038)
針對(duì)無(wú)壓穩(wěn)定滲流分析過(guò)程中常遇到的線形結(jié)構(gòu)問(wèn)題,提出將組合網(wǎng)格法應(yīng)用于無(wú)壓穩(wěn)定滲流計(jì)算。采用兩套獨(dú)立的網(wǎng)格,對(duì)整體區(qū)域采用尺寸較大的粗網(wǎng)格進(jìn)行模擬,不考慮線形結(jié)構(gòu)的影響,在線形結(jié)構(gòu)附近則采用尺寸比較小的細(xì)網(wǎng)格進(jìn)行模擬,考慮線形結(jié)構(gòu)的影響。計(jì)算在粗細(xì)兩套網(wǎng)格之間迭代調(diào)整,直到達(dá)到收斂精度。組合網(wǎng)格法大大簡(jiǎn)化離散工作量,且適應(yīng)于非規(guī)則網(wǎng)格。數(shù)值計(jì)算結(jié)果表明,組合網(wǎng)格法合理可行,計(jì)算程序可靠,這為線形結(jié)構(gòu)的無(wú)壓穩(wěn)定滲流分析提供了一種思路和方法。
無(wú)壓穩(wěn)定滲流;組合網(wǎng)格;線形結(jié)構(gòu);改進(jìn)初流量法;有限單元法
在飽和無(wú)壓穩(wěn)定滲流有限元數(shù)值模擬中,常常會(huì)遇到這樣一些結(jié)構(gòu),它們?cè)谀骋环较虻某叽缗c整體計(jì)算區(qū)域相比非常小,比如巖體中大量存在的斷層、蝕變帶、軟弱夾層等結(jié)構(gòu)面;土石壩中常用的土工膜、防滲面板、防滲心墻、防滲帷幕等局部構(gòu)造等等(為方便描述,下文中統(tǒng)稱(chēng)這些結(jié)構(gòu)為線形結(jié)構(gòu)或奇異區(qū)域)。當(dāng)這些結(jié)構(gòu)中的水流運(yùn)動(dòng)可以忽略不計(jì)的時(shí)候,可以采用結(jié)點(diǎn)雙編號(hào)[1]或內(nèi)部隔水邊界[2]對(duì)這些結(jié)構(gòu)進(jìn)行處理,否則對(duì)這些結(jié)構(gòu)需要進(jìn)行精細(xì)離散,但這些結(jié)構(gòu)在某一方向的尺寸可能達(dá)到厘米甚至毫米級(jí)別,與整體模擬區(qū)域的尺寸相差太遠(yuǎn),根據(jù)有限元誤差理論,為保持計(jì)算精度其它區(qū)域必然也需要精細(xì)的離散,從而導(dǎo)致離散所需結(jié)點(diǎn)數(shù)目和單元數(shù)目將會(huì)大大增加,且考慮到結(jié)構(gòu)的空間特性,整體計(jì)算區(qū)域的離散難度將會(huì)大幅度提升,計(jì)算所需時(shí)間也大大增加,甚至到達(dá)一般計(jì)算機(jī)難以承受的程度。
對(duì)于這類(lèi)局部特性強(qiáng)烈的線形結(jié)構(gòu),一般是采用非規(guī)則網(wǎng)格局部加密的方法對(duì)這些結(jié)構(gòu)及其附近區(qū)域進(jìn)行局部的離散加密,但局部加密網(wǎng)格法的要求比較高,不僅要求整個(gè)計(jì)算區(qū)域的網(wǎng)格質(zhì)量很好,而且還要求網(wǎng)格的尺寸過(guò)渡非常光滑,這些條件在二維情況下尚且能夠很好的滿(mǎn)足,但在三維情況下,對(duì)復(fù)雜區(qū)域的非規(guī)則網(wǎng)格劃分常常遇到下述問(wèn)題:直接劃分失敗,劃分成功但單元質(zhì)量差,或者劃分成功但單元尺寸過(guò)渡不光滑等。
為避免局部加密所遇到的問(wèn)題,一種可行的方案是采用非擬一致網(wǎng)格(Non-conform ing Grid)技術(shù),非擬一致網(wǎng)格技術(shù)的主要思想是將整體區(qū)域分解為相互之間只在區(qū)域表面有交集的不重疊小區(qū)域,區(qū)域采用獨(dú)立的網(wǎng)格離散,區(qū)域之間不需要滿(mǎn)足網(wǎng)格協(xié)調(diào)的要求,區(qū)域之間的連續(xù)性條件通過(guò)多點(diǎn)約束技術(shù)(Multiple-Point Constraints)實(shí)現(xiàn),通過(guò)運(yùn)用這種局部區(qū)域協(xié)調(diào)、整體區(qū)域非協(xié)調(diào)的網(wǎng)格離散方式,達(dá)到不同的局部區(qū)域上采用不同的網(wǎng)格離散或者不同的數(shù)值計(jì)算方法的目的,也就是說(shuō)可以在某個(gè)區(qū)域采用有限元方法進(jìn)行計(jì)算,而在其它區(qū)域采用有限體積或者有限差分進(jìn)行計(jì)算[3]。
非擬一致網(wǎng)格技術(shù)很好的將局部特性控制在小范圍內(nèi),使局部區(qū)域的特性不至于導(dǎo)致整體計(jì)算的急劇增長(zhǎng)。但有效的適應(yīng)局部特性目標(biāo)往往與計(jì)算求解過(guò)程相沖突,包括:(1)由于離散尺寸的差異變化而導(dǎo)致方程求解器會(huì)降低效率甚至求解失?。唬?)由于考慮了非結(jié)構(gòu)網(wǎng)格,數(shù)據(jù)結(jié)構(gòu)變得繁瑣;(3)計(jì)算機(jī)體系結(jié)構(gòu)可能會(huì)不能夠高效的處理非規(guī)則網(wǎng)格的影響(比如,向量化被抑制);(4)即使是在離散過(guò)程中,這個(gè)目標(biāo)也會(huì)帶來(lái)重重困難,對(duì)于有限差分格式而言,很難為非規(guī)則網(wǎng)格構(gòu)造出相應(yīng)的高精度差分格式;對(duì)于有限元計(jì)算而言,一組質(zhì)量好的非結(jié)構(gòu)化網(wǎng)格的自動(dòng)生成意味著高額的花銷(xiāo)[4-6]。
為避免非擬一致網(wǎng)格的困境,Mc Cormick[6]提出了一種有疊加的區(qū)域分解法——快速自適應(yīng)組合網(wǎng)格(FAC)方法。FAC方法首先形成一套對(duì)應(yīng)于整體區(qū)域的規(guī)則粗網(wǎng)格,然后對(duì)規(guī)則粗網(wǎng)格上相應(yīng)于關(guān)注的局部區(qū)域形成一套或者多套與整體粗網(wǎng)格嵌套的細(xì)網(wǎng)格,計(jì)算的時(shí)候先對(duì)粗網(wǎng)格求解,然后將粗網(wǎng)格的計(jì)算結(jié)果通過(guò)延拓算子延伸到細(xì)網(wǎng)格上,進(jìn)而獲得細(xì)網(wǎng)格上的誤差,將這些誤差再通過(guò)限定算子返回粗網(wǎng)格上,粗細(xì)網(wǎng)格之間循環(huán)迭代,反復(fù)對(duì)誤差進(jìn)行校正直至收斂[7-8]。
由于規(guī)則網(wǎng)的處理(離散與求解)相對(duì)不規(guī)則網(wǎng)而言要容易得多,從而使得快速自適應(yīng)組合網(wǎng)格(FAC)方法成為非常有效的方法,但是快速自適應(yīng)組合網(wǎng)格(FAC)方法高效性有2個(gè)基本要求:(1)“嵌套”和“結(jié)構(gòu)化或規(guī)則的網(wǎng)格”又往往成為其應(yīng)用的限制條件,由于工程實(shí)際問(wèn)題中區(qū)域的復(fù)雜性使得規(guī)則網(wǎng)格和“嵌套”的規(guī)則網(wǎng)格往往很難甚至無(wú)法得到;(2)FAC方法相關(guān)的數(shù)值計(jì)算格式以“嵌套”為基礎(chǔ),如果不滿(mǎn)足“嵌套”條件,那么幾乎所有相關(guān)的計(jì)算都無(wú)從下手。這些均使得FAC的應(yīng)用收到了限制。
為此,本文引入一種改進(jìn)的快速自適應(yīng)組合網(wǎng)格法——組合網(wǎng)格法(Composite Grid Method,CGM)。CGM繼承了FAC反復(fù)在粗細(xì)兩套網(wǎng)格進(jìn)行誤差校正的思想,并且可以使用非規(guī)則網(wǎng)格,細(xì)網(wǎng)格和粗網(wǎng)格相互獨(dú)立,不要求網(wǎng)格嵌套組合網(wǎng)絡(luò)法,該方法已廣泛應(yīng)用于水工結(jié)構(gòu)[9]、油田開(kāi)采數(shù)值模擬[10-11]、電磁場(chǎng)[12]、電磁-機(jī)械耦合[13]、焊接數(shù)值模擬[14-16]等領(lǐng)域。
根據(jù)微分體質(zhì)量守恒定律,在滲流區(qū)域D,三維地下水穩(wěn)定滲流微分方程為[17]:
式中:xi為笛卡爾直角坐標(biāo)系的坐標(biāo)軸;vi為沿坐標(biāo)軸xi方向的滲流速度,m/s;W為源或匯項(xiàng),1/s,以獲得水流為正。
同時(shí)滿(mǎn)足Darcy定律:
式中:kij為滲透系數(shù);H為總水頭,m,;p為孔隙水壓強(qiáng),Pa;γ為水的容重,N/m3。定解條件如下(如圖1所示)。
(1)定水頭邊界Γ1:
(2)流量邊界Γ2(以流入為正):
(3)自由面邊界Γ3:
圖1 滲流計(jì)算邊界條件
(4)滲流逸出邊界Γ4:
式中:H0為已知水頭函數(shù),m;q為已知流量,m3/s;ni為邊界外法線方向余弦。
無(wú)壓穩(wěn)定滲流控制方程屬于橢圓型偏微分方程,本文采用基于子域積分改進(jìn)的初流量法[18]進(jìn)行求解。
3.1 組合網(wǎng)格法基本思想 組合網(wǎng)格法是一種有疊加的區(qū)域分解算法,其基本思想是:采用粗網(wǎng)格和細(xì)網(wǎng)格兩套網(wǎng)格進(jìn)行分析計(jì)算,粗網(wǎng)格對(duì)應(yīng)整個(gè)區(qū)域,不考慮線形結(jié)構(gòu)的影響,細(xì)網(wǎng)格則為包含線形結(jié)構(gòu)在內(nèi)的局部區(qū)域。兩套網(wǎng)格之間反復(fù)迭代直至收斂,網(wǎng)格之間的信息交互通過(guò)網(wǎng)格空間位置插值實(shí)現(xiàn),網(wǎng)格之間無(wú)需嵌套,因此適應(yīng)性強(qiáng),規(guī)則網(wǎng)格和非規(guī)則網(wǎng)格均可使用[19]。組合網(wǎng)格法是文獻(xiàn)[20-21]所提出的方法的改進(jìn)從數(shù)學(xué)上證明了用兩套網(wǎng)格所求得的數(shù)值解可以在粗網(wǎng)格區(qū)域達(dá)到粗網(wǎng)格的精度,在細(xì)網(wǎng)格區(qū)域達(dá)到細(xì)網(wǎng)格的精度,從理論上證明了組合網(wǎng)格法的有效性。
圖2 計(jì)算區(qū)域
式中:u、v分別為細(xì)網(wǎng)格區(qū)域未知量及細(xì)網(wǎng)格以外區(qū)域的水頭值;為相應(yīng)區(qū)域的偏微分算子;為初始流速[22],Wp為相應(yīng)區(qū)域的源或匯項(xiàng);Γ=Ωc?Ωf為兩個(gè)區(qū)域的交界的邊界。
采用尺寸比較大的粗網(wǎng)格TH對(duì)整體區(qū)域Ω進(jìn)行網(wǎng)格剖分,得到有限元離散空間SH;采用尺寸比較小的細(xì)網(wǎng)格Th對(duì)區(qū)域Ωf也進(jìn)行網(wǎng)格剖分,得到有限元離散空間Sh。
根據(jù)虛位移原理,式(7)對(duì)應(yīng)的虛功方程為:
推導(dǎo)得到邊界Γ上滿(mǎn)足[19]:
3.3 程序編制 組合網(wǎng)格算法在三維穩(wěn)定無(wú)壓滲流計(jì)算中實(shí)現(xiàn)的具體流程如圖3[25]所示,其中a場(chǎng)代表粗網(wǎng)格區(qū)域,b場(chǎng)代表采用細(xì)網(wǎng)格控制參數(shù)的細(xì)網(wǎng)格區(qū)域,c場(chǎng)代表采用粗網(wǎng)格控制參數(shù)的細(xì)網(wǎng)格區(qū)域,Uf、Uc分別為b場(chǎng)和c場(chǎng)的水頭值。據(jù)此編制了相應(yīng)的三維有限元程序IIFM3DS-CGM[25]。
某無(wú)限長(zhǎng)土壩寬100.5m、高30.0m,中部有0.5 m的垂直心墻,心墻距離上游50.0 m,如圖4所示。心墻和壩殼均為各向同性材料,上游水位為30.0m,下游水位為6.0m。壩長(zhǎng)方向取單位長(zhǎng)度,采用組合網(wǎng)格法,用細(xì)網(wǎng)格模擬心墻及其上、下游17.0 m區(qū)域范圍,整體區(qū)域則用粗網(wǎng)格進(jìn)行離散,其中細(xì)網(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。為驗(yàn)證組合網(wǎng)格法的精度,文中同時(shí)采用較細(xì)的網(wǎng)格對(duì)整體區(qū)域進(jìn)行離散,整體細(xì)網(wǎng)格在心墻及其上、下游17 m的范圍內(nèi)采用與圖6一樣的網(wǎng)格剖分,其它區(qū)域用尺寸為1m×1m×0.25m的單元進(jìn)行離散。
表1給出了壩殼滲透系數(shù)k1與心墻滲透系數(shù)k2不同比值下,組合網(wǎng)格法計(jì)算得到的心墻下游出水
圖3 計(jì)算流程
圖4 計(jì)算模型
圖5 不考慮心墻的整體區(qū)域及相應(yīng)粗網(wǎng)格剖分(組合網(wǎng)格法)
圖6 考慮心墻的局部區(qū)域及相應(yīng)細(xì)網(wǎng)格剖分(組合網(wǎng)格法)
表1 心墻下游出水點(diǎn)高程比較
點(diǎn)高程與整體區(qū)域細(xì)網(wǎng)格的計(jì)算結(jié)果比較,從表1可以看到,組合網(wǎng)格法計(jì)算的結(jié)果與整體區(qū)域細(xì)網(wǎng)格的計(jì)算結(jié)果基本一致,最大誤差為2.3%。圖7給出了不同滲透系數(shù)比值條件下,組合網(wǎng)格法計(jì)算得到的心墻附近自由面與整體區(qū)域細(xì)網(wǎng)格計(jì)算結(jié)果的比較。從圖7可以看到,組合網(wǎng)格法的計(jì)算結(jié)果與整體區(qū)域細(xì)網(wǎng)格的結(jié)果非常接近。
無(wú)壓穩(wěn)定滲流有限元分析中,線形結(jié)構(gòu)增加了計(jì)算區(qū)域的離散難度。本文介紹了組合網(wǎng)格法的基本原理,并將其引入到穩(wěn)定無(wú)壓滲流有限元分析中,為線形結(jié)構(gòu)的模擬提供了一種新思路。
圖7 組合網(wǎng)格與整體細(xì)網(wǎng)格的自由面比較
由于組合網(wǎng)格的粗、細(xì)網(wǎng)格可以獨(dú)立生成,因此組合網(wǎng)格法對(duì)線形結(jié)構(gòu)的網(wǎng)格剖分具有較強(qiáng)的適應(yīng)性,能極大的簡(jiǎn)化無(wú)壓滲流計(jì)算的網(wǎng)格剖分工作。組合網(wǎng)格法在粗網(wǎng)格區(qū)域達(dá)到粗網(wǎng)格的精度,在細(xì)網(wǎng)格區(qū)域達(dá)到細(xì)網(wǎng)格的精度。
結(jié)合改進(jìn)的初流量法給出了穩(wěn)定無(wú)壓滲流組合網(wǎng)格算法格式和計(jì)算流程,并編制了相應(yīng)的三維有限元程序IIFM3DS-CGM。采用垂直心墻算例對(duì)算法和程序進(jìn)行了驗(yàn)證,結(jié)果證明算法是合理的,程序是有效的、可信的。
本文的研究成果還需結(jié)合實(shí)際工程進(jìn)行更深入的研究。
[1] 王浩然,朱國(guó)榮 .淄博市孝婦河源區(qū)地下水資源的開(kāi)發(fā)利用研究[J].高校地質(zhì)學(xué)報(bào),2000,4(2):198-204.
[2] 薛禹群,謝春紅.水文地質(zhì)學(xué)的數(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] 呂濤,石濟(jì)民,林振寶.區(qū)域分解算法[M].北京:科學(xué)出版社,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] 江思珉,朱國(guó)榮,王浩然,等 .FAC方法在地下水?dāng)?shù)值模擬中的應(yīng)用[J].水利學(xué)報(bào),2006,37(11):1389-1392.
[8] 陶文銓.數(shù)值傳熱學(xué)[M].第2版.西安:西安交通大學(xué)出版社,2001.
[9] 郭勝山,李德玉,唐菊珍.水工結(jié)構(gòu)動(dòng)接觸問(wèn)題的組合網(wǎng)格算法[J].水力發(fā)電,2009(5):24-26.
[10] 韓修廷,梁國(guó)平,吳恩成,等.油田套管損壞數(shù)值模擬專(zhuān)用軟件開(kāi)發(fā)及應(yīng)用[J].數(shù)字石油和化工,2006(4):33-36.
[11] 王伯軍,張士誠(chéng),張勁,等.考慮注采關(guān)系的三維地應(yīng)力場(chǎng)數(shù)值模擬研究[J].西南石油大學(xué)學(xué)報(bào):自然科學(xué)版,2006(6):33-35.
[12] 甘艷,阮江軍,張宇,等.組合網(wǎng)格法及其在電磁問(wèn)題中的應(yīng)用[J].電工技術(shù)學(xué)報(bào),2008,23(8):8-14.
[13] 張宇,阮江軍,劉兵,等.組合網(wǎng)格法在電磁-機(jī)械耦合問(wèn)題中的應(yīng)用[J].中國(guó)電機(jī)工程學(xué)報(bào),2007,27(37):42-47.
[14] 陳文平.組合網(wǎng)格法及其在焊接數(shù)值模擬中的應(yīng)用[D].福州:福建師范大學(xué),2009.
[15] 陳文平,馬昌鳳,蔣利華.組合網(wǎng)格法在攪拌摩擦焊接數(shù)值模擬中的應(yīng)用[J].安徽理工大學(xué)學(xué)報(bào):自然科學(xué)版,2009(2):57-61.
[16] 陳文平,馬昌鳳,蔣利華.激光焊接問(wèn)題的組合網(wǎng)格法[J].福建師范大學(xué)學(xué)報(bào):自然科學(xué)版,2009(3):29-32
[17] 王洪濤.多孔介質(zhì)污染物遷移動(dòng)力學(xué)[M].北京:高等教育出版社,2008.
[18] 張順福,丁留謙,劉昌軍,等.基于子域積分的初流量法改進(jìn)[J].水電能源科學(xué),2013,31(3):15-18.
[19] 王德生.組合網(wǎng)格法和非結(jié)構(gòu)化網(wǎng)格自動(dòng)生成[D].北京:中國(guó)科學(xué)院數(shù)學(xué)研究所,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].水利學(xué)報(bào),1988(8):18-26.
[23] 李新強(qiáng).無(wú)壓滲流有限元分析的改進(jìn)初流量法[J].水利學(xué)報(bào),2007,38(8):961-965.
[24] 王媛.求解有自由面滲流問(wèn)題的初流量法的改進(jìn)[J].水利學(xué)報(bào),1998(3):68-74.
[25] 張順福.子域積分法及組合網(wǎng)格法在無(wú)壓滲流中的應(yīng)用[D].北京:中國(guó)水利水電科學(xué)研究院,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
(責(zé)任編輯:李 琳)
2015-08-04
中國(guó)科技部國(guó)際合作項(xiàng)目(2010DFA74520)
張順福(1982-),男,廣西人,博士,工程師,主要從事滲流與控制研究。E-mail:zhangshunfu@live.cn