土質(zhì)心墻堆石壩非飽和滲流參數(shù)反演

尹 川，吳震宇，王 媛，卞 康

(1.四川大學水力學與山區(qū)河流開發(fā)保護國家重點實驗室 水利水電學院，四川 成都 610065；2.成都市市政工程設(shè)計研究院有限公司，四川 成都 610000)

因具備對地形地質(zhì)條件適應(yīng)性好，抗震性能優(yōu)良，能就地取材和充分利用建筑物開挖料等特點，土質(zhì)心墻堆石壩成為我國西部高壩廣泛采用的壩型之一，如瀑布溝(186m)、長河壩(240m)、雙江口(312m)等。土質(zhì)防滲體的工作性態(tài)是心墻堆石壩長期運行安全監(jiān)控的關(guān)鍵項目之一。水庫水位變化范圍內(nèi)的心墻滲流性態(tài)與防滲土料的非飽和滲流特性密切相關(guān)，在水庫水位快速變化條件下，土質(zhì)心墻可能發(fā)生滲透破壞。

對于防滲土料的非飽和滲流特性，眾多學者開展了相關(guān)研究。Rahardjo[1]通過非飽和數(shù)值模擬對比現(xiàn)場吸力，發(fā)現(xiàn)邊界條件對土體吸力的影響。Khanh[2]通過分析瞬態(tài)非飽和滲流條件下的邊坡穩(wěn)定性，研究降雨誘發(fā)邊坡破壞的潛在機制。胡云進[3]建立了單裂隙非飽和滲流模型，研究了滲透系數(shù)以及初始飽和度對裂隙滲透性的影響。數(shù)值模擬是預(yù)測心墻堆石壩滲流性態(tài)的有力工具，但其預(yù)測精度主要取決于合理的滲流材料模型及參數(shù)?；谠^數(shù)據(jù)進行模型識別和參數(shù)反演，以提高滲流數(shù)值模擬精度是目前工程中常用的方法。Strauss[4]基于貝葉斯理論進行滲透率統(tǒng)計反演，解決多孔介質(zhì)中單相滲流的滲透率病態(tài)問題。岑建[5]以中線式尾礦壩為例，改進加速遺傳算法反演巖土滲透系數(shù)，為分析該工程運行到最終壩頂標高時的滲流場提供依據(jù)。石春池針對糯扎渡高心墻堆石壩，利用粒子群優(yōu)化算法優(yōu)化徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)，基于實測滲壓反演壩體的滲透系數(shù)，并利用非飽和非穩(wěn)定滲流場測試反演結(jié)果。

近年來，為了使得反演結(jié)果更加真實地反映材料的實際滲流參數(shù)，人們不斷地引入優(yōu)化算法，進行各分區(qū)滲透系數(shù)的反演分析，應(yīng)用于工程實踐，且都取得了豐碩的成果。然而他們大多將滲透系數(shù)視為常數(shù)，認為材料的滲透特性不會隨環(huán)境量的波動而改變。實際上，飽和土的滲透系數(shù)的確可以認為是常量，但非飽和土中的滲透系數(shù)是含水量或基質(zhì)吸力的函數(shù)，會隨著水庫水位的變化而處于一個動態(tài)調(diào)整的過程。因此，進行非飽和滲流參數(shù)的反演，更能體現(xiàn)各材料的真實參數(shù)。本文以PG礫石土心墻堆石壩為例，基于滲流原觀數(shù)據(jù)進行滲流參數(shù)的反演，進而分析非飽和非穩(wěn)定滲流對心墻滲流場的影響。

1 非飽和非穩(wěn)定滲流參數(shù)反分析方法

1.1 非飽和非穩(wěn)定滲流基本微分方程

考慮介質(zhì)和土體壓縮性，非穩(wěn)定滲流微分方程式如下：

(1)

式中，kx、ky、kz—x、y、z方向的滲透速率分量；Ss—單位貯存量，Ss=ρg(α+nβ)；α—土體壓縮系數(shù)；β—水體壓縮系數(shù)；h—水頭；n—土體孔隙率；ρ—密度；g—重力加速度。

該式既適用于承壓含水層，也適用于無壓滲流。飽和-非飽和滲流認為滲流幾何區(qū)域全部處于水中，自由面是飽和區(qū)與非飽和區(qū)的分界線，邊界水位變化時，滲流場水頭的變化是通過飽和區(qū)與非飽和區(qū)的流量交換實現(xiàn)。非飽和非穩(wěn)定滲流微分方程如下：

(2)

