隧道围岩渗透系数反馈分析

蒋 进 ，关劲松

（1.重庆交通大学，重庆 400042；2.辽源市东辽县公路管理段，吉林 辽源 136200）

0 前言

岩土渗流有关参数的正确取用是渗流计算分析中的重要环节。数值计算或模拟试验的精度再高，若选用参数不当也将导致计算分析结果的不可靠，而使确定的渗流控制方案不安全或过于浪费。目前计算机的精度已高达小数点后若干位，而参数选用中的误差却会高达数量级，这种计算中的精度失调问题实在应该引起注意。然而由于介质的复杂性，无论是实验室试验还是现场试验，不仅需要花费大量的人力、物力和时间，而且往往得不到满意的结果，造成大量的安全隐患。随着数值解法和计算机的发展，使利用数学模型反求参数成为可能。

本文将针对这一事实专门阐述如何利用大型通用有限元计算软件ANSYS的优化设计功能反馈分析渗透系数，使人们能够根据已有经验或研究成果对地质勘探资料和建筑材料性能做出比较可靠的估计，或者对取样利用室内试验设备测定及野外现场试验测定取得的参数的准确性做出评价和修正以得到更可靠的有关参数。

1 渗流场与温度场相似对比

渗流场中存在水头梯度,水缓慢地从高水头的地方流向低水头的地方；同样温度场中的热能在温度梯度作用下从温度高的地方向温度低的地方传导。 由此看来温度场和渗流场的传播过程具有较好的相似性，同时渗流与热传导的微分方程形式、控制变量，以及相关参数也具有很好的相似性，因此利用温度场来模拟渗流场在理论上是具有可行性的。

1.1 控制方程的相似性

对于渗流，假定土体和水体是不可压缩且介质均匀各项同性，则应用能力守恒定律可推导其控制方程为：

当温度场中无热源时，qgen=0,此时两式各项相似，只是其物理意义有改变，T表示温度而H表示水头，k表示热传导系数和渗透系数。

1.2 模型的相似性

温度场模型和渗流场模型相比拟，还应满足几何相似和边界条件一致。

（1）几何相似。温度场模型的外部边界应和所研究渗流区域的外部边界在几何上相似。当渗流区域为均质岩层时，则模型也应是均质的；当渗流区域是非均质岩层时,则要求模型中不同导热介质的分界线应与非均质岩层的分界线亦保持相似，此时应注意渗流的微分控制方程有所改变。

（2）边界条件一致。即温度场模型的绝热边界与渗流区域的隔水边界相对应。导热边界和透水边界相对应，导热边界上的温度则和透水边界上的水头相对应。

1.3 ANSYS单元理论相似

对于二维四节点矩形有限元单元，同1.2节中的分析，可以推导出其传导矩阵为：

其中 C1=kx，C2=ky，C3=0，k 表示渗透系数。

对比渗流场和温度场不难发现，其二维四节点矩形有限元单元的传导矩阵的推到过程和结果都相同，这就更进一步表明，利用ansys温度场，模拟渗流场是可行的。

2 渗透系数反演的实现方法

2.1 反演在ANSYS中的实现

2.1.1 何谓反演

所谓的反演就是预先根据经验估计或现有的勘测结果确定各岩层的渗透系数，然后应用ANSYS求解渗透场的水头分布，最后将求得的水头值与相应的观测值进行比较，并根据其差值和地质钻孔资料来修改各岩层的轮廓和渗透系数，如此迭代，反复求解直至满足精度要求为止。需特别说明的是，这种反求参数的方法并不能代替地质勘探工作，而仅仅只是水文地质工作的一个补充延续而已。其目的是为了减少人力物力的投入并且得到更准确的参数。

最后求得的自由面上的水头值与相应的观测值之间应满足如下关系[2]：

α≤k1≤β

式中，函数E称为目标函数，ω为权因子，一般精度要求越高，相应的权因子值越大，H和H′，分别为节点水头的求解值和观测值，α、β分别为k的上限和下限。根据能量守恒定律可知，当各岩层的渗透系数满足一定的比例条件时，线性方程组 [K][H]=[f]非唯一解，为消除这个不定解的因素，预先精确地给定某一岩层的渗透系数或由观测资料给出渗流区域的渗流量。

2.1.2 迭代计算渗流自由面

渗流自由事先是不确定的，确定其位置是渗流计算的主要内容。渗流自由面上的水头压力等于大气压力，该面上任意一点水头h等于该点的位置高程。为保证式（3）的解存在唯一解，在渗流自由面上应满足条件：

对非稳定渗流自由面除应满足式（4）外，还应满足第二类边界条件的流量补给关系，即按下式计算渗流自由面降下时由渗流自由面流入介质内的单宽流量q：

式中，μ为渗流自由面变动范围内的介质有效孔隙率或给水度；θ为渗流自由面外法向与铅锤线的夹角。

用有限单元法计算一个物理问题，它要求计算区域的边界必须是完全确定的。然而隧道围岩渗流是有自由面的无压渗流，渗流自由面位置恰恰又是不确定的，渗流自由面位置是工程设计十分关心的问题，故需通过计算求出。渗流自由面需满足式（4）的条件，这样一来，渗流自由面位置需在计算中通过迭代来确定。迭代的计算步聚是[2]：

（1）首先根据渗流概念和经验假定一条渗流自由面，以确定有限单元法的计算区域。

（2）将假定的渗流自由面作为第二类边界，对于稳定的渗流由式 [K][H]=[f]计算出渗流自由面结点的水头值h。

（3）比较假定渗流自由面结点的计算水头值和其位置高程z视其是否满足式（4）的条件。若不满足，则用计算水头值去改变结点的z坐标，形成新的假定渗流自由面，同时还确定了新的有限单元法计算的区域。反复上述计算步骤直至渗流自由面上的结点全部满足给定的允许精度（|h-z|≤ε）。

