CN102799713B - 堆石坝心墙水力劈裂的数值模拟方法 - Google Patents

Abstract

一种堆石坝心墙水力劈裂的数值模拟方法，包括如下步骤：（1）根据室内平面应变试验获取实际心墙土料的应力应变曲线，建立颗粒离散元双轴数值模型；（2）根据颗粒离散元双轴试验确定的细观力学参数，建立心墙骨架颗粒模型；（3）利用连续介质模型建立流体模型，并对其控制方程进行离散；（4）利用心墙骨架颗粒模型和流体模型模拟水力劈裂。本发明根据室内试验获取细观参数，模拟实际心墙土料的力学性能，施加心墙的应力边界条件以及堆石坝上游的实际水压力边界条件，进行数值模拟，从细观角度判定心墙是否发生水力劈裂。根据模拟结果，从施工或设计方面提出改进措施，避免心墙水力劈裂的发生，从而对实际工程具有指导意义。

Description

堆石坝心墙水力劈裂的数值模拟方法

技术领域

本发明属于岩土工程中流固耦合分析技术领域，是将颗粒离散元与无网格法相结合的一种堆石坝心墙水力劈裂的数值模拟方法。

背景技术

在众多坝型中，土石坝因其具有广泛的适应性、可就地取材、施工技术简单且抗震性能良好等优点而成为高坝的首选坝型，在坝工建设中得到了极其广泛的应用。特别是近些年，随着我国西南地区水电开发的进一步深入，在建和拟建的200～300m级高土石坝的数量在逐渐地增加。然而在高土石坝的建设中，心墙堆石坝作为土石坝中常用的坝型之一，还存在着一些亟待解决的重要问题，其中心墙的水力劈裂问题是堆石坝设计和建设中倍受关注的焦点之一。

尽管目前已有不少专家和学者通过试验和有限元数值模拟等手段，从宏观的角度研究水力劈裂的发生、发展过程及其力学机理，但至今尚未形成一致的观点，例如应采用哪一个主应力或哪些主应力的组合与心墙前外水压力进行比较作为水力劈裂的判定准则以及心墙水力劈裂常发生的部位等观点，仍存在许多争议。

发明内容

由于目前从宏观的角度研究心墙水力劈裂尚未得到一致的结论，因此本发明采用颗粒离散元与无网格相结合的方法，从细观角度研究水力劈裂的发生、发展过程，从而实现对心墙水力劈裂机理的分析和研究。

本发明为堆石坝心墙水力劈裂的数值模拟方法，其步骤如下：

(1)根据室内平面应变试验获取实际心墙土料的应力应变曲线，建立颗粒离散元双轴数值模型；

(2)根据颗粒离散元双轴试验确定的细观力学参数，建立心墙骨架颗粒模型；

(3)利用连续介质模型建立流体模型，并对其控制方程进行离散；

(4)利用心墙骨架颗粒模型和流体模型模拟水力劈裂

(4.1)将颗粒与流体间的相互耦合力、颗粒之间的作用力之和作用于所有颗粒，利用离散元法计算并得到心墙骨架颗粒模型内部颗粒的速度、应力分布及渗透特性；初始计算时，颗粒与流体间的相互耦合力为零；

(4.2)提取步骤(4.1)计算得到的渗透特性，利用全局弱形式的无单元EFG法计算并得到渗流流体的速度、压力及颗粒与流体间的相互耦合力；

所述步骤(4.2)的利用全局弱形式的无单元EFG法对流体进行计算，还包括如下子步骤：

(4.2.1)选取每一固体颗粒的中心点作为场节点，并利用这组场节点对流体域进行离散；

(4.2.2)通过离散的节点在流体域内形成用于数值积分的背景网格，将整体计算域划分为若干个积分子域；

(4.2.3)采用无单元EFG法对离散后的控制方程进行数值积分

$\underset{{\mathrm{\Ω}}_j}{\∫}\mathrm{\Φ}[\frac{{\hat{u}}_i^{t+\mathrm{\Δt}}-{\hat{u}}_i^t}{\mathrm{\Δt}}-Q_i]d{\mathrm{\Ω}}_j=0$

式中：为u_i的移动最小二乘近似，Φ是MLS形函数，Ω_j为第j个积分子域；

(4.2.4)通过循环方式，在每一个背景网格内进行数值积分，从而获得整体系统方程组的系数矩阵；

(4.2.5)根据实际的上游水位，给心墙骨架颗粒模型施加已知的压力边界条件；

