基于響應面法的巖溶區堤壩注漿帷幕設計優化研究

謝夢珊，曹 磊，郭曉靜，劉 旭，張 聰，蔣買勇，李金友

(1.中南林業科技大學土木工程學院，湖南 長沙 410004；2.湖南水利水電職業技術學院，湖南 長沙 410131)

20世紀50—70年代，中國修建了大量堤壩工程，受當時勘察設計水平、機械化施工水平以及經濟條件等限制，部分巖溶區堤壩存在諸多先天不足，而后隨著服役年限的增長，多數堤壩出現了不同程度的病險，已嚴重影響了其安全運營[1-2]。為此，近十年來全國多次開展大規模病險堤壩除險加固工程，也隨之形成了大量的堤壩除險加固技術[3-4]。其中，注漿加固技術因具有施工便捷、安全環保、經濟效益高等優點，已逐步成為堤壩除險加固工程的首選[5-6]。然而，大量工程實踐發現部分巖溶區堤壩在注漿除險加固完成會再次出現滲漏險情，分析原因為巖溶區地質條件十分復雜，現有地質勘察技術無法實現對地層的精準探測，而不合理的注漿帷幕設計又進一步導致了該類工程問題的發生。因此，開展巖溶區堤壩注漿帷幕設計優化研究具有十分重要的現實意義和實用價值。

目前，國內外學者對于巖溶區注漿帷幕設計優化方面的研究尚少，研究主要集中在非巖溶區注漿帷幕設計優化方面，且研究手段以理論推導和數值模擬為主。在注漿帷幕設計理論研究方面，王恩志等[7]利用三維滲流計算，對地層滲透性、帷幕的防滲性能及帷幕深度和長度進行了敏感性分析；崔文娟等[8]分別采用遺傳算法和復合形法，從灌漿帷幕整體效應最大出發尋求最優帷幕參數的組合；黨發寧等[9]提出了基于平衡防滲原理的防滲帷幕優化設計方法；Chai等[10]從最小滲透壓力的角度提出了帷幕灌漿的最佳厚度理論。盡管這些理論對于研究注漿帷幕參數以及注漿帷幕設計優化具有重要的參考價值，但其推導過程大多數均以現有的注漿帷幕參數經驗公式為基礎，且假定的地層條件較單一，缺乏在復雜地質條件下推廣應用的普適性。為解決理論推導的諸多不足，多數學者嘗試采用數值模擬方法來研究上述問題，如杜風雷等[11]采用GMS軟件，進行了防滲隔污帷幕的設計參數敏感性分析；基于Modflow數值模擬軟件，陳富東等[12]針對影響止水帷幕隔水效果的因素開展了正交試驗，Yang等[13]通過詳細的現場調查數據和數值模擬，提出了一套防滲帷幕優化設計程序。上述研究成果對于帷幕設計參數優化有一定的促進作用，但大多數研究僅從單因素設計角度提出帷幕優化設計，未考慮到多種設計因素交互影響下帷幕服役性能的變化。此外，為有效模擬巖溶水對注漿帷幕的滲水壓力與水化學溶蝕作用，數值仿真模型還需考慮到巖溶水環境下注漿帷幕微觀結構的改變。

鑒于此，本文首先采用多場耦合軟件構建巖溶區堤壩帷幕服役模型，提出巖溶區堤壩帷幕服役性能評價標準，并應用于某巖溶區堤壩帷幕。其次，通過Box-Behnken響應面試驗設計，建立巖溶區帷幕滲透系數、帷幕高度、帷幕寬度與帷幕體內Ca(OH)2溶蝕率之間的回歸模型。最后，考慮到帷幕的服役性能要求與工程的經濟性，基于Matlab中的fmincon函數建立巖溶區注漿帷幕設計參數優化模型，求解某巖溶區險情堤壩帷幕設計參數的最優組合。研究結果可彌補現有巖溶區注漿帷幕設計方法的不足，且對于合理選取帷幕設計參數、提高帷幕施工質量和控制帷幕工程造價具有顯著的工程價值與意義。

1 巖溶區堤壩帷幕服役性能評估

