一維管單元模型確定裂隙巖體滲透張量的直接優(yōu)化分析

王俊奇，薛振曉，穆孟婧

(華北電力大學(xué)，北京 102206)

裂隙巖體中的滲流影響到水利水電工程、石油工程、土木工程等穩(wěn)定，還涉及核廢料污染的防治、地下結(jié)構(gòu)工程的排水、煤田瓦斯危害的防治、填埋垃圾的污水下滲和地下水資源的開發(fā)等。因此，研究裂隙巖體的滲透性規(guī)律具有非常重要的實(shí)踐意義[1]。目前在巖體滲流的研究中，廣泛采用的模型有：孔隙-裂隙雙重介質(zhì)模型、等效連續(xù)介質(zhì)模型和離散裂隙網(wǎng)絡(luò)模型。離散裂隙網(wǎng)絡(luò)模型由于忽略了巖塊本身的孔隙系統(tǒng)及裂隙巖體之間的水交替過程，用于巖體滲流分析更加簡(jiǎn)便，因此得到迅速發(fā)展[2]。

二維離散裂隙網(wǎng)絡(luò)模型構(gòu)造較為簡(jiǎn)單，發(fā)展成熟，已有學(xué)者在任意裂隙網(wǎng)絡(luò)裂隙-孔隙滲流[3]、二維無壓滲流[4]等方面做出嘗試。而三維模型雖能相對(duì)精確地描述實(shí)際情況，但發(fā)展緩慢，目前在面狀無限延展裂隙、任意裂隙滲流分析方面的研究已有一些成果[5- 6]。

即使離散型裂隙網(wǎng)絡(luò)得到快速發(fā)展，但用以解決工程問題仍有很多困難，在解決實(shí)際工程問題時(shí)滲透張量表征滲透能力的連續(xù)介質(zhì)有限元法仍是常用手段。確定滲透張量一般有物理試驗(yàn)法、理論解析法和數(shù)值方法。數(shù)值法方面，楊天鴻等[7]就二維裂隙網(wǎng)絡(luò)嘗試確定滲透張量；Baghbanan等[8]基于UDEC計(jì)算平臺(tái)，考慮隙寬和跡長(zhǎng)關(guān)聯(lián)性研究二維滲透張量，確定REV。LEE I[9]在TOUGH2/ECO2N模型基礎(chǔ)上，發(fā)展出一種離散裂隙網(wǎng)絡(luò)DFN和非結(jié)構(gòu)網(wǎng)格生成UMG模型，用以分析結(jié)構(gòu)體中CO2的輸運(yùn)；吳錦亮等[10]生成尺寸不同的三維裂隙網(wǎng)絡(luò)，采用復(fù)合單元法計(jì)算各向等效滲透系數(shù)，得到模型巖體的三維滲透張量和REV；李新強(qiáng)等[11]利用邊界元法研究三維裂隙滲透張量；王俊奇等[12]基于一維管單元用遺傳算法反分析嘗試確定三維裂隙網(wǎng)絡(luò)滲透張量。數(shù)值試驗(yàn)法在實(shí)際研究過程中成本較低，操作簡(jiǎn)單，可重復(fù)性高，是目前計(jì)算滲透張量的常用方法。

本文采用數(shù)值模擬方法，建立三維滲流模型，生成圓盤形三維裂隙網(wǎng)絡(luò)，刪除孤立和不連通的裂隙，將圓盤形裂隙網(wǎng)絡(luò)簡(jiǎn)化為三維空間上的管單元模型網(wǎng)絡(luò)。管單元的直徑通過模型校正管徑和隙寬比M的辦法確定，即通過優(yōu)化方法尋找適當(dāng)?shù)墓軉卧睆剑褂?jì)算獲得的滲透張量值能夠滿意地?cái)M合實(shí)測(cè)值。

1 離散裂隙滲流網(wǎng)絡(luò)管單元模型與滲流解法

1.1 三維離散裂隙網(wǎng)絡(luò)及研究域的生成

