基于水泥浆液触变性和VOF求解扩散锋面的帷幕灌浆数值模拟方法与流程

本发明属于水利水电工程中大坝基岩帷幕灌浆数值模拟领域，更具体的说，是涉及一种基于水泥浆液触变性和VOF求解扩散锋面的帷幕灌浆数值模拟方法。

背景技术：

帷幕灌浆是大坝基岩防渗处理的主要手段，对保证水工建筑物的安全运行起着重要作用。由于防渗帷幕完全位于地下，具有隐蔽性，加之地质条件复杂，导致浆液在地层中的运动过程也极为复杂，为精确描述浆液运动规律带来了困难，对施工后帷幕的连续性和完整性灌浆效果难以进行评估。因此对帷幕灌浆过程进行数值模拟和分析，揭示地质情况和浆液性能之间的关系，确定帷幕灌浆浆液锋面搭接闭合情况，对指导灌浆设计和施工具有重要的理论意义和工程应用价值。

国外关于灌浆数值模拟的研究较少^[1]，其方法多是采用UDEC软件模拟隧道裂隙注浆。UDEC是一个处理不连续介质的二维离散元程序，故数值分析仅局限于二维，同时也不能考虑灌浆孔附近的地质因素和水泥浆液的流变性。

国内灌浆模拟研究中，李粮纲等^[2]采用基于UDEC离散单元程序的数值模拟方法，对深矿钻孔帷幕灌浆压力的扩散分布及浆液的流动扩散进行模拟；袁敬强等^[3]基于散体介质理论的颗粒流方法，运用PFC^2D软件对软弱地层注浆过程进行了细观力学模拟，研究了不同注浆要素对浆液扩散半径和注浆类型的影响规律；刘健等^[4]采用comsol流体力学模块研究了水泥浆液在静水和动水条件下平面裂隙中的扩散规律；韩磊等^[5]用FLAC中两相流模块对某土石坝一5m灌浆孔段的灌浆过程进行了二维模拟分析。上述方法均是二维数值计算，未考虑地质条件和浆液流变性质。王晓玲等^[6]采用混合模型对宾汉姆流体浆液在采空区中的紊流扩散进行了模拟计算。邓韶辉等^[7]将水泥浆液看作单一均质宾汉姆流体对浆液在大坝基岩中的渗透灌浆进行了模拟研究。本课题组虽实现了灌浆三维数值研究，但将浆液看作宾汉姆流体，不如触变性模型能对浆液性质进行更符合的流变性质表达，且均为涉及浆液扩散锋面的详细计算方法。

综上所述，国内外现有研究多数是二维数值计算，且未考虑地质条件和流体性质；极个别三维数值研究中将水泥浆液看作宾汉姆流体，不如触变性模型更精确地表达水泥浆液流变性能，且均为涉及浆液扩散锋面的计算方法。目前基于水泥浆液触变性模型和VOF计算确定浆液扩散锋面的考虑复杂地质条件的三维灌浆数值研究还未见报道。

[1]Yong Xu,Weihong Peng.Research on Multiple Holes Grouting of Fractured Rock Mass[J].Applied Mechanics and Materials,2013,547-551.

[2]李粮纲，唐平，何维山，等.深矿钻孔帷幕灌浆的数值模拟与应用[J].探矿工程，2010，37(12)：36-40.

[3]袁敬强，陈卫忠，谭贤君，等.软弱地层注浆的细观力学模拟研究[J].岩土力学，2011，32(S2)：653-659.

[4]刘健，刘人太，张霄，等.水泥浆液裂隙注浆扩散规律模型试验与数值模拟[J].岩石力学与工程学报，2012，31(12):2445-2452.

[5]韩磊，陈建生，陈亮.帷幕灌浆扩散半径及数值模拟的研究[J].岩土力学，2012，33(7):2235-2240.

[6]王晓玲,王青松,周正印,敖雪菲.南水北调工程采空区三维宾汉姆流体紊流模拟[J].水利学报,2013,11:1295-1302.