式中，Q—單位時間內(nèi)在垂直方向從單位體積含水層中流入或流出的水量；mw—土中水的質(zhì)量；γw—水容重；y—位置水頭。

1.2 BC&BC非飽和滲流本構(gòu)模型

土的透水性和土的飽和度密切相關(guān)，人們將土的飽和度Se隨基質(zhì)吸力的變化定義為土水特征曲線，將土的滲透系數(shù)隨基質(zhì)吸力的變化定義為滲透系數(shù)函數(shù)。基質(zhì)吸力指的是土體中的毛細作用，當土體飽和時基質(zhì)吸力為零，滲透系數(shù)Ks為常數(shù)，隨著基質(zhì)吸力的增加，滲透系數(shù)呈非線性減小。BC&BC模型是最為常用的非飽和滲流本構(gòu)，其土水特征曲線和滲流系數(shù)函數(shù)的數(shù)學方程見公式(3)—(4)所示。

(3)

(4)

式中，h—吸力水頭；hb—進氣壓力水頭；λ—土的孔徑分布指數(shù)。

1.3 滲流參數(shù)反演方法

對于非飽和滲流模型的參數(shù)反分析，通過建立響應(yīng)面代理模型，代替數(shù)值模擬中大壩的滲透參數(shù)與有限元計算滲壓效應(yīng)量之間的響應(yīng)關(guān)系，然后基于代表性時刻水位原觀滲壓監(jiān)測數(shù)據(jù)與變形響應(yīng)面模型的關(guān)系，構(gòu)建目標函數(shù)，進行滲流參數(shù)尋優(yōu)。本文采用多目標優(yōu)化求解目標函數(shù)，通過心墻上區(qū)和心墻下區(qū)的兩個目標函數(shù)，考慮了滲壓與自身區(qū)域材料特性的關(guān)系，使目標監(jiān)測點均能很好地逼近實測值?；谧钚《朔ǖ脑?，可構(gòu)建參數(shù)反分析的目標函數(shù)：

(5)

響應(yīng)面法是數(shù)學理論和統(tǒng)計技術(shù)的結(jié)合，用于建模分析受多個變量影響的問題，其目的是優(yōu)化響應(yīng)。通過構(gòu)建響應(yīng)面方程能有效的替代大量復(fù)雜繁瑣的數(shù)值模擬計算，可大幅減少工作量，使目標函數(shù)的求解過程相對簡化。本文以心墻滲流參數(shù)和代表性時刻水位對應(yīng)滲壓主成分作為變量擬合響應(yīng)面方程。采用不含交叉項的3次多項式進行構(gòu)建，其形式如下：

(6)

式中，H(x)—正交試驗樣本對應(yīng)有限元模擬的滲壓；ai、bi、ci、e—響應(yīng)面方程系數(shù)，通過統(tǒng)計回歸擬合得到；xi—第i個反分析參數(shù)的取值；n—待反演的目標參數(shù)個數(shù)。

1.4 基于NSGA-Ⅱ多目標遺傳算法

基本遺傳算法是簡單地認為個體適應(yīng)度高更適合遺傳優(yōu)良的基因，而基于帶精英策略的非支配排序(NSGA-Ⅱ)多目標遺傳算法的評優(yōu)機制則不同。多目標函數(shù)中每個目標函數(shù)均可得到一個適應(yīng)度值，但不能同時最優(yōu)，因此其解呈現(xiàn)為Pareto解集，其中任意兩個解之間沒有絕對優(yōu)劣。

基于NSGA-Ⅱ多目標優(yōu)化的主要步驟如下，具體流程圖如圖1所示。

圖1 NSGA-Ⅱ算法的基本流程

(1)初始化：設(shè)定待優(yōu)化參數(shù)集的搜索范圍，結(jié)合正交設(shè)計與隨機補充的個體形成初始種群。

(2)適應(yīng)度：計算每個新個體的適應(yīng)度，父代個體與后代種群合并以保證優(yōu)良基因不丟失。

(3)非支配排序：遍歷所有個體支配關(guān)系，進行非支配排序。

(4)擁擠距離計算：計算同一層中各個個體的擁擠距離。

(5)停止準則：若達到預(yù)設(shè)迭代數(shù)，則輸出當前非支配前沿，并停止反之繼續(xù)。

(6)選擇準則：選擇非支配序小、擁擠度大的個體。

(7)交叉準則：通過SBX交叉算子遺傳父代基因給子代。

(8)變異準則：為避免局部最優(yōu)和過早收斂，將部分父代基因突變來增加搜索范圍。

