基于响应面法的岩溶区堤坝注浆帷幕设计优化研究

谢梦珊，曹 磊，郭晓静，刘 旭，张 聪，蒋买勇，李金友

(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的溶蚀率减小；随着帷幕高度的增加，可溶蚀Ca(OH)2含量增大，此时帷幕体内Ca(OH)2溶蚀率也随之增大，当帷幕增加到某一高度，帷幕内部Ca(OH)2的增加量远大于溶蚀的Ca(OH)2增加量，致使Ca(OH)2溶蚀率减小；随着帷幕宽度的增加，帷幕内部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函数建立了岩溶区注浆帷幕设计参数优化模型，并将优化模型应用于某岩溶区险情堤坝的帷幕设计优化。提出的优化设计方法能弥补经验性设计原则的不足，为岩溶区堤坝注浆帷幕设计参数优化提供一种新思路，且对于提高帷幕施工质量和控制帷幕工程造价具有显著的工程价值与意义。