蔣 進(jìn) ，關(guān)勁松

（1.重慶交通大學(xué)，重慶 400042；2.遼源市東遼縣公路管理段，吉林 遼源 136200）

0 前言

巖土滲流有關(guān)參數(shù)的正確取用是滲流計(jì)算分析中的重要環(huán)節(jié)。數(shù)值計(jì)算或模擬試驗(yàn)的精度再高，若選用參數(shù)不當(dāng)也將導(dǎo)致計(jì)算分析結(jié)果的不可靠，而使確定的滲流控制方案不安全或過(guò)于浪費(fèi)。目前計(jì)算機(jī)的精度已高達(dá)小數(shù)點(diǎn)后若干位，而參數(shù)選用中的誤差卻會(huì)高達(dá)數(shù)量級(jí)，這種計(jì)算中的精度失調(diào)問(wèn)題實(shí)在應(yīng)該引起注意。然而由于介質(zhì)的復(fù)雜性，無(wú)論是實(shí)驗(yàn)室試驗(yàn)還是現(xiàn)場(chǎng)試驗(yàn)，不僅需要花費(fèi)大量的人力、物力和時(shí)間，而且往往得不到滿意的結(jié)果，造成大量的安全隱患。隨著數(shù)值解法和計(jì)算機(jī)的發(fā)展，使利用數(shù)學(xué)模型反求參數(shù)成為可能。

本文將針對(duì)這一事實(shí)專門闡述如何利用大型通用有限元計(jì)算軟件ANSYS的優(yōu)化設(shè)計(jì)功能反饋分析滲透系數(shù)，使人們能夠根據(jù)已有經(jīng)驗(yàn)或研究成果對(duì)地質(zhì)勘探資料和建筑材料性能做出比較可靠的估計(jì)，或者對(duì)取樣利用室內(nèi)試驗(yàn)設(shè)備測(cè)定及野外現(xiàn)場(chǎng)試驗(yàn)測(cè)定取得的參數(shù)的準(zhǔn)確性做出評(píng)價(jià)和修正以得到更可靠的有關(guān)參數(shù)。

1 滲流場(chǎng)與溫度場(chǎng)相似對(duì)比

滲流場(chǎng)中存在水頭梯度,水緩慢地從高水頭的地方流向低水頭的地方；同樣溫度場(chǎng)中的熱能在溫度梯度作用下從溫度高的地方向溫度低的地方傳導(dǎo)。 由此看來(lái)溫度場(chǎng)和滲流場(chǎng)的傳播過(guò)程具有較好的相似性，同時(shí)滲流與熱傳導(dǎo)的微分方程形式、控制變量，以及相關(guān)參數(shù)也具有很好的相似性，因此利用溫度場(chǎng)來(lái)模擬滲流場(chǎng)在理論上是具有可行性的。

1.1 控制方程的相似性

對(duì)于滲流，假定土體和水體是不可壓縮且介質(zhì)均勻各項(xiàng)同性，則應(yīng)用能力守恒定律可推導(dǎo)其控制方程為：

當(dāng)溫度場(chǎng)中無(wú)熱源時(shí)，qgen=0,此時(shí)兩式各項(xiàng)相似，只是其物理意義有改變，T表示溫度而H表示水頭，k表示熱傳導(dǎo)系數(shù)和滲透系數(shù)。

1.2 模型的相似性

溫度場(chǎng)模型和滲流場(chǎng)模型相比擬，還應(yīng)滿足幾何相似和邊界條件一致。

（1）幾何相似。溫度場(chǎng)模型的外部邊界應(yīng)和所研究滲流區(qū)域的外部邊界在幾何上相似。當(dāng)滲流區(qū)域?yàn)榫|(zhì)巖層時(shí)，則模型也應(yīng)是均質(zhì)的；當(dāng)滲流區(qū)域是非均質(zhì)巖層時(shí),則要求模型中不同導(dǎo)熱介質(zhì)的分界線應(yīng)與非均質(zhì)巖層的分界線亦保持相似，此時(shí)應(yīng)注意滲流的微分控制方程有所改變。

（2）邊界條件一致。即溫度場(chǎng)模型的絕熱邊界與滲流區(qū)域的隔水邊界相對(duì)應(yīng)。導(dǎo)熱邊界和透水邊界相對(duì)應(yīng)，導(dǎo)熱邊界上的溫度則和透水邊界上的水頭相對(duì)應(yīng)。

1.3 ANSYS單元理論相似

對(duì)于二維四節(jié)點(diǎn)矩形有限元單元，同1.2節(jié)中的分析，可以推導(dǎo)出其傳導(dǎo)矩陣為：