[7]邓韶辉,王晓玲,敖雪菲,任炳昱,李瑞金.大坝基岩帷幕宾汉姆浆液灌浆的三维数值模拟[J].水利学报,2016,02:165-172+179.

技术实现要素：

本发明的目的是为了克服现有技术中的不足，提供一种基于水泥浆液触变性和VOF求解扩散锋面的帷幕灌浆数值模拟方法，实现了复杂地质条件下基于触变性模型和VOF计算确定浆液扩散锋面的三维灌浆数值模拟，准确获悉复杂地质条件下基岩帷幕灌浆的灌后宏观效果，为大坝帷幕灌浆的灌前预测或灌后分析提供一种精确的数值手段。

本发明的目的是通过以下技术方案实现的。

本发明的基于水泥浆液触变性和VOF求解扩散锋面的帷幕灌浆数值模拟方法，包括以下步骤：

(1)建立三维地质模型和灌浆孔模型，耦合三维地质模型和灌浆孔模型构建三维统一模型；

(2)将步骤(1)中的三维统一模型导入CFD计算软件，划分计算网格；

(3)建立水泥浆液触变性流体模型，VOF法捕捉浆液扩散锋面；

(4)实现复杂地质条件下的考虑浆液触变性和VOF计算浆液扩散锋面的大坝基岩帷幕灌浆三维数值模拟。

所述步骤(1)中三维地质模型的建立：对地表进行地形建模，利用勘测的原始资料，以TIN模型方式将这些离散的地形等高线按一定的规则划分形成互联的不规则三角形网络，以TIN模型为中介，用NURBS技术对地表TIN模型进行三维地形重构，得到三维地形可视化模型。

所述步骤(1)中灌浆孔模型的建立：通过建立中尺度灌浆孔的数据变量集，依据“点‐线‐面‐体”建模方法，实现灌浆孔三维模型的构建。

所述步骤(1)中三维统一模型的建立：将三维地质模型和灌浆孔模型进行布尔计算，在三维地质模型的基础上将灌浆数据向三维空间模型转化，使三维地质模型和灌浆孔模型正确耦合。

所述步骤(3)中水泥浆液触变性流体模型为：

式中，为表观粘度，η₀为触变性流体的塑性粘度，τ₀为触变性流体的屈服应力，λ₀是材料的结构参数，为剪切应变率，t为时间，α为表示材料属性的无量纲参数。

所述步骤(3)中VOF法捕捉浆液扩散锋面具体为VOF方法通过求解各时刻计算网格单元中流体体积比函数来确定浆液扩散锋面的位置和形状，所述流体体积比函数为：

式中，F为流体体积比函数，x_i表示x,y,z三个不同方向，u_i表示x,y,z三个不同方向的速度。

所述步骤(4)中实现复杂地质条件下的考虑浆液触变性和VOF计算浆液扩散锋面的大坝基岩帷幕灌浆三维数值模拟采用连续性方程、耦合VOF法的气液两相流N-S方程、气液两相混合密度方程和混合速度方程、触变性流体本构方程。

与现有技术相比，本发明的技术方案所带来的有益效果是：

(1)本发明将三维地质模型和灌浆孔模型进行耦合，为整个灌浆工程分析提供一个新的视角，对灌浆孔的空间分布和地质条件形成整体认识，为灌浆数值分析提供精细的三维计算几何模型；

(2)本发明通过三维地质模型、灌浆孔模型数据与CFD模型数据之间的接口转化，建立CFD计算网格模型，实现真实复杂地质条件在CFD软件中的精确表达，弥补了以往CFD建模中通过坐标绘制网格模型而使网格精确性不足的局限；

(3)本发明用触变性模型表达水泥流变性质，考虑了灌浆过程中的剪切率和时间两相因素，比宾汉姆模型更能反映水泥浆液的时变性，采用VOF法捕捉帷幕浆液扩散锋面，具有计算工作量小、数据存储量小及对计算易收敛的特点；

