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

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

0 前言

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

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

1 滲流場與溫度場相似對比

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

1.1 控制方程的相似性

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

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

1.2 模型的相似性

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

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

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

1.3 ANSYS單元理論相似

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

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

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

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

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

2.1.1 何謂反演

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

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

α≤k1≤β

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

2.1.2 迭代計算滲流自由面

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

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

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

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

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

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

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

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

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

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

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

3.1 工程概況

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

3.1.1 原始資料

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

表1 實測水頭值統(tǒng)計表

3.1.2 模型參數(shù)

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

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

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

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

表3 ANSYS反演計算的結(jié)果和實測的結(jié)果對比表（單位：m）

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

4 結(jié)論及問題

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

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

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