(4.2.6)求解系统方程组，获得场节点的流速和压力。

(4.3)判断模型内部是否达到平衡状态，如否，回到步骤(4.1)；如是，根据颗粒模型内部是否产生贯穿性裂缝，输出模型内部在指定的心墙上游水压力作用下是否发生水力劈裂的结论，结束模拟过程。

所述步骤(4.2)为经步骤(4.1)循环计算m次后，提取第m次计算得到的渗透特性，利用全局弱形式的无单元EFG法计算并得到渗流流体的速度、压力及颗粒与流体间的相互耦合力；所述m根据计算精度及模拟对象的需要确定。

本发明堆石坝心墙水力劈裂的数值模拟方法的优点是：在采用有限元模拟心墙水力劈裂时，需要根据裂缝的扩展重新划分模型网格，计算工作量十分繁重，而本发明采用颗粒离散元的模拟方法不仅可以避免这一问题，而且还可以方便地模拟心墙水力劈裂发生与扩展的整个过程。心墙水力劈裂模拟涉及到流固耦合问题，本发明根据固体和流体部分各自的计算特点，分别采用了不同的数值模拟方法。针对心墙的固体骨架部分采用颗粒离散元模拟，即用圆形，椭圆形或其他组合形状的颗粒模拟心墙骨架。而渗入心墙内部的流体部分采用无网格法，与有限元相比该方法无需预定义网格信息，克服了对网格的依赖性，节省了划分网格所花费的大量时间和精力。另外，考虑到水力劈裂发生的根本原因是局部存在的高水力梯度，而无网格法又恰好适用于各种具有高梯度、奇异性等特殊问题的数值模拟。因此，采用无网格法模拟心墙水力劈裂中的流体部分具有一定的优势。本发明结合了颗粒离散元和无网格法各自的优点，为堆石坝心墙水力劈裂的数值模拟提供了新的研究思路。

采用颗粒离散元可以较好地模拟心墙骨架的变形特性，采用无网格法计算流体部分克服了对网格的依赖性，节省了划分网格所花费的大量时间和精力。

此外，本发明根据室内试验获取细观参数，模拟实际心墙土料的力学性能。然后，在此基础上，施加心墙的应力边界条件以及堆石坝上游的实际水压力边界条件，进行数值模拟，从细观角度判定心墙是否发生水力劈裂。如发生水力劈裂，则该发明可以实时观测裂纹的扩展过程。最后，根据模拟结果，从施工或设计方面提出改进措施，避免心墙水力劈裂的发生，从而对实际工程具有一定的指导意义。

附图说明

图1为带有初始裂缝的心墙骨架颗粒模型示意图。

图2为流体模型中场节点布置示意图。

图3为数值积分背景网格示意图。

图4为产生贯穿性裂缝的心墙骨架颗粒模型示意图。

图5为本发明一种实施例的流程图。

图6为本发明又一实施例的流程图。

图中：1为墙体，2为心墙骨架圆形颗粒，3为初始裂缝或缺陷部位。

具体实施方式

如图1-6所示，一种堆石坝心墙水力劈裂的数值模拟方法，其步骤如下：

(1)根据室内平面应变试验获取实际心墙土料的应力应变曲线，建立颗粒离散元双轴数值模型；

(1.1)取实际工程的心墙土料作为试样，对其进行室内平面应变试验，绘制应力应变曲线；

(1.2)建立颗粒离散元数值模型，根据室内试验获取细观力学特征参数(如：颗粒大小，刚度，摩擦系数，粘结强度等)使数值模型的宏观应力应变特性与室内试验的结果在整体上具有一致的发展趋势，以保证数值模拟的结果能够对实际工程具有一定的指导意义；

(2)根据颗粒离散元双轴试验确定的细观力学参数，建立心墙骨架的颗粒模型；

(2.1)在模型四周生成墙体，其的作用主要有两方面：一方面四周的墙体将随后生成的颗粒限制在模型内部，另一方面墙体作为刚性加载板对模型施加荷载，在模型内部达到指定的应力状态；

(2.2)心墙堆石坝的防渗土料中通常掺有一定比例的砾石料，以减小心墙和两侧堆石体间的不均匀变形，减小心墙发生水力劈裂的可能性，因此采用颗粒离散元法生成堆石坝心墙模型时，在墙体所围成的计算区域内按一定颗粒级配和孔隙率随机生成骨架颗粒。由于计算机容量的限制，根据相似模拟的原理，将心墙土料的级配转换成数值试验所对应的颗粒级配，保证数值模型与原型料几何相似并具有相同的孔隙率；