基于溝槽流假定，三維離散裂隙網(wǎng)絡(luò)滲流發(fā)生在相互交錯(cuò)的裂隙圓盤內(nèi)。根據(jù)裂隙圓盤由七個(gè)特征參數(shù)定義：裂隙圓盤中心的三維坐標(biāo)x、y、z，裂隙圓盤的半徑r，傾向α，傾角β，隙寬b。用蒙特卡羅法生成隨機(jī)裂隙，形成三維結(jié)構(gòu)面網(wǎng)絡(luò)系統(tǒng)。為了保證滲流計(jì)算的合理性及減小計(jì)算規(guī)模，在生成域中選取一研究域立方體，祛除無用裂隙，再用深度優(yōu)先搜索算法求得連通分量，最終確定有效研究域，如圖1所示。

表1 裂隙巖體實(shí)測(cè)數(shù)據(jù)表

圖1 研究域生成示意圖

1.2 一維管單元模型

根據(jù)平面滲流假定，水在等厚度的整個(gè)圓盤平面內(nèi)流動(dòng)，有限元計(jì)算時(shí)需要?jiǎng)澐制矫婢W(wǎng)格，分析較大規(guī)模的實(shí)際問題就不太現(xiàn)實(shí)，研究表明[7]裂隙之間水的流動(dòng)形態(tài)是以溝槽的形式流動(dòng)，據(jù)此可以將三維裂隙網(wǎng)絡(luò)簡(jiǎn)化成一維管單元進(jìn)行分析。根據(jù)文獻(xiàn)[6]建立管單元模型，研究域中相交裂隙產(chǎn)生交線，將交線中點(diǎn)與圓盤圓心相連形成一個(gè)管單元，滲流即表現(xiàn)為管單元通道中的水流。因管單元所屬裂隙面的不同組對(duì)管單元分組，如圖2所示。

圖2 管單元模型

建立了一維管單元模型后，根據(jù)文獻(xiàn)[8]，可以對(duì)滲流區(qū)域內(nèi)所有節(jié)點(diǎn)建立有限元方程組，表現(xiàn)成矩陣形式為：

[K]{H}={q}

(1)

式中，[K]—總滲透矩陣；{H}—節(jié)點(diǎn)水頭的列向量；{q}—節(jié)點(diǎn)流量的列向量。

2 算例

2.1 裂隙網(wǎng)絡(luò)的生成

根據(jù)文獻(xiàn)[13]，查詢獲取實(shí)測(cè)的巖體裂隙數(shù)據(jù)，見表1。

本文自編滲流計(jì)算的有限元程序，輸入已知的實(shí)測(cè)數(shù)據(jù)，生成100m×100m×100m的生成域，在生成域中心選取20m×20m×20m的區(qū)域作為研究域進(jìn)行分析計(jì)算。研究域生成裂隙結(jié)構(gòu)面1978條，其中第1組裂隙的結(jié)構(gòu)面為1226條，第2組為402條，第3組為350條，將這三組不同的裂隙組分別用不同的顏色表示，如圖3所示。將其簡(jiǎn)化為管單元后，研究域中有4104個(gè)單元，3303個(gè)節(jié)點(diǎn)，在研究域邊界沿x方向，y方向，z方向分別施加定水頭H1=30m，H2=10m，其余邊界就為梯度水頭邊界，給定模型邊界條件示意圖如圖4所示(施加ii方向的水頭)，水頭變化趨勢(shì)如圖5所示。

圖3 研究域三組裂隙的結(jié)構(gòu)面生成圖

圖4 計(jì)算模型的邊界條件示意圖

圖5 裂隙網(wǎng)絡(luò)的節(jié)點(diǎn)水頭圖

2.2 滲透張量的計(jì)算

模擬同一個(gè)野外實(shí)際入滲試驗(yàn)的巖體裂隙網(wǎng)絡(luò)的滲流，施加如圖4所示的水頭邊界條件，巖石的管徑隙寬比M=D/b，其中D為管單元直徑，b為裂隙張開寬度。每取一個(gè)M值，可以模擬出一個(gè)滲透張量矩陣：

(2)

設(shè)校準(zhǔn)的的滲透張量矩陣為：

(3)

由于滲透張量矩陣的對(duì)稱性，獨(dú)立的數(shù)據(jù)只有6個(gè)。

設(shè)置目標(biāo)函數(shù)：

(4)

為了得到適當(dāng)?shù)南秾挶仁沟糜?jì)算出滲透張量盡可能地與標(biāo)準(zhǔn)滲透張量吻合，根據(jù)文獻(xiàn)[14]經(jīng)驗(yàn)，先取1～20的M值，求出f的值，將計(jì)算結(jié)果見表2。