(4)本发明建立的基于复杂地质条件考虑浆液触变性和VOF计算浆液锋面的帷幕灌浆数值模拟方法，实现了坝基帷幕灌浆过程中水泥浆液的多孔分序灌浆精细模拟，模拟条件更加符合现场实际灌浆条件，可对施工后帷幕的连续性和完整性灌浆效果进行有效评估。

附图说明

图1是耦合三维地质模型和灌浆孔模型构建三维统一模型的方法流程图；

图2是利用CFD技术实现触变性浆液大坝基岩帷幕灌浆三维数值模拟技术流程图。

具体实施方式

下面结合附图对本发明作进一步的描述。

本发明提出一种基于复杂地质条件考虑浆液触变性和VOF计算浆液锋面的帷幕灌浆三维数值模拟方法。首先耦合三维地质模型和灌浆孔模型构建三维统一模型；然后通过数据接口，将三维统一模型导入CFD软件，采用多面体网格划分技术，对三维统一模型划分精细网格；最后对水泥浆液建立能真实反映其流变性能的触变性流体模型，利用VOF法计算捕捉浆液扩散锋面，实现复杂地质条件下的考虑浆液触变性和VOF计算浆液扩散锋面的大坝基岩帷幕灌浆三维数值模拟。本发明的基于水泥浆液触变性和VOF求解扩散锋面的帷幕灌浆数值模拟方法的技术流程如图1和图2所示，具体地讲涉及包括耦合三维地质模型和灌浆孔模型的构建方法、触变性模型模拟表达水泥浆液、VOF求解浆液扩散锋面方法和基于以上方法的大坝基岩帷幕灌浆数值模拟，最佳实施方式如下：

一、建立三维地质模型和灌浆孔模型，耦合三维地质模型和灌浆孔模型构建三维统一模型，包括：

(1)建立三维地质模型，其实现过程可简述如下：首先对地表进行地形建模，利用等高线等原始资料，以TIN模型方式将这些离散的地形等高线按一定的规则划分形成互联的不规则三角形网络，以TIN模型为中介，用NURBS技术对地表TIN模型进行三维地形重构，得到三维地形可视化模型，三维地形可视化模型也作为所有地质体运算操作的空间约束和基础。①根据每一个地质体自身的原始地质数据，利用上述的三维建模技术，分别建立地质体的上、下结构面S1和S2；②利用结构面S1和S2的空间位置关系，分别对受体—地表三维模型进行布尔操作，形成地质体T1③最终对所有不良地质进行空间拓扑分析。

(2)建立灌浆孔模型：中尺度灌浆工程结构物是指实际灌浆中分布在岩体空间内的灌浆载体，通过建立中尺度灌浆孔的数据变量集，依据“点-线-面-体”建模方法，实现灌浆孔三维模型的构建。综合灌浆孔孔口位置信息(x,y,z)，段长信息L、孔口直径D信息、倾角β信息，实现灌浆孔的三维可视化建模。①用(x,y,z)坐标生成对应的三维空间点P1；②P1为圆心，D/2为半径，利用犀牛Circle命令(函数ON_Circle::ON_Circle)生成曲线C1；③Z方向作为轴线，L为长度，利用犀牛ExtrudeCurveStraight命令(函数RhinoExtrudeCurveStraight)实现曲线C1的拉伸放样，从而形成曲面S1；④用犀牛Cap命令(函数Cap)使曲面S1形成空间Z方向约束，完成灌浆孔从数字型数据向三维空间模型的转化，以上步骤①-④用程序可自动完成。

(3)耦合三维地质模型和灌浆孔模型构建三维统一模型：在三维地质模型的基础上，实现灌浆数据向三维空间模型的转化，也就是将三维地质模型和灌浆孔模型进行布尔计算，在三维地质模型的基础上将灌浆数据向三维空间模型转化，使三维地质模型和灌浆孔模型正确耦合，将为整个灌浆工程分析提供一个新的视角。采用该方法建立的三维统一模型，为本研究的大坝基岩帷幕灌浆CFD分析提供精细的计算几何模型基础。

