熊彬等：有限单元法分析堤坝管涌拟流场特征与一种快速检测渗漏的测量新方式
1 引言
截至2020年底，全国已建成5级及以上江河堤防32.8万km，累计达标堤防24.0万km，堤防达标率为73.0%，27.0%的堤防存在渗漏隐患。堤防工程中常见的隐患包括天然地质缺陷、施工中的质量缺陷、生物破坏造成的洞穴、裂缝及修补时的人为薄弱环节等，每年汛期，这些隐患容易诱发管涌、漫顶、散浸、滑坡、崩岸、坍塌等险情。为了及时查明隐患，预先排除险情是堤防工程中迫切要解决的问题。目前，堤坝渗漏检测的主要技术还是依据各种物探手段，如面波法、瞬变电磁法、高密度电法等。流场法是中南大学何继善院士提出的一种堤坝管涌渗漏探测新方法，其原理是在背水面的堤垸和迎水面的水中同时发送一种人工信号——特殊波电流场——去拟合并强化异常水流场的分布，通过分析电流场的分布来查明水流场的动向。由于流场法侧重于探测管涌渗漏的入水口部位，因此它的探测结果可以直接用于汛期堤防管涌渗漏的封堵，为在汛期从根本上消除管涌渗漏的危害提供经济、快速的技术支持。戴前伟等利用伪随机流场法与双频激电法对深圳某水库进行综合勘探，有效排除了金属管对解释带来的困难。白广明等建立拟流场数学模型，并用有限单元法计算得到拟流场电位、电场强度及电流密度分布图，为数据解释提供了指导建议。房瑞等通过研究发现，采用接近拐角频率的1024Hz伪随机信号能够获得了更稳定的数据，且异常响应明显，抗干扰能力强。何继善院士采用广域电磁法和流场法组合，准确获得了矿区地下采空区和水害分布，为矿区安全生产提供了可靠的技术保障。近年来，学术界对于堤防渗漏隐患的研究已经非常丰富，但结合具体的研究成果来看，这些研究内容更多地专注于对后期的数据处理方面。本文在研究过程中，利用有限元方法对不同模型进行分析，总结出汛期快速检测渗漏入水口的测量方式，即在水库底部设置测线测量电流密度z分量，通过分析电流密度z分量的异常响应推断渗漏入水口处位置，并应用在深圳某水库堤防的渗漏通道入水口检测的实例中，圈定了疑似渗漏入水口区域，经过工程验证，证实其为渗漏入水口。式中，u为电位，单位为V；σ为地下介质电导率，单位为S/m；I为供电电流强度，单位是A；n为界面外法线方向；Γ∞当无穷远边界；r为点电流源A到Γ∞的距离，单位为m；Ωs为计算区域；Γs为地面。将边值问题转化为变分问题：
将计算区域剖分成若干个四面体单元，如图1所示为节点局部编号示意图，将泛函问题离散化。图1 剖分的四面体单元节点编号示意图
其中，Ne为单元个数；Ωe为单元区域；Γe为无穷远边界单元区域。
单元内节点以外的电位使用线性插值得到：
式中，Ni(i=1，2，3，4)为插值形函数；ui(i=1，2，3，4)为节点单元电位。将插值函数代入离散公式，化简后得：
其中，ke为单元系数矩阵，组装单元系数矩阵得到总体系数矩阵方程如下：
本文采用压缩行方式(CRS，Compressed Row Storage)存储稀疏矩阵以节约空间，使用不完全柯列斯基分解(Cholesky)共轭梯度法求解式(6)大型线性方程组获得所有网格节点上的电位u。由下式计算电流密度响应：
对正常堤坝电流密度背景场特点和规律进行分析，有助于后续对电流密度各分量异常的分析。图2为正常堤防模型截面图。
图2 正常堤防模型截面
假设堤防模型由堤身、基岩和水体组成，堤顶高12m，堤底部宽52m，水体深度为10m，各项参数具体如表1所示。
表1 堤坝模型各项参数
利用有限单元法计算模型电流密度各分量响应曲线如图3所示。其中，L01，L02，L03，L04为过堤防迎水坡坡脚，分别位于水底、距水底3m、距水底6m、水面的测线。在正常堤防模型中，电流密度x分量存在一个双极性异常，这个双极性异常可以由电法勘探相关理论解释，也从侧面证明算例结果是可信的。电流密度y分量存在一个正的极大值，且在整个侧线上电流密度方向不发生改变。而电流密度z分量在点源投影点上有一个负的极大值异常，这是因为受到堤坝的影响，电阻率最短路径位于水中，因此在点源投影点上出现负的极大值异常。
图3 无渗漏模型电流密度各分量曲线
通过分析电流密度各分量在水中不同深度的曲线可知，电流密度x分量对深度变化没有明显响应。电流密度y分量从水底至水面幅值增大，而电流密度z分量则完全相反。因此x分量不受水深影响，y分量在水面幅值大，而z分量在水底幅值大。如图4截面图所示，在正常堤坝模型上设计了一条沿y轴从坝体迎水面底部到坝体背水面底部渗漏通道，渗漏入水口坐标为(100，200，0)，出水口坐标为(100，148，0)，渗漏通道电阻率如表1所示。
图4 渗漏通道模型截面
利用有限单元法分析上述模型，得到模型的拟流场法数值模拟结果。图5为有渗漏模型的电流密度各分量响应。分析图5可知，电流密度各分量在渗漏入水口处均有异常响应，其中x分量在渗漏入水口处出现一个双极性异常，y分量和z分量在入水口处出现一个正的极大值异常。但x分量异常幅度较小，y分量异常容易和正常堤坝模型电流密度响应混淆。并且各分量异常响应离水底越近，其异常响应越明显。因此，在实际测量中，一般在水底布设测线，测量电流密度z分量异常。
图5 有渗漏模型电流密度各分量曲线
考虑到渗漏通道倾斜对堤坝电流密度响应的影响，设计了一个渗漏通道倾斜模型，如图6所示为渗漏通道倾斜模型xy平面图。渗漏通道沿x轴方向倾斜，从坝体迎水面底部到坝体背水面底部，渗漏入水口坐标为(90，200，0)，出水口坐标为(100，148，0)。
图6 倾斜渗漏通道堤坝模型
图7(a)为图4的模型y=200m平面电流密度矢量图，图7(b)为图6的模型y=200m平面电流密度矢量图。可以明显发现，在虚框内也就是渗漏入水口附近，电流密度矢量出现方向改变的异常，越靠近渗漏入水口，越接近水底，其异常越明显，异常出现的位置与渗漏入水口位置保持一致。这进一步说明了实际测量应在水底布设测线，其异常响应更加明显。
图7 电流密度矢量图
以上3个模型考虑的是渗漏通道从堤坝迎水面到堤坝背水面的情况，但实际上，有些渗漏通道是从水体的底部基岩到堤坝后面的基岩。为考虑基岩渗漏模型对电流密度z分量的影响，构建基岩渗漏通道模型。如图8所示为基岩渗漏通道模型截面图，渗漏通道为基岩渗漏，从坝前水域底部到坝后底部，渗漏入水口坐标为(100，140，0)，渗漏出水口坐标为(100，220，0)。
图8 基岩渗漏通道模型
利用有限单元法分析上述模型，得到模型的电流密度z分量响应，如图9所示。其中，D01、D02、D03、D04是分别位于Y=210m、Y=220m、Y=225m、Y=230m的水底测线。可以明显发现，基岩渗漏模型的电流密度z分量曲线与堤坝渗漏模型的电流密度z分量曲线不一致，堤坝渗漏模型z分量在远离渗漏入水口处的方向为负向，在渗漏入水口处出现正的极大值点，而基岩渗漏模型整体电流密度z分量均为正向，在渗漏入水口处出现正的极大值点。其相同点均在渗漏通道入水口处出现极大值点，且越靠近渗漏通道，其异常越明显。因此，在实际观测中，可以通过在坝体前面水域底下沿坝体走向方向布设不同测线，测量其电流密度z分量，当z分量出现极大值点，可以确为定渗漏入水口处。深圳某水库发现二期泄洪隧洞出口洞脸处山体排水管中有水流出，水体外观呈透明清澈。由于周边排水系统不完善，长期被水体侵扰，不利于管理。长久的渗漏，很可能会导致更为严重的地质灾害。为了解决这个难题，决定对水库堤坝进行渗漏探测。测量区域如图10所示，将堤坝管涌发射机放置在4号坝与1号坝连接处附近的4号哨所内，一端连接铜电极A，将铜电极A放置在渗透出水区域，其中一根电极直接放置在有明显水流的出水口处，发射机另一端连接铝电极B，铝电极B放在一号坝北面尽头离岸20m的水中。发射电压65V，发射电流200mA，将堤坝管涌接收机和RTK定位仪器放置在船上，由实验人员在船上拖动探头，根据水深调整探头入水的深度，按照提前布置的测线进行扫描测量。其电流密度z分量平面图如图11所示。可以明显发现，在启闭塔附近很多区域出现了较高的异常值，在柱子的四个角疑似为起闭塔未封严实造成，而在起闭塔右边局部的高异常值疑似为渗漏区域。本文采用基于非结构化网格有限单元法实现了拟流场法三维有限元正演，并对多种堤坝管涌模型进行正演计算，通过对模型数值结果进行分析，发现当堤坝出现渗漏隐患时，其电流密度各分量都会出现一些异常。电流密度x分量会在渗漏处出现一个双极性异常。而电流密度y分量在渗漏处出现一个异常极值点。但电流密度x分量和电流密度y分量异常幅度较小，而电流密度z分量的异常幅度很大，出现z方向正向的极值点，且堤坝渗漏时，远离渗漏通道入水口处的电流密度z分量为负值，而基岩渗漏时，远离渗漏通道入水口处的电流密度z分量为正值。因此，在实际测量中，应多加考虑电流密度z分量的异常。通过分析不同测线的电流密度，发现越靠近渗漏入水口处，其异常幅值越高。由此可以考虑在水域底部布设侧线，能更好地反应渗漏入水口位置。通过实例分析，进一步证实了电流密度z分量在实际测量中能很好地反映出渗漏点所在位置。
