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

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

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

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

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

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

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

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

(1)

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

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

(2)

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

1.2 BC&BC非飽和滲流本構模型

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

(3)

(4)

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

1.3 滲流參數反演方法

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

(5)

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

(6)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1 工程概況

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

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

2.2 有限元模型

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

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

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

2.3 心墻滲流參數反演

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

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

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

圖5 2014年PG堆石壩庫水位歷時曲線

圖6 P83滲壓響應面模型擬合

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

表2 參數反分析結果

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

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

表3 參數反分析有限元結果誤差分析

對比非飽和非穩(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 非飽和滲流特性分析

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

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

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

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

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

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

3 結語

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