表2 裂隙管徑與隙寬之比和目標(biāo)函數(shù)的對(duì)應(yīng)關(guān)系

圖6 M-f擬合曲線圖(M=1～14)

將表2的數(shù)據(jù)繪制成M-f擬合曲線圖可以看出，當(dāng)M的取值在10到12之間的某一個(gè)數(shù)值時(shí)，f取得最小值，為了精準(zhǔn)地求出f的極小值，在10到12之間再插入10個(gè)數(shù)據(jù)求出f的值見表3。

表3 隙寬比為10～12間對(duì)應(yīng)的目標(biāo)函數(shù)值

根據(jù)表3中所提供的數(shù)據(jù)畫出以M為橫坐標(biāo)，以f為縱坐標(biāo)的擬合曲線如圖7所示。

圖7 M-f擬合曲線圖(M=10～12)

得到擬合函數(shù)的趨勢(shì)線為：

f=10-14(M3-40M2+400M-1000)

(5)

對(duì)求導(dǎo)可得：

f'=10-14(3M2-80M+400)

(6)

令f’=0，解得M=11.22，由此可得當(dāng)管徑與隙寬的比值為11.22時(shí)，函數(shù)取得極小值，即模擬得出的滲透系數(shù)與現(xiàn)場(chǎng)測(cè)量得到的滲透張量最為接近。權(quán)且稱此法為直接優(yōu)化的管單元模型，得到的滲透張量矩陣為：

(7)

繼而得到如下相應(yīng)的主滲透系數(shù)和主滲透方向：

(8)

根據(jù)文獻(xiàn)[9]得到現(xiàn)場(chǎng)勘測(cè)的數(shù)據(jù)，采用邊界元方法初步確定主滲透系數(shù)和滲透主方向，接著用現(xiàn)場(chǎng)壓水試驗(yàn)進(jìn)行校核，認(rèn)為具有足夠精度，最終得到以下相對(duì)精確的結(jié)果，本文采用管單元模型計(jì)算得到的滲透張量即是以這個(gè)計(jì)算結(jié)果作為擬合的標(biāo)準(zhǔn)進(jìn)行比對(duì)驗(yàn)證：

(9)

用直接優(yōu)化管單元模型計(jì)算出的結(jié)果式(8)與工程實(shí)例得到的結(jié)果式(9)相比，誤差(ε為相對(duì)誤差，θ為絕對(duì)誤差)如下：

(10)

通過式(10)計(jì)算結(jié)果比較可以看出：求得的滲透主方向與實(shí)測(cè)巖體滲透主方向相差較小，不到30°；而求得的滲透主值與實(shí)測(cè)的滲透主值最大誤差為0.61也在合理的誤差范圍內(nèi)。

文獻(xiàn)[14]具體操作中先做主方向的對(duì)比判斷，再根據(jù)主值進(jìn)行擬合，目標(biāo)函數(shù)取為：

(11)

采用遺傳算法反分析得到的管徑隙寬比為M=11.5069，滲透張量計(jì)算結(jié)果如下：

(12)

(13)

遺傳算法優(yōu)化管單元模型計(jì)算出的結(jié)果式(13)與工程實(shí)例得到的結(jié)果式(9)相比，誤差為：

(14)

3 結(jié)論

(1)一維管單元模型較于三維裂隙滲流網(wǎng)絡(luò)結(jié)構(gòu)更加簡(jiǎn)單，滲流計(jì)算所需的節(jié)點(diǎn)數(shù)和單元數(shù)更少，因而計(jì)算規(guī)模減小，為通過裂隙網(wǎng)絡(luò)滲流確定滲透張量提供可行方法。

(2)通過給定一系列的管徑隙寬比M，用直接優(yōu)化方法計(jì)算出當(dāng)比值M=11.22時(shí)，得到的滲透張量與實(shí)測(cè)的滲透張量相差不大，該比值從統(tǒng)計(jì)意義上認(rèn)為可以用于確定裂隙網(wǎng)絡(luò)巖體滲透張量。

(3)通過與遺傳算法優(yōu)化管單元模型得到的結(jié)果進(jìn)行對(duì)比發(fā)現(xiàn)，在擬合滲透主值方面目標(biāo)函數(shù)優(yōu)化法更精確，在擬合滲透主方向方面遺傳算法誤差更小。