二、将三维统一模型导入CFD计算软件，划分计算网格具体为：

(1)将三维统一模型以CFD软件可以接受兼容的.stl格式导出，通过CFD(计算流体力学)软件的数据接口，将包含真实地质地形和灌浆孔数据的.stl格式文件导入到CFD软件中，实现三维统一模型在CFD软件中的真实表达。

(2)基于导入的三维统一模型，对不同地质区域采用多面体网格划分技术建立模型：在全局选用较大的尺寸，生成与计算区域边界重合的、疏密程度不均匀的网格，使得网格的边界与计算区域边界一一对应；在灌浆孔局部区域进行设置较小尺寸，使网格加密，尽可能保留维持灌浆孔的几何形状，为帷幕灌浆的CFD计算提供精细网格模型。

三、建立水泥浆液触变性流体模型，VOF法捕捉浆液扩散锋面具体为：

(1)水泥浆液触变性流体模型：新鲜水泥浆表现为非牛顿粘性流体。在瞬变流情况下，表观粘度是剪切率和时间的函数。在灌浆工程中，灌注浆液随着灌浆时间性质发生变化，可用瞬变流情况下的触变模型表示。

触变性流体的本构关系为：

式中：τ为偏应力张量，Pa；为表观粘度，是λ和的函数，Pa·s；λ是材料的结构参数，为剪切应变率；

为避免数值耗散，对触变性水泥浆液表观粘度改进为：

式中：为表观粘度，Pa·s；η₀为触变性流体的塑性粘度，Pa·s；τ₀为触变性流体的屈服应力，Pa；λ₀是材料的结构参数；为剪切应变率；t为时间；α为表示材料属性的无量纲参数。

将以上触变性流体表观粘度(式(2))进行编程，开发改进到CFD软件中进行计算，保证了流体性能的真实表达。

(2)浆液扩散锋面是气液两相流的分界面，VOF模型追踪气液交界面，习惯上用一个单场模型，即气液两相共用一个速度场和密度场。基本思路为：通过求解一套动量方程和连续方程模拟两种互不掺混的流体的运动，追踪每种流体所占的体积，在每一控制体中所有流态的体积函数为整数1。VOF方法通过求解各时刻计算网格单元中流体体积比函数(流体体积量的变化和网格单元本身体积的比值函数)来确定浆液扩散锋面的位置和形状，所述流体体积比函数为：

式中，F为流体体积比函数，x_i表示x,y,z三个不同方向，u_i表示x,y,z三个不同方向的速度。

对于某时刻任意网格单元中的指定相流体而言，若F＝0，表明网格单元被其它相流体所充满，即不包含指定相流体；若0<F<1，表明网格单元内存在多相物质的交界面；若F＝1，表明网格单元被指定相流体全部充满，即是指定相流体单元。

四、实现复杂地质条件下的考虑浆液触变性和VOF计算浆液扩散锋面的大坝基岩帷幕灌浆三维数值模拟具体为：

基本方程包括连续性方程、耦合VOF法的气液两相流N-S方程、气液两相混合密度方程和混合速度方程、触变性流体本构方程；混合速度方程和混合密度方程用以表达渗流场中气液两相的混合速度和混合密度，混合速度方程和混合密度方程均通过渗流阻力源项耦合于气液两相流N-S方程，触变性流体的表观粘度η关联N-S方程的运动粘度ν用以表达水泥浆液的触变性能，耦合VOF法的N-S方程用以求解浆液扩散锋面。其中：

连续性方程：

气液两相流N-S方程：

