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

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

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

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

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

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—重力加速度。

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

(2)

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

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

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

(3)

(4)

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

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

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

(5)

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

(6)

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

1.4 基于NSGA-Ⅱ多目標(biāo)遺傳算法

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

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

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

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

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

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

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

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

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

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

(8)變異準(zhǔn)則：為避免局部最優(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ù)該斷面滲壓測點(diǎn)位置在有限元模型中布置相應(yīng)結(jié)點(diǎn)，以便于提取滲壓計(jì)算結(jié)果。有限元模型的模擬范圍為：自建基面向下延伸2倍壩高，上下游方向自坡腳分別延伸2倍壩高。整個(gè)有限元模型共劃分為11555個(gè)單元。壩體分為8個(gè)區(qū)域，其中心墻以歷史最低水位為分界，分為心墻上、下區(qū)。有限元分析模型如圖3所示。壩體和壩基各類材料的設(shè)計(jì)滲透系數(shù)見表1。

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

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

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

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

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

圖4 PG心墻堆石壩心墻特征點(diǎn)示意圖

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

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

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

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

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

圖7 P83滲壓反演與實(shí)測值對(duì)比

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

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

圖8 P20滲壓反演與實(shí)測值對(duì)比

2.4 非飽和滲流特性分析

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

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

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

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

計(jì)算典型時(shí)刻堆石壩心墻和防滲墻的滲透坡降并繪制云圖，如圖10所示。結(jié)果表明心墻孔壓和滲透坡降分布均勻，滲透坡降較大值出現(xiàn)在心墻中下部和防滲墻中上部，心墻出逸點(diǎ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-Ⅱ多目標(biāo)遺傳算法，進(jìn)行滲透系數(shù)和非飽和特性參數(shù)的反演，其結(jié)果都在合理范圍內(nèi)。反演結(jié)果表明，反演得到的滲流參數(shù)對(duì)于心墻上區(qū)滲壓變化的擬合情況較為準(zhǔn)確，能整體反映大壩心墻滲壓場的變化情況?？紤]心墻材料的非飽和特性，引入土-水特征曲線及滲透系數(shù)函數(shù)模型，更正了非飽和區(qū)滲透系數(shù)恒定不變從而導(dǎo)致非穩(wěn)定滲流求解的不準(zhǔn)確性，使心墻非穩(wěn)定滲流的模擬更加符合實(shí)際情況。