(2.3)根据数值试验的级配，获得不同粒径范围的颗粒含量。生成颗粒时，首先生成最大粒径的颗粒组，逐一生成较小的粒组，最后生成最小粒径的粒组；

(2.3.1)用两个二维数组Dmax(i)和Dmin(i)分别存储每一粒组粒径范围的上下限；

(2.3.2)在生成某一粒组的颗粒时，分别指定最大颗粒半径R_max＝Dmax(i)/2.0，最小颗粒半径R_min＝Dmin(i)/2.0，然后在[R_min，R_max]范围内，按均匀分布生成该粒组的所有颗粒，需生成的颗粒数目N由下式确定

$N=\frac{A(1-n)}{\mathrm{\π}{\stackrel{\‾}{R}}^2}$

$\stackrel{\‾}{R}=\frac{R_{\min }+R_{\max }}{2.0}$

式中：A为模型面积，n为模型孔隙率；

(2.3.3)为了在指定空间内成功生成所有颗粒，在生成颗粒时，先将所要生成的颗粒半径均缩小a(a＞1.0)倍，即指定生成颗粒的范围为[R_min/a，R_max/a]，颗粒数为N。所生成的颗粒位置是随机指定的，如果该颗粒与已生成颗粒间的重合量小于允许值，则生成该颗粒，否则，重新随机指定生成位置，直到满足要求；

(2.3.4)当完成所有粒组的颗粒生成之后，再将全体颗粒的半径统一放大a倍，即所有颗粒粒径还原为初始指定粒径；

(2.3.5)将步骤(1.2)中获得的细观参数赋值给颗粒及墙体；

(2.4)循环运行心墙骨架颗粒模型，使颗粒的最大平衡力与颗粒间的最大接触力之比小于指定的允许值，即模型达到平衡状态；

(2.4.1)颗粒粒径还原为初始粒径后，相邻颗粒间必然会有较大的重叠量，根据颗粒离散元的接触刚度模型，可以分别得到颗粒的法向和切向接触力

Fⁿ＝KⁿUⁿ

式中：K_n——接触点处的法向刚度，属于割线模量，与总位移Uⁿ和力Fⁿ相对应；

切向作用力增量与位移增量的关系可写为：

ΔF^s＝-k^sΔU^s

式中：k^s为接触点切向刚度，属于切线模量，与切向位移增量ΔU^s和切向力增量ΔF^s相对应；

(2.4.2)由于颗粒间接触处的重叠，颗粒在各接触点处受到接触力的作用，根据牛顿第二运动定律，颗粒会在不平衡力的作用下发生运动，颗粒离散元程序通过控制最大不平衡力与最大接触力间的比值，最终使模型达到初始平衡状态；

(2.5)通过伺服机制控制模型四周墙体的速度，使模型内部达到指定的初始应力状态；

(2.5.1)根据作用在墙体上的力，计算得到墙体所承受的应力值，并由当前墙体的应力值与目标应力值之差，确定墙体速度(公式如下)，直到墙体应力达到指定的各向同性应力

${\stackrel{\·}{u}}^{\left(w\right)}=G({\mathrm{\σ}}^{\left(w\right)}-{\mathrm{\σ}}^{\left(t\right)})=\mathrm{G\Δ\σ}$

式中：为墙体的运动速度，σ^(w)为墙体应力，σ^(t)为模型内部的目标应力，G为控制墙体速度的参数，其值可由稳定条件确定；

(2.5.2)由于在蓄水之前堆石坝心墙内的最大主应力方向接近铅直方向，因此给模型施加竖向轴压，使其达到指定的σ₁，与此同时，通过调整左右两侧墙体的速度，保持水平向的应力不变；

(2.6)在相互接触的颗粒间设置接触粘结，以模拟粘土心墙，当颗粒间的接触力超过粘结强度，则相应的粘结发生破坏；

(2.7)现有的研究成果表明，心墙发生水力劈裂的必备条件之一是，心墙内部存在初始裂缝或缺陷，而对于完全均质且无任何缺陷的心墙则不会发生水力劈裂现象，因此在生成颗粒流模型时，需在模型内部设置初始裂缝；

(2.7.1)在施工过程中，心墙是成层碾压，层间较其他方向容易形成初始缺陷，因此在模型的左侧中部设置一水平初始裂缝；