(9)迭代循環(huán)：轉(zhuǎn)至第②步。

2 PG礫石土心墻堆石壩非飽和滲流參數(shù)反演

2.1 工程概況

PG水電站是一座以發(fā)電為主，兼有防洪等綜合利用效益的大型水電工程。攔河壩型為礫石土心墻堆石壩，主要包含礫石土心墻、反濾層、過渡層、堆石區(qū)。在心墻上游和下游分別設(shè)置雙層反濾層，反濾層與上、下游壩殼堆石之間設(shè)有過渡層。壩頂高程856.00m，最大壩高186.00m。水庫正常蓄水位850.00m，校核洪水位853.78m。壩基防滲采用墻幕結(jié)合的形式，設(shè)主、副兩道混凝土防滲墻，副墻位于主墻上游側(cè)，防滲墻底設(shè)防滲帷幕，防滲墻與心墻及基巖防滲帷幕共同構(gòu)成主防滲平面。其典型河床0+240斷面如圖2所示。

圖2 PG水電站河床0+240斷面示意圖

2.2 有限元模型

根據(jù)壩體結(jié)構(gòu)輪廓尺寸、材料分區(qū)和壩址地形地質(zhì)條件等資料，本文選取大壩0+240斷面，建立二維有限元模型。根據(jù)該斷面滲壓測點位置在有限元模型中布置相應(yīng)結(jié)點，以便于提取滲壓計算結(jié)果。有限元模型的模擬范圍為：自建基面向下延伸2倍壩高，上下游方向自坡腳分別延伸2倍壩高。整個有限元模型共劃分為11555個單元。壩體分為8個區(qū)域，其中心墻以歷史最低水位為分界，分為心墻上、下區(qū)。有限元分析模型如圖3所示。壩體和壩基各類材料的設(shè)計滲透系數(shù)見表1。

圖3 PG礫石土心墻堆石壩0+240斷面各分區(qū)有限元網(wǎng)格

表1 壩體及壩基各區(qū)材料滲透系數(shù) 單位：cm/s

2.3 心墻滲流參數(shù)反演

通過敏感性分析，本文選取5個參數(shù)進行反分析，分別為心墻上、下區(qū)和弱風化基巖的滲透系數(shù)以及兩個非飽和土體參數(shù)hb和λ。選用2014年的上游庫水位輸入有限元進行非穩(wěn)定非飽和滲流場的連續(xù)計算，分別得到每組樣本對應(yīng)的各特征點的滲壓值。

將上半年作為擬合段，均勻選取時間控制點進行響應(yīng)面擬合。下半年作為預(yù)測段，用于驗證反分析結(jié)果的適用性。擬合段和預(yù)測段各有6個時間控制點。采用不含交叉項的3次多項式，對每個測點滲壓在各個時間節(jié)點進行響應(yīng)面方程擬合，其復(fù)相關(guān)系數(shù)均大于0.98，表明所構(gòu)建的響應(yīng)面代理模型對該工程的滲壓效應(yīng)量擬合效果良好，可用于之后的滲流參數(shù)反演。選取滲壓計特征點位置如圖4所示，上游水位歷時曲線及時間控制點如圖5所示。測點擬合選用P83進行展示，擬合情況如圖6所示。

圖4 PG心墻堆石壩心墻特征點示意圖

圖5 2014年P(guān)G堆石壩庫水位歷時曲線

圖6 P83滲壓響應(yīng)面模型擬合

根據(jù)滲壓測點響應(yīng)面代理模型及對應(yīng)的原觀監(jiān)測滲壓數(shù)據(jù)，采用歸一化的方法，通過每個測點不同時間控制節(jié)點對應(yīng)的實測滲壓值和響應(yīng)面方程模擬滲壓值的相對誤差平方和構(gòu)建目標函數(shù)。采用多目標優(yōu)化對目標函數(shù)求解，得到包含2000組解的Pareto解集。通過Pareto解集與反演參數(shù)設(shè)計值隸屬度函數(shù)尋優(yōu)獲得PG堆石壩心墻非飽和滲流參數(shù)反分析結(jié)果，見表2。

表2 參數(shù)反分析結(jié)果

將解代回有限元模型進行數(shù)值模擬計算，擬合效果選取心墻上區(qū)測點P83和心墻下區(qū)測點P20進行展示，如圖7—8所示，與實測值對比的誤差分析見表3。結(jié)果表明，代回有限元模型后可以較好地擬合實測滲壓的變化趨勢，擬合精度高，各測點的平均絕對誤差均值為10.1kPa，平均相對誤差均值為1.71%。