其中 C1=kx，C2=ky，C3=0，k 表示滲透系數(shù)。

對(duì)比滲流場(chǎng)和溫度場(chǎng)不難發(fā)現(xiàn)，其二維四節(jié)點(diǎn)矩形有限元單元的傳導(dǎo)矩陣的推到過(guò)程和結(jié)果都相同，這就更進(jìn)一步表明，利用ansys溫度場(chǎng)，模擬滲流場(chǎng)是可行的。

2 滲透系數(shù)反演的實(shí)現(xiàn)方法

2.1 反演在ANSYS中的實(shí)現(xiàn)

2.1.1 何謂反演

所謂的反演就是預(yù)先根據(jù)經(jīng)驗(yàn)估計(jì)或現(xiàn)有的勘測(cè)結(jié)果確定各巖層的滲透系數(shù)，然后應(yīng)用ANSYS求解滲透場(chǎng)的水頭分布，最后將求得的水頭值與相應(yīng)的觀測(cè)值進(jìn)行比較，并根據(jù)其差值和地質(zhì)鉆孔資料來(lái)修改各巖層的輪廓和滲透系數(shù)，如此迭代，反復(fù)求解直至滿足精度要求為止。需特別說(shuō)明的是，這種反求參數(shù)的方法并不能代替地質(zhì)勘探工作，而僅僅只是水文地質(zhì)工作的一個(gè)補(bǔ)充延續(xù)而已。其目的是為了減少人力物力的投入并且得到更準(zhǔn)確的參數(shù)。

最后求得的自由面上的水頭值與相應(yīng)的觀測(cè)值之間應(yīng)滿足如下關(guān)系[2]：

α≤k1≤β

式中，函數(shù)E稱為目標(biāo)函數(shù)，ω為權(quán)因子，一般精度要求越高，相應(yīng)的權(quán)因子值越大，H和H′，分別為節(jié)點(diǎn)水頭的求解值和觀測(cè)值，α、β分別為k的上限和下限。根據(jù)能量守恒定律可知，當(dāng)各巖層的滲透系數(shù)滿足一定的比例條件時(shí)，線性方程組 [K][H]=[f]非唯一解，為消除這個(gè)不定解的因素，預(yù)先精確地給定某一巖層的滲透系數(shù)或由觀測(cè)資料給出滲流區(qū)域的滲流量。

2.1.2 迭代計(jì)算滲流自由面

滲流自由事先是不確定的，確定其位置是滲流計(jì)算的主要內(nèi)容。滲流自由面上的水頭壓力等于大氣壓力，該面上任意一點(diǎn)水頭h等于該點(diǎn)的位置高程。為保證式（3）的解存在唯一解，在滲流自由面上應(yīng)滿足條件：

對(duì)非穩(wěn)定滲流自由面除應(yīng)滿足式（4）外，還應(yīng)滿足第二類邊界條件的流量補(bǔ)給關(guān)系，即按下式計(jì)算滲流自由面降下時(shí)由滲流自由面流入介質(zhì)內(nèi)的單寬流量q：

式中，μ為滲流自由面變動(dòng)范圍內(nèi)的介質(zhì)有效孔隙率或給水度；θ為滲流自由面外法向與鉛錘線的夾角。

用有限單元法計(jì)算一個(gè)物理問(wèn)題，它要求計(jì)算區(qū)域的邊界必須是完全確定的。然而隧道圍巖滲流是有自由面的無(wú)壓滲流，滲流自由面位置恰恰又是不確定的，滲流自由面位置是工程設(shè)計(jì)十分關(guān)心的問(wèn)題，故需通過(guò)計(jì)算求出。滲流自由面需滿足式（4）的條件，這樣一來(lái)，滲流自由面位置需在計(jì)算中通過(guò)迭代來(lái)確定。迭代的計(jì)算步聚是[2]：

（1）首先根據(jù)滲流概念和經(jīng)驗(yàn)假定一條滲流自由面，以確定有限單元法的計(jì)算區(qū)域。

（2）將假定的滲流自由面作為第二類邊界，對(duì)于穩(wěn)定的滲流由式 [K][H]=[f]計(jì)算出滲流自由面結(jié)點(diǎn)的水頭值h。

（3）比較假定滲流自由面結(jié)點(diǎn)的計(jì)算水頭值和其位置高程z視其是否滿足式（4）的條件。若不滿足，則用計(jì)算水頭值去改變結(jié)點(diǎn)的z坐標(biāo)，形成新的假定滲流自由面，同時(shí)還確定了新的有限單元法計(jì)算的區(qū)域。反復(fù)上述計(jì)算步驟直至滲流自由面上的結(jié)點(diǎn)全部滿足給定的允許精度（|h-z|≤ε）。