1.1 巖溶區堤壩帷幕服役模型

1.1.1滲流溶蝕耦合方程

a)滲流方程。假設壩基和帷幕體為均質多孔連續介質，地下水在多孔介質內流動服從達西定律，則巖溶水運動方程可見式(1)[14]：

(1)

式中ρ——巖溶水密度；ε——壩基孔隙率；μ——巖溶水動力黏度；p——壓力；D0——位置水頭；k——壩基滲透系數；Qm——源匯項。

b)鈣離溶解遷移方程。考慮現有注漿材料主要為水泥基漿材，故注漿帷幕中水泥結石體(主要成分為Ca(OH)2)在巖溶水環境下會發生溶蝕，溶蝕產生的鈣離子受擴散作用與滲流作用影響會從高濃度處向低濃度處擴散，并從高水壓處向低水壓處輸送。為此，假定帷幕溶質不與巖溶水中的離子發生化學反應，且將溶蝕產生的溶質通量作為源匯項并入溶質運移方程中。以Ca2+為目標離子，溶蝕過程中鈣離子的遷移方程[15]可見式(2)：

(2)

式中θ——孔隙度；i——大壩壩基不同部位；C——溶蝕Ca2+的濃度；DL、DT——縱向、橫向彌散系數；ux、uy——地下水水平、垂直方向實際流速；R——壩基中Ca2+的增長速率，其數值等于Ca(OH)2溶解速率。

根據文獻[16]可建立Ca(OH)2溶解速率與溶液中離子濃度的關系，見式(3)：

(3)

式中Rs——Ca(OH)2的溶解反應速率，其數值與壩基中Ca2+的增大速率相同；A——反應動力學系數；CCa2+、COH-——溶液中Ca2+、OH-濃度；n——動力學指數；Ksp——Ca(OH)2的溶度積常數。

c)地層物性參數演化方程。已有研究表明，水泥加固體在完全溶蝕后，其孔隙度增大引起的體積變化與硬化水泥漿體中Ca(OH)2的體積變化基本相同，C-S-H析鈣不會對加固體孔隙度產生較大影響[17]。因此，本文僅考慮Ca(OH)2溶蝕對帷幕孔隙度的影響，則溶蝕過程孔隙度的增量Δθ的計算見式(4)：

(4)

式中Vm——Ca(OH)2的摩爾體積；M——Ca(OH)2的溶蝕率；CCH——固相Ca(OH)2的濃度；S——溶蝕Ca2+的分布區域；H——帷幕高度；D——帷幕寬度；MCH——Ca(OH)2的摩爾質量；ρCH——固相Ca(OH)2的密度。

1.1.2模型構建

湖南某水庫至今已運行50多年，由于區內發育大量溶洞、溶槽及巖溶塌陷等，下部巖溶通道滲漏嚴重，多年來險情不斷，安全度汛形勢極為嚴峻。現欲采用帷幕注漿技術降低地層滲透性，阻斷地下水滲漏通道，從而保證水庫運行及周圍環境的安全。為優化注漿帷幕設計方案，采用多場耦合軟件Comsol Multiphysics構建注漿帷幕服役狀態模型?？紤]模型尺寸效應，研究區域取壩踵向左延伸80 m，壩腳向右延伸100 m，深度取5倍壩高，帷幕深度15 m、厚度5 m；堤壩上游水位10.97 m、下游水位0.20 m，設為水頭邊界；網格劃分采用自由三角形網格，對帷幕邊界區域進行加密處理，共生成單元4 109個網格，壩基水文地質參數見表1，模型示意見圖1。

表1 壩基水文地質參數

圖1 數值計算模型示意

選用軟件中Porous Media and Transmission of Rare Materials模塊(Darcy Law)與Chemical Transfer模塊(Transfer of Rare Matter in Porous Media)建立帷幕滲流與鈣溶解遷移多場耦合計算模型，模型計算參數見表2。

表2 模型計算參數

1.2 帷幕服役性能評價標準