式(4)和式(5)中，φ为多孔介质的孔隙率；u是微元断面的速度矢量，m/s，u与孔隙流体真实速度u′的关系是u＝φu′；t为时间，单位是s；f为单位质量力，单位是N；p为流体微元体上的压力，单位是Pa；ρ为流体密度，单位是m³/s；ν为流体运动粘度，m²/s，ν与流体表观粘度η的关系为ν＝η/ρ；S为渗流阻力源项，等于惯性损失项S_i与粘性损失项S_v之和。对于流速较低的流动可忽略惯性损失项，即认为S_i＝0，粘性损失项S_v的表达式为：

式(6)中，为粘性阻力系数，单位是1/m²，其表达式为：

式(7)中，K为多孔介质的渗透系数，单位是m/s。

耦合VOF法后，气液两相共用一套方程，式(4)至式(7)中的ρ和ν不再是常量，而是采用气液两相的混合密度方程和混合速度方程确定的变量，气液两相混合密度方程和混合速度方程的表达式分别为：

混合密度方程ρ＝F_gρ_g+(1-F_g)ρ_a (8)

混合速度方程ν＝F_gν_g+(1-F_g)ν_a (9)

式(8)和式(9)中：ρ_g、ρ_a分别为浆液、气两种流体的密度，单位是kg/m³；ν_g、ν_a分别为浆液、气两种流体的运动黏性系数，单位是m²/s；F_g为浆液的体积分数，无因次。

引入触变性流体的本构关系，即

式中：τ为偏应力张量，Pa；为表观粘度，是λ和的函数，Pa·s；λ是材料的结构参数，为剪切应变率；

为避免数值耗散，对触变性水泥浆液表观粘度改进为：

式中：为表观粘度，Pa·s；η₀为触变性流体的塑性粘度，Pa·s；τ₀为触变性流体的屈服应力，Pa；λ₀是材料的结构参数；为剪切应变率；t为时间；α为表示材料属性的无量纲参数。

计算条件包括：

(1)进口条件：根据实际灌浆工程的灌浆过程记录数据计算确定水泥浆液进口的速度分布，压力分布以及相应的体积浓度分布，各变量取均匀分布，进口边界条件为灌浆孔边界的浆液的灌入速度，及进口的相应浆液VOF值；根据工程设计方案，灌浆孔分为Ⅰ、Ⅱ、Ⅲ序孔，自编程序实现按时间变化的进口边界条件，进行Ⅰ、Ⅱ、Ⅲ序孔的分序灌浆。

(2)出口条件：考虑在灌浆过程中坝基底部可能存在抬动现象，因此取基岩上表面为压力出口边界。在计算域的边界上，压力应满足Neumann条件，即第二类边界条件。

(3)固壁边界条件：坝基上下游两个侧边界和底部边界均为无滑动壁面边界，按固壁定律处理，在壁面上采用无滑移条件：u＝v＝w＝0，其中u、v、w为x、y、z方向的速度。

(4)根据细砂岩、粉砂岩、粉砂质泥岩和泥质粉砂岩等地质层的地质参数，采用不同的孔隙率和阻力系数。

根据上述确定的连续性方程、耦合VOF法的气液两相流N-S方程、气液两相混合速度方程和混合密度方程、触变性浆液本构方程，以及上述设置的计算条件，通过有限体积法进行离散，利用SIMPLE(Semi-Implicit Method for Pressure Linked Equations)算法进行求解，利用触变性模型对水泥浆液流体性质进行表达，利用VOF法精细捕捉浆液扩散锋面，用CFD技术自编程序实现了坝基帷幕灌浆过程中触变性水泥浆液的多孔分序灌浆精细模拟，分析了灌浆完成后帷幕搭接闭合情况及触变性水泥浆液扩散渗透过程；同时进行了灌浆量和灌浆时间的模拟值与现场实测值的对比研究。

尽管上面结合附图对本发明的功能及工作过程进行了描述，但本发明并不局限于上述的具体功能和工作过程，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可以做出很多形式，这些均属于本发明的保护之内。