在迭代过程中为了不改变网格的划分和避免单元在压缩或伸展的过程中出现畸形，可以采用丢弃单元（节点）的方法。即是当渗流自由面结点T的计算水头值小于或等于其下结点的z坐标时，将T结点丢弃，改S结点为渗流自由面结点，形成新的有限单元法汁算区域，而与T结点有关的单元则不再参加计算。

2.2 ANSYS优化设计

ANSYS优化设计有两种方法：第一种是通用的函数逼近优化方法[3]。这种方法的本质是采用最小二乘逼近，求取一个函数面来拟合解空间，然后再对函数面求极值。这种方法是一种普适的优化方法，不容易陷入局部极值点，但优化精度不是很高。另一种方法是对第一种方法的改进方法，叫做梯度寻优。如果说第一种方法是C0阶、大范围普适的粗优方法，那么第二种就是C1阶、局部寻优的精优方法。一般来说，一个比较复杂的问题需要同时采用两种优化方法，先用第一种初步求得优解的基本范围，然后再采用第二种方法对优解的范围进行更精确的确定。但是第二种方法不仅耗时间，而且还容易陷入局部最小点，因此通常的问题利用C0阶函数逼近优化就足够了。ANSYS中优化设计模块有三大变量：（1）设计变量DV:即自变量。（2）状态变量SV：用来体现优化的边界条件，是设计变量的函数。（3）目标变量OBJ：最终的优化目的。它必须是设计变量的函数，而且只能求其最小值。概括说来ANSYS中实现优化设计的过程就是：在状态变量SV的边界条件下，在设计变量的取值范围内不断调整DV的取值组合，以使得目标函数OBJ逼近最小值。

3 渗透系数反演的工程实例

3.1 工程概况

为了验证以上数学模型在ANSYS数值模拟环境中分析的准确性，以及在实际工程中实用性，现以八台山隧道为实例，在ANSYS中进行结构分析。隧道处于岩体强风化带、碎屑岩类的构造裂隙发育区、地表岩溶分布区为地下水的补给区，主要接受大气降水的补给。

3.1.1 原始资料

根据原始资料，选取某断面四组实测水头值，其值如表1所列。

表1 实测水头值统计表

3.1.2 模型参数

模型计算范围选取埋深60m，隧道底部以下32m，轴线左右各60m，故整个模型为120m×100 m。左边界水头为64.8 m，右边界水头值为9.4 m（该值为相对值，相对于模型的原点）,采用ansys plane 55单元。该模型由于是反求隧道围岩的渗透系数，所以没有模拟衬砌和加固圈，并且由于实际计算中对超过三维的渗透系数进行反分析的理论还不成熟，而且工程当中又有一些渗透系数比较稳定，比较容易确定。例如衬砌的渗透系数，在实验室中就可以比较准确地测定出来。显然这类系数就没有反演的必要性。另外为保证反分析解的唯一性，并且由于干燥状态下围岩的渗透系数对渗流量的影响很小，跟设计时的出入不大，所以选取其设计值，不参与优化设计，而作为优化结果的控制条件。具体的参数取值如表2所列，其中干燥状态下的围岩渗透系数为一个确定的值，而饱和状态下的围岩渗透系数为一个设计优化区域。

表2 渗透系数一览表

3.2 优化设计

经过10次优化设计参数搜索最终得到10组优化设计组合，其中带“*”号的为系统推荐的最优组合，最优设计又包含了33次迭代计算，并得最优组合*SET10*中涌水量为7.402 3m3/d,围岩饱和渗透系数为0.451 79m/d。将ANSYS反演计算的结果和实际测量的结果对比可以得到以下数据和结论（见表3）。

表3 ANSYS反演计算的结果和实测的结果对比表（单位：m）

由表3可以看出，隧道围岩饱和渗透系数采用优化设计最优解时，计算水头和实际勘测水头相对残差均值为-0.405，最大残差为-0.873，而标准残差为0.508，这说明计算水头值与实测水头值拟合很好，能够满足精度要求。另外状态变量（涌水量）的计算值为7.406，而实测量为7.42，也表明了所求渗透系数的合理性。

4 结论及问题

通过ANSYS优化设计模块，成功地反演出隧道围岩的渗透系数，可以方便地检验实际勘测的渗透系数的准确性，在一定程度上减少了野外实验在人力、物力方面的消耗。对同等地质条件下的工程提供了一个参考思路。虽然没有对隧道加固圈进行分析，不过同样提供解决加固圈厚度的有限元计算方法—优化设计，以加固圈厚度为设计变量，分析加固圈厚度与涌水量的关系，利用ANSYS温度场模块可以较好地实现对均匀介质渗流场的模拟。不过实际工程问题中，特别是岩体中是没有均匀介质的，只有对模型做一定的简化，才可以应用。该例正是因为隧道断面处风化严重、岩层破碎厉害，并且没有呈有规律方向的裂缝，所以才简化为均匀介质。渗透系数的反演在一定程度上解决了在野外测定渗透系数的困难。本文在一定程度上说明了，可以利用ANSYS的温度场求解渗流场的问题，不过由于渗流场介质的复杂性，在应用时还需考虑诸多因素对刚度矩阵的影响。ANSYS优化设计模块的巧妙应用对实际工程来说，无论在安全上还是在经济上都是大有益处的。

[1]王崧,刘丽娟，董春敏，等译.有限元分析—ANSYS理论与应用（第三版）[M].北京：电子工业出版社，2011：296～316.

[2]毛昶熙.渗流计算分析与控制 (第二版)[M].北京：中国水力水电出版社，2003：306～316.