根據文獻[18]可知，帷幕中結石體強度將在CaO累計溶出率大于25%時急劇下降，與此同時帷幕防滲能力也顯著降低。本文為考慮巖溶水對帷幕的滲流溶蝕作用，采用COMSOL內置后處理工具計算服役一定年限后帷幕溶蝕出的鈣離子總量，再與帷幕體內Ca(OH)2總量相比計算得到帷幕體內Ca(OH)2的溶蝕率，其Ca(OH)2的溶蝕率即為巖溶區堤壩帷幕服役性能優劣的評價指標。此外，本文以巖溶區帷幕體內Ca(OH)2溶蝕率超過25%作為帷幕失效判據，帷幕此時已不具備防滲功能。

1.3 服役性能評估

帷幕服役性能的劣化與壩基滲流壓力和其體內Ca(OH)2的溶蝕率有密切聯系，壩基揚壓力和溶蝕鈣離子濃度分布見圖2。

a)揚壓力分布(Pa)

由圖2可知，幕后揚壓力明顯小于幕前揚壓力，帷幕對水流具有明顯的阻滯效果；溶蝕出的鈣離子在壩基揚壓力的作用下發生對流擴散，引起帷幕上下游鈣離子濃度增大；隨著帷幕服役年限增加，帷幕附近鈣離子的分布范圍變大、濃度增大；計算得到10、30、100 a后帷幕體內Ca(OH)2溶蝕率分別為3.02%、8.70%、28.36%，說明鈣離子溶蝕量隨服役年限逐漸增大，且100 a時帷幕體內Ca(OH)2溶蝕率超過25%，即帷幕此時已失效，不具備防滲功能。

2 響應面法注漿帷幕優化設計

2.1 響應面試驗設計

大量工程實踐表明，帷幕滲透系數、帷幕深度以及帷幕寬度是影響帷幕服役性能的主要因素，也是帷幕設計過程中普遍關注的設計參數，故選取帷幕滲透系數、帷幕高度以及帷幕寬度為因素開展Box-Behnken響應面試驗設計。其中，各因素的水平取值參考相關文獻[19]、標準規范[20]以及工程實踐確立，見表3。

表3 因素及其水平

2.2 響應面結果分析

采用構建的巖溶區堤壩帷幕服役模型，對響應面樣本點設計的16組工況進行數值計算，不同工況下100 a后帷幕體內Ca(OH)2溶蝕率計算結果見表4。

表4 響應面樣本點設計及計算結果

2.2.1模型回歸分析

采用二次多項式回歸方程對表4中的試驗結果進行擬合分析，溶蝕率與帷幕設計參數的回歸模型見式(5)。

M=-317.55811+1.30305×108K+45.58989H+3.21466D-4.88889×106KH-5.79125×105KD+0.12533HD-1.55813H2-0.62568D2

(5)

式中M——帷幕體內Ca(OH)2的溶蝕率；K——帷幕滲透系數；H——帷幕高度；D——帷幕寬度。

為檢驗回歸模型的有效性，對響應面回歸模型開展方差分析與顯著性檢驗(表5)。由表5可知，回歸項K、H、D與平方項K2、H2均顯著(P<0.05)，表明滲透系數、帷幕寬度均對溶蝕率影響顯著，且對于M的影響程度由大到小依次為K、D、H。此外，在試驗設計范圍內模型的P<0.000 1，表明響應面模型為極顯著水平，且擬合精度好，近似模型有效。

表5 響應面回歸模型的方差分析與顯著性檢驗結果

2.2.2交互作用分析

根據響應面回歸模型得到各帷幕設計參數對帷幕體中Ca(OH)2溶蝕率交互作用的響應曲面，見圖3。

a)高度-初始滲透系數