圖7 P83滲壓反演與實測值對比

表3 參數(shù)反分析有限元結(jié)果誤差分析

對比非飽和非穩(wěn)定滲流與飽和非穩(wěn)定滲流，非飽和特性對心墻下區(qū)滲壓影響較小，兩種滲流的心墻下區(qū)測點的平均絕對誤差分別為1.03、1.89m，但是對心墻非飽和區(qū)滲壓影響較大，在水庫水位下降過程中，兩種滲流的平均絕對誤差分別為2.87、9.77m，飽和非穩(wěn)定滲流的最大絕對誤差達到33.43m。因此非飽和滲流有限元模型不僅能反映飽和區(qū)滲流變化，也能反映非飽和區(qū)滲流變化，使心墻滲流的模擬更符合實際情況。

圖8 P20滲壓反演與實測值對比

2.4 非飽和滲流特性分析

根據(jù)反分析獲得的參數(shù)進行飽和滲流與非飽和滲流的滲流場孔隙水壓對比，選取2個典型時刻繪制對應(yīng)時刻的孔壓云圖如圖9所示。

圖9 典型時刻的飽和非穩(wěn)定滲流與非飽和非穩(wěn)定滲流孔壓

不考慮非飽和特性的非穩(wěn)定滲流，浸潤面變化基本同步于水庫水位變化，在心墻下游面迅速折減。考慮非飽和特性后，上游水位下降時，孔隙水壓力不能及時消散，心墻浸潤線變化滯后于水位的降低，出現(xiàn)心墻浸潤線高于水庫水位。水位上升時，存在由非飽和狀態(tài)至飽和狀態(tài)的過渡，心墻浸潤線變化滯后于水位的上升，且心墻內(nèi)浸潤線升至與水庫水位相同穩(wěn)定浸潤線需要經(jīng)過一段時間。

出現(xiàn)上述現(xiàn)象是由于心墻上區(qū)受到土體基質(zhì)吸力影響?；|(zhì)吸力隨土體飽和度降低而增加，心墻滲透系數(shù)隨吸力的增加，呈非線性的降低，上游水位下降時心墻內(nèi)水分由于基質(zhì)吸力來不及消散而滯留在心墻內(nèi)，堆石區(qū)不考慮非飽和特性且滲透系數(shù)較大，當堆石區(qū)浸潤線與水位保持同步時，心墻內(nèi)仍然維持相對較高的浸潤線，呈現(xiàn)“上凸”的拋物線形狀。這種情況容易引起心墻材料容重改變，帶來的非穩(wěn)定滲流力造成土體顆粒之間有效應(yīng)力減小，進一步導致土的抗剪強度降低，最終危機壩體穩(wěn)定性。水位降落速度越快，心墻浸潤線相對位置越高，其非穩(wěn)定滲流力越大，對穩(wěn)定性影響也越大。

計算典型時刻堆石壩心墻和防滲墻的滲透坡降并繪制云圖，如圖10所示。結(jié)果表明心墻孔壓和滲透坡降分布均勻，滲透坡降較大值出現(xiàn)在心墻中下部和防滲墻中上部，心墻出逸點附近滲透坡降為2～3，出現(xiàn)的心墻最大滲透坡降為2.72，小于心墻允許滲透坡降3。主、副防滲墻滲透坡降基本相近，防滲墻最大坡降出現(xiàn)在主防滲墻上部，出現(xiàn)的防滲墻最大滲透坡降為88.55，小于防滲墻允許滲透坡降120。

圖10 非飽和非穩(wěn)定滲流坡降云圖

3 結(jié)語

本文依托PG礫石土心墻堆石壩工程，構(gòu)建非飽和非穩(wěn)定滲流模型?；诙咽瘔螡B流原觀數(shù)據(jù)和響應(yīng)面代理模型，采用NSGA-Ⅱ多目標遺傳算法，進行滲透系數(shù)和非飽和特性參數(shù)的反演，其結(jié)果都在合理范圍內(nèi)。反演結(jié)果表明，反演得到的滲流參數(shù)對于心墻上區(qū)滲壓變化的擬合情況較為準確，能整體反映大壩心墻滲壓場的變化情況?？紤]心墻材料的非飽和特性，引入土-水特征曲線及滲透系數(shù)函數(shù)模型，更正了非飽和區(qū)滲透系數(shù)恒定不變從而導致非穩(wěn)定滲流求解的不準確性，使心墻非穩(wěn)定滲流的模擬更加符合實際情況。