技术领域
[0001] 本
发明属于
水电工程中渗流场数值模拟领域，涉及基于三维精细地质模型的土石坝渗流场数值模拟方法。具体讲，涉及包括不良地质体和裂隙的精细地质模型的构建方法和基于精细地质模型的土石坝渗流场数值模拟方法，以及土石坝浸润面求解方法。
背景技术
[0002] 在水利工程中，水对工程设施的破坏是影响工程安全最主要的方面。其中
地下水对水利工程的作用不易观察，而且一经发现就难以挽回，经常会给工程带来重大的安全问[1]题和经济损失 。因此水利工程对地下水渗流破坏的
预防非常关键。
[0003] 土石坝作为世界坝工建筑中应用最广泛，发展最快的一种形式，结构简单，取材方便。土石坝筑坝的结构形式和材料特性，使得土石坝成为在各种形式的档水
建筑物中，与渗流关系最为密切的一种。渗流问题直接关系到土石坝工程的稳定和安全，国内外许多溃决或破坏的土石坝工程，都是由渗流引发的，所以对土石坝渗流场合理准确地模拟分析是很有必要的。本发明主要解决复杂地质条件下准确对土石坝渗流场模拟分析方面的问题。
[0004] 国内外在渗流场模拟研究中多采用达西定律分析求解渗流场，或利用
温度场与渗流场等效的原理进行渗流计算，应用限制条件多，且无法有效准确的确定浸润面的
位置。土石坝渗流属于有浸润面的自由渗流，如何简捷、准确地确定浸润面的位置是自由渗流场分析的关键及难点。目前广泛采用的方法是有限元固定网格法，包括Bathe和[2] [3] [4,5]Khoshgoftaar 提出的单元渗透矩阵调整法，Desai 提出的剩余流量法，张有天等 提[6]
出的初流量法以及Jiang等 采用的变分不等式方法等，该类方法虽对复杂区域的适应性较好，但由于浸润面位置待定，需要
迭代确定，有限元固定网格法影响了渗流场的计算精[7]
度。同时，其计算工作量大，容易产生迭代计算不收敛现象。Jie等 提出运用无网格自然单元法求解带浸润面的渗流场，该法虽然计算简便，
精度较高，但求解效率较低且容易产生[8] [9] [10] [11]
积分误差。杨海英 、钟小彦 、宋永占 、张瑜 等曾运用有限体积法和
流体体积函数(volume of fluid，VOF)法分别对闸坝坝基、
辐射井、重
力坝坝基、堤防等的渗流场进行分析。然而目前大坝渗流研究中大多对地质条件和坝体细节进行了较大程度的简化，尤其是在坝基方面，未能全面考虑复杂地质条件下的不同岩层地质信息、不良地质体信息以及渗控措施等各种因素对渗流场的影响；此外，运用有限体积法和VOF法对土石坝三维复杂自由渗流问题的研究鲜见报道。
[0005] 其次，排水孔的模拟一直是渗流场模拟中的技术难题，目前常采用的方法有王镭[12] [13] [14]等 提出的排水子结构法，Zhan等 提出的近似解析解结合有限元
算法，王恩志 提[15]
出的“以管代孔”、“以缝代井列”法，胡静 提出的空气单元法等。但这些方法大多都是基于等效原理提出来的，在理论和应用方面都有一定的不足。
[0006] 综上诉述，国内外研究的重点在于渗流的计算方法和具体工程应用等方面。在土石坝渗流模拟过程中，鉴于地质条件和裂隙的复杂及不确定性，很难能够建立准确的地质及裂隙模型，因此常对实际地质条件进行比较大的简化处理，这显然会大大降低模拟结果的精度。此外，目前渗流计算中常采用的方法是
有限元分析方法，将温度场等效为渗流场，或采用达西定律求解渗流场，应用限制条件多，且无法有效准确的确定浸润面的位置。对于无压渗流的浸润自由面的模拟多采用有限元固定网格法，该方法虽对复杂区域的适应性较好，但由于浸润自由面位置是待求的，事先未知，需要迭代确定，有限元固定网格法影响了渗流场的计算精度。为此，本发明提出了建立包括不良地质体与裂隙的枢纽区三维地质精细模型的方法，以及基于有限体积法和VOF法的土石坝渗流计算的数学模型，实现了对土石坝复杂自由渗流场的数值模拟。
[0007] [参考文献]
[0008] [1]富强.基于先进数值模拟方法的复杂条件下土石坝渗流研究[D].浙江：浙江大学，2010。
[0009] [2]Bathe K J,Khoshgoftaar M R.Finite element free surface analysis without mesh iteration[J].Int J Numer Anal Met,1979,3:13-22。
[0010] [3]Desai C S.Finite element residual schemes for unconfined flow[J].Int J Numer Meth Eng,1976,10:1415-1418。
[0011] [4]张有天，陈平，王镭.有自由面渗流分析的初流量法[J].水利学报，1988，(8):18-26。
[0012] [5]潘树来，王全凤，俞缙.利用初流量法分析有自由面渗流问题之改进[J].岩土工程学报，2012，34(2):202-209。
[0013] [6]Jiang Q H,Ye Z Y,Yao C,et al.A new variational inequality formulation for unconfined seepage flow through fracture networks[J].Sci China Tech Sci,2012,55:3090-3101。
[0014] [7]Jie Y X,Liu L Z,Xu W J,et al.Application of NEM in seepage analysis with a free surface[J].Mathematics and Computers in Simulation,2013,89:23-37。
[0015] [8]杨海英.采用有限体积法进行闸坝地基渗流分析与防渗措施研究[D].西安：西安理工大学,2005。
[0016] [9]钟小彦.基于多孔介质模型和VOF法的渗流场数值模拟[D].西安：西安理工大学,2010。
[0017] [10]宋永占.基于有限体积法的砼重力坝坝基防渗排水措施分析研究[D].兰州：兰州理工大学,2011。
[0018] [11]张瑜.基于Fluent的堤防工程渗流场数值分析[D].太原：太原理工大学,2013。
[0019] [12]王镭，刘中，张有天.有排水孔幕的渗流场分析[J].水利学报，1992，(4):15-20。
[0020] [13]Zhan M L,Su B Y.New method of simulating concentrated drain holes in seepage control analysis[J].Journal of Hydrodynamics.Ser.B,1999,3:27～35。
[0021] [14]王恩志，王洪涛，王慧明.“以缝代井列”——排水孔幕模拟方法探讨[J].
岩石力学与工程学报，2002，21(1):98-101。
[0022] [15]胡静，陈胜宏.渗流分析中排水孔模拟的空气单元法[J].岩土力学，2003，24(2):281-283。
发明内容
[0023] 为克服
现有技术的不足，准确地获悉复杂地质条件下土石坝枢纽区的渗流场情况，本发明提出一种既能精确反映地质信息又能准确模拟出土石坝枢纽区渗流场的数值模拟方法。
[0024] 本发明一种基于三维精细地质模型预测土石坝枢纽区渗流场的方法，包括以下步骤：
[0025] 步骤一、建立包括不良地质体与裂隙的枢纽区三维地质精细模型，包括：
[0026] 将水电工程地质数据分为确定性数据和统计性数据；对于确定性地质数据，采用三维混合数据结构、地质构造曲面和地质体的NURBS(Non-Uniform Rational B-Splines)构建枢纽区岩层、不良地质体和渗控结构的三维地质模型；对于统计性地质数据，采用蒙特卡罗(Monte Carlo)法构建岩体三维
裂隙网络模型；针对三维地质模型与岩体三维裂隙网络模型，基于布尔运算分析模型之间的空间位置关系，对重叠的地质体进行布尔差运算，实现岩层、不良地质体以及裂隙网络模型在空间位置上的耦合，从而获得三维精细地质模型；
[0027] 步骤二、基于三维精细地质模型建立三维渗流计算
网格模型：通过计算
流体力学CFD(Computational Fluid Dynamics)
软件，将包含
地层信息、不良地质体信息和裂隙信息的地质数据耦合转化为计算流体力学模型数据，同时，采用贴体网格和局部加密网格划分法将三维精细地质模型转化为三维渗流计算网格模型，用以实现实际地质信息在计算流体力学CFD软件计算模型中的精确表达；
[0028] 步骤三、基于三维精细地质模型的土石坝渗流场数值模拟，包括：
[0029] 3-1、建立用以土石坝渗流场数值模拟的数学模型，该数学模型包括连续性方程、耦合VOF(Volume of Fluid)法的水气两相流N-S方程、水气两相混合
密度方程和混合速度方程；混合速度方程和混合密度方程用以表达渗流场中水气两相的混合速度和混合密度，混合速度方程和混合密度方程均通过渗流阻力源项耦合于水气两相流N-S方程；其中：连续性方程：
[0030]
[0031] 水气两相流N-S方程：
[0032]
[0033] 式(1)和式(2)中：u是微元断面的平均速度，u与孔隙流体真实速度u'的关系是u＝nu'，n为多孔介质的孔隙率，单位是m/s；t为时间，单位是s；f为单位
质量力，单位3
是N；p为流体微元体上的压力，单位是Pa；ρ为流体密度，单位是m/s；v为流体运动粘滞
2
系数，单位是m/s；S为渗流阻力源项，等于惯性损失项Si与粘性损失项Sv之和；设Si＝0，粘性损失项Sv的表达式为：
[0034]
[0035] 式(3)中： 为粘性阻力系数，单位是1/m2，其表达式为：
[0036]
[0037] 式(4)中：K为多孔介质的渗透系数，单位是m/s；耦合VOF法后，水气两相共用一套方程，式(1)至式(4)中的ρ和ν不再是常量，而是采用水气两相的混合密度方程和混合速度方程确定的变量，水气两相混合密度方程和混合速度方程的表达式分别为：
[0038] 混合密度方程：ρ＝Fwρw+(1-Fw)ρa (5)
[0039] 混合速度方策：ν＝Fwνw+(1-Fw)νa (6)
[0040] 式(5)和式(6)中：ρw、ρa分别为水、气两种流体的密度，单位是kg/m3；νw、νa2
分别为水、气两种流体的运动黏性系数，单位是m/s；Fw为水的体积分数，单位是无因次；
[0041] 3-2、确定三维渗流计算网格模型的边界条件，包括：进口边界条件：上游库水淹没部分，取为定
水头压力流入边界，根据实际工程确定坝前上游水位，从而确定出上游定水头压力，根据上游库水位设置进口水气两相组分比；出口边界条件：下游边界取为流出边界，将位于下游水位以下部分取为定水头压力流出边界，下游水位以上部分取为自由出水边界，根据下游库水位设置出口水气两相组分比；固体壁面边界条件：按固壁定律处理，所有固壁处的
节点均采用无滑移条件处理；
[0042] 3-3、根据步骤3-1确定的连续性方程，耦合VOF法的水气两相流N-S方程，水气两相混合速度方程和混合密度方程，步骤3-2设置的边界条件，以及根据实际工程确定的渗透参数，通过有限体积法进行离散，利用SIMPLE(Semi-Implicit Method for Pressure Linked Equations)算法进行求解，实现基于三维精细地质模型的土石坝渗流场数值模拟，从而预测土石坝枢纽区渗流场。
[0043] 本发明中提供的一种确定土石坝坝体浸润面的方法是利用上述基于三维精细地质模型预测土石坝枢纽区渗流场的方法获得的渗流场，并采用VOF法确定土石坝坝体浸润面，方法如下：
[0044] VOF法引入单流
体模型来处理
多相流问题，对于水气两相流流场，设在同一单元中的水、气或者二者的混合体具有相同的速度，即服从同一组动量方程，将水气两相的体积函数在整个渗流场中都作为单独变量；在任意一个单元中，水和气的体积分数之和等于1，如果Fw表示水的体积分数，则气的体积分数为1-Fw；当Fw＝1，表示该单元全部为水相所占据；当Fw＝0，表示该单元全部为气相所占据；当0
[0045] 体积率函数：
[0046]
[0047] 式(7)中：t为时间，单位是s；ρ为流体密度，单位m3/s；ν为流体运动粘滞系数，2
单位是m/s；
[0048] 通过求解体积率函数可以得到空间各处水气两相各自的体积分数，在界面区域内采用分段线性插值的几何重建方法得到水气分界面；土石坝中的渗流是具有浸润面的自由渗流，而渗流浸润面本质上是水气两相流的分界面，因此可以采用VOF法来获得碾压
混凝土坝浸润面。
[0049] 与现有技术相比，本发明的有益效果是：
[0050] (1)采用NURBS构造技术建立包括各类不良地质体和渗控结构的三维精细地质模型，而且采用蒙特卡罗技术对三维裂隙网络进行构建，并通过布尔运算分析实现地质模型与裂隙网络的耦合，实现了对坝基复杂地质条件的精细模拟。
[0051] (2)通过地质模型数据与CFD模型数据之间的耦合转化，建立CFD模型，实现了真实复杂地形在CFD软件中的精确表达，弥补了以往建模中网格精确性不够的局限。
[0052] (3)基于三维精细地质模型，采用计算流体力学方法(CFD)，建立等效连续性多孔介质模型和三维渗流场数值模拟模型，实现了复杂地质条件下的土石坝枢纽区三维渗流场的准确模拟。
[0053] (4)采用VOF法模拟确定坝体渗流场浸润面，实现了直接模拟渗流浸润面，模拟结果更加精确。
附图说明
[0054] 图1是本发明基于三维精细地质模型的土石坝渗流场数值模拟方法
流程图；
[0055] 图2是本发明中工程地质三维精细建模总体结构；
[0056] 图3是本发明中SIMPLE算法流程图。
具体实施方式
[0057] 下面结合附图对本发明技术方案作进一步详细描述。
[0058] 本发明构建了一种能够耦合三维精细地质模型的土石坝渗流场数值模拟方法，其既能精确反映地质信息又能准确模拟出土石坝枢纽区渗流场，达到准确预测土石坝枢纽区渗流场的目的，为土石坝渗流控制提供指导。
[0059] 本发明的基于三维精细地质模型预测土石坝枢纽区渗流场的方法主要包括包含地层、不良地质体和裂隙的枢纽区三维地质精细模型建模方法，基于三维精细地质模型的土石坝渗流场数值模拟方法，以及运用VOF法求解土石坝坝体浸润面方法。总体技术流程如图1所示，主要包括以下三大部分：(1)三维地质建模：首先，根据水利枢纽坝区多源地质数据和大坝设计资料，建立坝体、地层、不良地质体、渗控结构、裂隙等五大模型；然后，基于布尔运算分析地质模型与裂隙模型之间的空间位置关系，对重叠的地质体进行布尔差运算，得到三维精细地质统一模型；最后，基于三维精细地质统一模型，运用贴体网格和局部加密网格划分法建立三维渗流计算网格模型。(2)三维渗流场数值模拟：基于三维渗流计算网格模型，根据工程实际情况设置模型的渗透参数及边界条件，采用三维渗流计算数学模型和三VOF数学模型，计算得到土石坝枢纽区渗流场。然后通过对土石坝枢纽区渗流场的分析得到枢纽区的水头、水力梯度和扬压力的分布，以及渗漏量。(3)坝体浸润面求解分析：根据部分(2)得到的土石坝枢纽区渗流场，采用VOF法计算水气两相的体积分数，得到水气两相的交界面，进而得到土石坝坝体的浸润面。
[0060] 本发明的最佳实施方式如下：
[0061] 步骤一、建立包括不良地质体与裂隙的枢纽区三维地质精细模型，[0062] 包括：
[0063] (1)将水电工程地质数据分为确定性数据和统计性数据；对于确定性地质数据，采用三维混合数据结构、地质构造曲面和地质体的NURBS(Non-Uniform Rational B-Splines)构造技术构建枢纽区岩层、不良地质体和渗控结构的三维地质模型；其中，不良地质体包括
挤压带、
破碎带、溶蚀带和软弱夹层中的一种或几种。
[0064] (2)对于统计性地质数据，采用蒙特卡罗(Monte Carlo)法构建岩体三维裂隙网络模型；三维裂隙网络模型构建主要包括以下7个步骤：1)确定
采样区域，根据实际采样区域内裂隙的产状进行有效的分组；2)确定每组裂隙的数量和密度；3)确定每组裂隙迹长、间距、产状等几何参数所服从的概率模型；4)利用蒙特卡罗方法生成一系列模拟参数，并验证生成的结果满足预期的概率模型；5)生成每组裂隙的三维裂隙网络模型；6)剖切裂隙模型，对每组裂隙迹长、间距的参数进行校核，保证其与预设值向一致；7)生成初步三维网络模型。
[0065] (3)针对三维地质模型与岩体三维裂隙网络模型，基于布尔运算分析地质模型与裂隙模型之间的空间位置关系，对重叠的地质体进行布尔差运算，实现岩层、不良地质体以及裂隙网络模型在空间位置上的耦合，从而获得三维精细地质模型。工程地质与裂隙的三维精细建模总体结构如图2所示，依据原始勘测数据、2D地质剖面数据和工程设计数据，采用混合数据结构面和向对象技术，进行几何建模，建模过程主要分为以下三个部分：(1)建立坝体、渗控结构、廊道等人工对象的模型，然后将它们组合在一起，得到工程建筑物模型；(2)建立地形、地层、裂隙、不良地质体等自然地质对象的模型，然后将它们组合在一起，得到三维地质模型，最后再运用原始勘测数据、2D地质剖面数据和工程设计数据进行改进的地质趋势面分析，从而对建立的三维地质模型进行补充验证；(3)基于布尔运算分析地质模型与工程建筑物模型之间的空间位置关系，对重叠部分进行布尔差运算，从而得到土石坝枢纽区三维统一模型。
[0066] 步骤二、基于三维精细地质模型建立三维渗流计算网格模型：
[0067] 基于建立的包含复杂地层信息的三维地质模型，通过计算流体力学CFD(Computational Fluid Dynamics)软件，将包含地层信息、不良地质体信息和裂隙信息的地质数据耦合转化为计算流体力学模型数据，同时，采用贴体网格和局部加密网格划分法将三维精细地质模型转化为三维渗流计算网格模型，用以实现实际地质信息在计算流体力学CFD软件计算模型中的精确表达。
[0068] 步骤三、基于三维精细地质模型的土石坝渗流场数值模拟：
[0069] 包括：
[0070] 3-1、建立用以土石坝渗流场数值模拟的数学模型，该数学模型包括连续性方程、耦合VOF(Volume of Fluid)法的水气两相流N-S方程、水气两相混合密度方程和混合速度方程；混合速度方程和混合密度方程用以表达渗流场中水气两相的混合速度和混合密度，混合速度方程和混合密度方程均通过渗流阻力源项耦合于水气两相流N-S方程；
[0071] 其中：
[0072] 连续性方程：
[0073]
[0074] 水气两相流N-S方程：
[0075]
[0076] 式(1)和式(2)中：u是微元断面的平均速度，u与孔隙流体真实速度u'的关系是u＝nu'，n为多孔介质的孔隙率，单位是m/s；t为时间，单位是s；f为单位质量力，一般仅3
为重力，单位是N；p为流体微元体上的压力，单位是Pa；ρ为流体密度，单位是m/s；ν为
2
流体运动粘滞系数，单位是m/s；S为渗流阻力源项，等于惯性损失项Si与粘性损失项Sv之和；对于流速较低的
层流状态可忽略惯性损失项，即认为Si＝0；对于流速较高的流动可忽略粘性损失项，即认为Sv＝0。由于本发明中研究的碾压混凝土坝渗流场中水的渗透速度较低，属层流状态，因此可以忽略惯性阻力项，设Si＝0，粘性损失项Sv的表达式为：
[0077]
[0078] 式(3)中： 为粘性阻力系数，单位是1/m2，其表达式为：
[0079]
[0080] 式(4)中：K为多孔介质的渗透系数，单位是m/s；
[0081] 耦合VOF法后，水气两相共用一套方程，式(1)至式(4)中的ρ和ν不再是常量，而是采用水气两相的混合密度方程和混合速度方程确定的变量，水气两相混合密度方程和混合速度方程的表达式分别为：
[0082] 混合密度方程 ρ＝Fwρw+(1-Fw)ρa (5)
[0083] 混合速度方策 ν＝Fwνw+(1-Fw)νa (6)
[0084] 式(5)和式(6)中：ρw、ρa分别为水、气两种流体的密度，单位是kg/m3；νw、νa分别为水、气两种流体的运动黏性系数，单位是m2/s；Fw为水的体积分数，单位是无因次。
[0085] 3-2、确定三维渗流计算网格模型的边界条件，包括：
[0086] 进口边界条件：上游库水淹没部分，取为定水头压力流入边界，根据实际工程确定坝前上游水位，从而确定出上游定水头压力，根据上游库水位设置进口水气两相组分比；
[0087] 出口边界条件：下游边界取为流出边界，将位于下游水位以下部分取为定水头压力流出边界，下游水位以上部分取为自由出水边界，根据下游库水位设置出口水气两相组分比；
[0088] 固体壁面边界条件：按固壁定律处理，所有固壁处的节点均采用无滑移条件处理。
[0089] 3-3、根据步骤3-1确定的连续性方程，耦合VOF法的水气两相流N-S方程，水气两相混合速度方程和混合密度方程，步骤3-2设置的边界条件，以及根据实际工程确定的渗透参数(包括坝体、地层及不良地质体的透水率和渗透系数)，通过有限体积法进行离散，利用SIMPLE(Semi-Implicit Method for Pressure Linked Equations)算法进行求解，实现基于三维精细地质模型的土石坝渗流场数值模拟，从而预测土石坝枢纽区渗流场。SIMPLE算法流程如图3所示，主要包括以下9个步骤：(1)边界条件设置完毕，准备计算；
(2)假设一个速度分布u0，v0，用于计算首次迭代时的动量离散方程的系数和常数项(3)假设一个压力场，即给定压力猜测值p*；(4)根据当前速度场和压力场，计算动量离散方程等方程中的系数和常数项 (5)求解动量离散方程，得到速度u*,v*；(6)根据速度u*,v*，求解压力修正方程，得到修正后的压力p；(7)根据修正后的压力改进速度，得到修正后速度u、v；(8)利用修正后的速度场求解所有其他的离散化输运方程，得到动量离散方程等方程中的系数和常数项 (9)判定收敛与否，若收敛，则计算结束；若未收敛，则令p*＝p，u*＝u，V*＝v， 返回步骤(4)中重新计算，循环迭代，直到收敛为止。
[0090] 本发明中利用上述基于三维精细地质模型预测土石坝枢纽区渗流场的方法获得的渗流场，并采用VOF法确定土石坝坝体浸润面，
[0091] 由Hirt和Nichols提出的VOF法是一种处理自由表面的有效方法。土石坝坝体中的渗流是具有浸润面的无压渗流，而渗流浸润面本质上是水气两相流的交界面，因此可以采用VOF法来获得浸润面。该方法如下：
[0092] VOF法引入单流体模型来处理多相流问题，对于水气两相流流场，设在同一单元中的水、气或者二者的混合体具有相同的速度，即服从同一组动量方程，将水气两相的体积函数在整个渗流场中都作为单独变量。在任意一个单元中，水和气的体积分数之和等于1。在一个单元中，如果Fw表示水的体积分数，则气体的体积分数为1-Fw。Fw有3种情况：Fw＝0、0
[0093] 当Fw＝0，表示该单元全部为气相所占据；
[0094] 当0
[0095] 当Fw＝1，表示该单元全部为水相所占据；
[0096] 其它水气共有的变量，如水气两相共有的流速、流量、压力都可以采用体积函数的加权平均值来表示。
[0097] 体积率函数：
[0098]
[0099] 式(7)中：t为时间，单位是s；ρ为流体密度，单位m3/s；ν为流体运动粘滞系数，单位是m2/s。
[0100] 通过求解体积率函数可以得到空间各处水气两相各自的体积分数，在界面区域内采用分段线性插值的几何重建方法得到水气分界面，即浸润面。土石坝中的渗流是具有浸润面的自由渗流，而渗流浸润面本质上是水气两相流的分界面，因此可以采用VOF法来获得碾压混凝土坝浸润面。
[0101] 尽管上面结合附图对本发明进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨的情况下，还可以做出很多
变形，这些均属于本发明的保护之内。