在迭代過(guò)程中為了不改變網(wǎng)格的劃分和避免單元在壓縮或伸展的過(guò)程中出現(xiàn)畸形，可以采用丟棄單元（節(jié)點(diǎn)）的方法。即是當(dāng)滲流自由面結(jié)點(diǎn)T的計(jì)算水頭值小于或等于其下結(jié)點(diǎn)的z坐標(biāo)時(shí)，將T結(jié)點(diǎn)丟棄，改S結(jié)點(diǎn)為滲流自由面結(jié)點(diǎn)，形成新的有限單元法汁算區(qū)域，而與T結(jié)點(diǎn)有關(guān)的單元?jiǎng)t不再參加計(jì)算。

2.2 ANSYS優(yōu)化設(shè)計(jì)

ANSYS優(yōu)化設(shè)計(jì)有兩種方法：第一種是通用的函數(shù)逼近優(yōu)化方法[3]。這種方法的本質(zhì)是采用最小二乘逼近，求取一個(gè)函數(shù)面來(lái)擬合解空間，然后再對(duì)函數(shù)面求極值。這種方法是一種普適的優(yōu)化方法，不容易陷入局部極值點(diǎn)，但優(yōu)化精度不是很高。另一種方法是對(duì)第一種方法的改進(jìn)方法，叫做梯度尋優(yōu)。如果說(shuō)第一種方法是C0階、大范圍普適的粗優(yōu)方法，那么第二種就是C1階、局部尋優(yōu)的精優(yōu)方法。一般來(lái)說(shuō)，一個(gè)比較復(fù)雜的問(wèn)題需要同時(shí)采用兩種優(yōu)化方法，先用第一種初步求得優(yōu)解的基本范圍，然后再采用第二種方法對(duì)優(yōu)解的范圍進(jìn)行更精確的確定。但是第二種方法不僅耗時(shí)間，而且還容易陷入局部最小點(diǎn)，因此通常的問(wèn)題利用C0階函數(shù)逼近優(yōu)化就足夠了。ANSYS中優(yōu)化設(shè)計(jì)模塊有三大變量：（1）設(shè)計(jì)變量DV:即自變量。（2）狀態(tài)變量SV：用來(lái)體現(xiàn)優(yōu)化的邊界條件，是設(shè)計(jì)變量的函數(shù)。（3）目標(biāo)變量OBJ：最終的優(yōu)化目的。它必須是設(shè)計(jì)變量的函數(shù)，而且只能求其最小值。概括說(shuō)來(lái)ANSYS中實(shí)現(xiàn)優(yōu)化設(shè)計(jì)的過(guò)程就是：在狀態(tài)變量SV的邊界條件下，在設(shè)計(jì)變量的取值范圍內(nèi)不斷調(diào)整DV的取值組合，以使得目標(biāo)函數(shù)OBJ逼近最小值。

3 滲透系數(shù)反演的工程實(shí)例

3.1 工程概況

為了驗(yàn)證以上數(shù)學(xué)模型在ANSYS數(shù)值模擬環(huán)境中分析的準(zhǔn)確性，以及在實(shí)際工程中實(shí)用性，現(xiàn)以八臺(tái)山隧道為實(shí)例，在ANSYS中進(jìn)行結(jié)構(gòu)分析。隧道處于巖體強(qiáng)風(fēng)化帶、碎屑巖類的構(gòu)造裂隙發(fā)育區(qū)、地表巖溶分布區(qū)為地下水的補(bǔ)給區(qū)，主要接受大氣降水的補(bǔ)給。

3.1.1 原始資料

根據(jù)原始資料，選取某斷面四組實(shí)測(cè)水頭值，其值如表1所列。

表1 實(shí)測(cè)水頭值統(tǒng)計(jì)表

3.1.2 模型參數(shù)

模型計(jì)算范圍選取埋深60m，隧道底部以下32m，軸線左右各60m，故整個(gè)模型為120m×100 m。左邊界水頭為64.8 m，右邊界水頭值為9.4 m（該值為相對(duì)值，相對(duì)于模型的原點(diǎn)）,采用ansys plane 55單元。該模型由于是反求隧道圍巖的滲透系數(shù)，所以沒有模擬襯砌和加固圈，并且由于實(shí)際計(jì)算中對(duì)超過(guò)三維的滲透系數(shù)進(jìn)行反分析的理論還不成熟，而且工程當(dāng)中又有一些滲透系數(shù)比較穩(wěn)定，比較容易確定。例如襯砌的滲透系數(shù)，在實(shí)驗(yàn)室中就可以比較準(zhǔn)確地測(cè)定出來(lái)。顯然這類系數(shù)就沒有反演的必要性。另外為保證反分析解的唯一性，并且由于干燥狀態(tài)下圍巖的滲透系數(shù)對(duì)滲流量的影響很小，跟設(shè)計(jì)時(shí)的出入不大，所以選取其設(shè)計(jì)值，不參與優(yōu)化設(shè)計(jì)，而作為優(yōu)化結(jié)果的控制條件。具體的參數(shù)取值如表2所列，其中干燥狀態(tài)下的圍巖滲透系數(shù)為一個(gè)確定的值，而飽和狀態(tài)下的圍巖滲透系數(shù)為一個(gè)設(shè)計(jì)優(yōu)化區(qū)域。