(2.7.2)按一定级配生成的颗粒流模型，其内部颗粒呈无序排列，因此需指定某一范围内颗粒间的接触粘结强度为0，以此模拟实际的初始裂缝(见附图3)，并且标记出紧邻裂缝两侧的颗粒；

(3)利用连续介质模型建立流体模型，并对其控制方程进行离散；

(3.1)采用经典土力学中的渗流模型，假设研究的渗流区域内不存在土颗粒骨架，渗透水流被视为连续介质可充满研究域的全部空间，并将土颗粒对水流的作用概化为作用于水流的力；

(3.2)流体计算过程中，采用流体的连续性方程和基于局部平均化的Navier-Stokes(N-S)方程作为渗透水流的控制方程，其分量形式如下

$\frac{\∂n}{\∂t}+\frac{\∂\left({\mathrm{nu}}_i\right)}{{\∂x}_i}=0$

$\frac{\∂\left({\mathrm{nu}}_i\right)}{\∂t}+\frac{\∂\left({\mathrm{nu}}_i{u}_j\right)}{{\∂x}_j}=-\frac{n}{{\mathrm{\ρ}}_f}(\frac{\∂p}{{\∂x}_i}-\frac{f_{\mathrm{inti}}}{n})$

式中：n为土体骨架的孔隙率，u_i，u_j为流体的流速分量，x_i为坐标分量，t为时间，p为流体压力，ρ_f为流体密度，f_int为单位体积内土颗粒与流体之间的相互作用力，其计算公式如下

f_intj≡β_intj(v_j-u_j)

式中：β_intj为流体和固体颗粒间的摩擦系数，其计算公式如下

${\mathrm{\β}}_{\mathrm{intj}}=150\frac{{(1-n)}^2}{{\mathrm{nd}}_p^2}{\mathrm{\μ}}_f+1.75\frac{(1-n){\mathrm{\ρ}}_f}{d_p}|v_j-u_j|,(n\≤0.8)$

${\mathrm{\β}}_{\mathrm{intj}}=\frac{3}{4}{\mathrm{\ρ}}_f{C}_D\frac{{(1-n)n}^{-1.7}}{d_p}|v_j-u_j|,(n>0.8)$

式中：d_p为颗粒的直径，μ_f为流体动力粘度，v_j为固体颗粒的速度分量，u_j流体的速度分量，C_D为拖曳系数，是雷诺数R_e的函数，其计算公式如下

$C_D=\frac{24}{R_e}(1+0.15R_e^{0.687}),(R_e<1000)$

C_D＝0.44         (R_e≥1000)

其中

$R_e=\frac{n|v_j-u_j|d_p}{v_f}$

(3.3)将连续性方程代入局部平均化的N-S方程中得

$\frac{{\mathrm{du}}_i}{\mathrm{dt}}=-\frac{1}{{\mathrm{\&rho;}}_f}(\frac{\&PartialD;p}{{\&PartialD;x}_i}-\frac{f_{\mathrm{inti}}}{n})$

令 $Q_i=-\frac{1}{{\mathrm{\&rho;}}_f}(\frac{\&PartialD;p}{{\&PartialD;x}_i}-\frac{f_{\mathrm{inti}}}{n})$

则 $\frac{{\mathrm{du}}_i}{\mathrm{dt}}=Q_i$

(3.4)对控制方程中的时间进行离散

采用有限差分法对时间求导

$\frac{{\mathrm{du}}_i}{\mathrm{dt}}\left(\mathrm{\&epsiv;}\right)=\frac{u_i^{t+\mathrm{\&Delta;t}}-u_i^t}{\mathrm{\&Delta;t}}$

式中：Δt为时间步长，ε为[t,Δt]内任意时刻；

(4)利用心墙骨架颗粒模型和流体模型模拟水力劈裂

(4.1)将颗粒与流体间的相互耦合力、颗粒之间的作用力之和作用于所有颗粒，利用离散元法计算并得到心墙骨架颗粒模型内部颗粒的速度、应力分布及渗透特性；初始计算时，颗粒与流体间的相互耦合力为零；

(4.2)固体骨架和流体部分的计算交错循环进行，进一步地，在满足计算精度的前提下，为了节省计算成本，经步骤(4.1)循环计算m次后，提取第m次计算得到的渗透特性，利用全局弱形式的无单元EFG法计算并得到渗流流体的速度、压力及颗粒与流体间的相互耦合力；所述m根据计算精度及模拟对象的需要确定，此处取100；具体包括如下步骤：