由圖3a可知，選取不同的帷幕滲透系數進行分析，隨著帷幕高度增加，帷幕體內Ca(OH)2溶蝕率呈現出先增大再減小的趨勢；選取不同的帷幕寬度進行分析，隨著帷幕滲透系數的減小，帷幕體內Ca(OH)2溶蝕率大幅減小。由圖3b可知，選取不同的寬度進行分析，隨著帷幕滲透系數的減小，帷幕體內Ca(OH)2溶蝕率顯著減??；選取不同的帷幕滲透系數進行分析，隨著帷幕寬度的增加，帷幕體內Ca(OH)2溶蝕率呈現減小的趨勢。由圖3c可知，選取不同的帷幕高度進行分析，隨著帷幕寬度的增加，帷幕體內Ca(OH)2溶蝕率成緩慢增加的趨勢；選取不同的帷幕寬度進行分析，隨著帷幕高度的增加，帷幕體內Ca(OH)2溶蝕率，呈現出先增大再減小的趨勢。

分析其原因，隨著帷幕滲透系數的減小，巖溶水滲流穿過帷幕內部孔隙骨架的難度增大，進而減緩了巖溶水對帷幕的溶蝕作用，導致Ca(OH)2的溶蝕率減?。浑S著帷幕高度的增加，可溶蝕Ca(OH)2含量增大，此時帷幕體內Ca(OH)2溶蝕率也隨之增大，當帷幕增加到某一高度，帷幕內部Ca(OH)2的增加量遠大于溶蝕的Ca(OH)2增加量，致使Ca(OH)2溶蝕率減?。浑S著帷幕寬度的增加，帷幕內部Ca(OH)2含量增加，其增加量大于Ca(OH)2溶蝕量的增加量，故Ca(OH)2溶蝕率逐漸減小。

2.3 響應面參數優化

最優的注漿帷幕設計參數需使注漿帷幕同時滿足安全性和經濟性雙層要求。本文以前述依托工程為例，以滿足帷幕服役性能條件下最小帷幕注漿量為優化原則，采用Matlab中的fmincon函數建立巖溶區注漿帷幕設計參數優化模型。

2.3.1目標函數與約束函數

根據設計要求，工程帷幕滲透系數K需小于1.00×10-8m/s，服役年限需達百年工程要求，即保證百年服役期間注漿帷幕具備一定的服役性能。因此，以單位長度帷幕所需要的注漿量為目標函數、100 a時帷幕體內Ca(OH)2的溶蝕率為約束函數，進行帷幕寬度和帷幕高度的最優值求解。目標函數見式(6)，約束函數見式(7)。

S=HD

(6)

式中S——單位長度帷幕所需要的注漿量，L/m。

M=-317.55811+1.30305×108K+45.58989H+3.21466D-4.88889×106KH-5.79125×105HD+0.12533HD-1.55813H2-0.62568D2≤25

(7)

2.3.2優化參數求解

經過運算，本工程的注漿帷幕參數的優化組合為H=15.80 m、D=4.51 m。

3 結論

a)以某病險巖溶區堤壩為例，采用COMSOL Multiphysics多物理場耦合軟件建立了巖溶區注漿帷幕滲流-溶蝕耦合模型，并提出巖溶區堤壩帷幕服役性能評定標準，為評估帷幕服役性能優劣、預測帷幕服役壽命提供了一種有效方法，但模型方程中沒有具體考慮到Ca(OH)2除外的其他組分以及其他離子化學反應，如何構建更加合理、科學的模型值得進一步研究。

b)基于Box-Behnken響應面試驗設計建立了巖溶區帷幕滲透系數、帷幕高度、帷幕寬度與帷幕體內Ca(OH)2溶蝕率之間的回歸模型，回歸模型達到極顯著水平，即可適用于預測數值模擬計算結果；通過響應曲面分析了各因素對帷幕體內Ca(OH)2溶蝕率的影響規律，結果表明對帷幕服役性影響從小到大為帷幕高度、帷幕寬度、帷幕滲透系數。

c)以滿足帷幕服役性能條件下最小帷幕注漿量為優化原則，采用Matlab中的fmincon函數建立了巖溶區注漿帷幕設計參數優化模型，并將優化模型應用于某巖溶區險情堤壩的帷幕設計優化。提出的優化設計方法能彌補經驗性設計原則的不足，為巖溶區堤壩注漿帷幕設計參數優化提供一種新思路，且對于提高帷幕施工質量和控制帷幕工程造價具有顯著的工程價值與意義。