表2 滲透系數(shù)一覽表

3.2 優(yōu)化設(shè)計(jì)

經(jīng)過(guò)10次優(yōu)化設(shè)計(jì)參數(shù)搜索最終得到10組優(yōu)化設(shè)計(jì)組合，其中帶“*”號(hào)的為系統(tǒng)推薦的最優(yōu)組合，最優(yōu)設(shè)計(jì)又包含了33次迭代計(jì)算，并得最優(yōu)組合*SET10*中涌水量為7.402 3m3/d,圍巖飽和滲透系數(shù)為0.451 79m/d。將ANSYS反演計(jì)算的結(jié)果和實(shí)際測(cè)量的結(jié)果對(duì)比可以得到以下數(shù)據(jù)和結(jié)論（見表3）。

表3 ANSYS反演計(jì)算的結(jié)果和實(shí)測(cè)的結(jié)果對(duì)比表（單位：m）

由表3可以看出，隧道圍巖飽和滲透系數(shù)采用優(yōu)化設(shè)計(jì)最優(yōu)解時(shí)，計(jì)算水頭和實(shí)際勘測(cè)水頭相對(duì)殘差均值為-0.405，最大殘差為-0.873，而標(biāo)準(zhǔn)殘差為0.508，這說(shuō)明計(jì)算水頭值與實(shí)測(cè)水頭值擬合很好，能夠滿足精度要求。另外狀態(tài)變量（涌水量）的計(jì)算值為7.406，而實(shí)測(cè)量為7.42，也表明了所求滲透系數(shù)的合理性。

4 結(jié)論及問(wèn)題

通過(guò)ANSYS優(yōu)化設(shè)計(jì)模塊，成功地反演出隧道圍巖的滲透系數(shù)，可以方便地檢驗(yàn)實(shí)際勘測(cè)的滲透系數(shù)的準(zhǔn)確性，在一定程度上減少了野外實(shí)驗(yàn)在人力、物力方面的消耗。對(duì)同等地質(zhì)條件下的工程提供了一個(gè)參考思路。雖然沒有對(duì)隧道加固圈進(jìn)行分析，不過(guò)同樣提供解決加固圈厚度的有限元計(jì)算方法—優(yōu)化設(shè)計(jì)，以加固圈厚度為設(shè)計(jì)變量，分析加固圈厚度與涌水量的關(guān)系，利用ANSYS溫度場(chǎng)模塊可以較好地實(shí)現(xiàn)對(duì)均勻介質(zhì)滲流場(chǎng)的模擬。不過(guò)實(shí)際工程問(wèn)題中，特別是巖體中是沒有均勻介質(zhì)的，只有對(duì)模型做一定的簡(jiǎn)化，才可以應(yīng)用。該例正是因?yàn)樗淼罃嗝嫣庯L(fēng)化嚴(yán)重、巖層破碎厲害，并且沒有呈有規(guī)律方向的裂縫，所以才簡(jiǎn)化為均勻介質(zhì)。滲透系數(shù)的反演在一定程度上解決了在野外測(cè)定滲透系數(shù)的困難。本文在一定程度上說(shuō)明了，可以利用ANSYS的溫度場(chǎng)求解滲流場(chǎng)的問(wèn)題，不過(guò)由于滲流場(chǎng)介質(zhì)的復(fù)雜性，在應(yīng)用時(shí)還需考慮諸多因素對(duì)剛度矩陣的影響。ANSYS優(yōu)化設(shè)計(jì)模塊的巧妙應(yīng)用對(duì)實(shí)際工程來(lái)說(shuō)，無(wú)論在安全上還是在經(jīng)濟(jì)上都是大有益處的。

[1]王崧,劉麗娟，董春敏，等譯.有限元分析—ANSYS理論與應(yīng)用（第三版）[M].北京：電子工業(yè)出版社，2011：296～316.

[2]毛昶熙.滲流計(jì)算分析與控制 (第二版)[M].北京：中國(guó)水力水電出版社，2003：306～316.