(4.2.1)选取每一固体颗粒的中心点作为场节点，并利用这组场节点对流体域进行离散；

(4.2.2)通过离散的节点(此处的节点可不同于上述的场节点)在流体域内形成用于数值积分的背景网格，将整体计算域划分为若干个积分子域；

(4.2.3)采用无单元EFG法对离散后的控制方程进行数值积分

$\underset{{\mathrm{\&Omega;}}_j}{\&Integral;}\mathrm{\&Phi;}[\frac{{\hat{u}}_i^{t+\mathrm{\&Delta;t}}-{\hat{u}}_i^t}{\mathrm{\&Delta;t}}-Q_i]d{\mathrm{\&Omega;}}_j=0$

式中：为u_i的移动最小二乘(MLS)近似，Φ是MLS形函数，Ω_j为第j个积分子域；

(4.2.4)通过循环方式，在每一个背景网格内进行数值积分，从而获得整体系统方程组的系数矩阵；

(4.2.5)堆石心墙坝的心墙上游侧为堆石体，其渗透性很大，可认为在模型的左侧边界及初始裂缝内直接作用上游水荷载，因此根据实际的上游水位，给心墙骨架颗粒模型施加已知的压力边界条件；

(4.2.6)求解系统方程组，获得场节点的流速和压力。

(4.3)判断模型内部是否达到平衡状态，如否，回到步骤(3.1)；如是，根据颗粒模型内部是否产生贯穿性裂缝，输出在指定的心墙上游水压力作用下模型内部是否发生水力劈裂的结论，结束模拟过程。

Claims (2)

1.一种堆石坝心墙水力劈裂的数值模拟方法，其特征在于包括如下步骤：

(1)根据室内平面应变试验获取实际心墙土料的应力应变曲线，建立颗粒离散元双轴数值模型；

(2)根据颗粒离散元双轴试验确定的细观力学参数，建立心墙骨架颗粒模型；

(3)利用连续介质模型建立流体模型，并对其控制方程进行离散；

(4)利用心墙骨架颗粒模型和流体模型模拟水力劈裂；

(4.1)将颗粒与流体间的相互耦合力、颗粒之间的作用力之和作用于所有颗粒，利用离散元法计算并得到心墙骨架颗粒模型内部颗粒的速度、应力分布及渗透特性；初始计算时，颗粒与流体间的相互耦合力为零；

(4.2)提取步骤(4.1)计算得到的渗透特性，利用全局弱形式的无单元EFG法计算并得到渗流流体的速度、压力及颗粒与流体间的相互耦合力，具体子步骤如下：

(4.2.1)选取每一固体颗粒的中心点作为场节点，并利用这组场节点对流体域进行离散；

(4.2.2)通过离散的节点在流体域内形成用于数值积分的背景网格，将整体计算域划分为若干个积分子域；

(4.2.3)采用无单元EFG法对离散后的控制方程进行数值积分

$\underset{{\mathrm{\&Omega;}}_j}{\&Integral;}\mathrm{\&Phi;}[\frac{{\hat{u}}_i^{t+\mathrm{\&Delta;t}}-{\hat{u}}_i^t}{\mathrm{\&Delta;t}}-Q_i]d{\mathrm{\&Omega;}}_j=0$

式中：Φ是MLS形函数，u_i为流体的流速分量，为u_i的移动最小二乘近似，Δt为时间步长，Ω_j为第j个积分子域；

(4.2.4)通过循环方式，在每一个背景网格内进行数值积分，从而获得整体系统方程组的系数矩阵；

(4.2.5)根据实际的上游水位，给心墙骨架颗粒模型施加已知的压力边界条件；

(4.2.6)求解系统方程组，获得场节点的流速和压力；

(4.3)判断模型内部是否达到平衡状态，如否，回到步骤(4.1)；如是，根据颗粒模型内部是否产生贯穿性裂缝，输出模型内部在指定的心墙上游水压力作用下是否发生水力劈裂的结论，结束模拟过程。

2.如权利要求1所述的堆石坝心墙水力劈裂的数值模拟方法，其特征在于：所述步骤(4.2)为经步骤(4.1)循环计算m次后，提取第m次计算得到的渗透特性，利用全局弱形式的无单元EFG法计算并得到渗流流体的速度、压力及颗粒与流体间的相互耦合力；所述m根据计算精度及模拟对象的需要确定。