CN113627052B - 一种考虑水力耦合效应的堆石坝流变数值模拟方法 - Google Patents

Abstract

本发明公开了一种考虑水力耦合效应的堆石坝流变数值模拟方法，首先开展考虑渗流作用的堆石坝流变实验，然后采用抛物Signorini型变分不等式方法模拟堆石坝和地基的非稳定渗流过程。为了克服渗流过程中潜在渗出面的溢出点在自由面迭代步中不稳定，采用带自适应罚Heaviside函数的抛物Signorini型变分不等式方法消除渗流溢出点的奇异性，并实现对非稳定渗流溢出点和自由面点的准确定位；通过变形场和渗流场，建立变形和渗流耦合控制方程，建立孔隙率的变化与体积应变关系，最后得到考虑水力耦合效应的堆石坝流变数值模拟结果。本发明揭示渗流作用下面板堆石坝流变机理，为面板堆石坝的建设提供理论依据和技术支撑。

Description

一种考虑水力耦合效应的堆石坝流变数值模拟方法

技术领域

本发明属于堆石坝流变数值模拟技术领域，具体涉及一种考虑水力耦合效应的堆石坝流变数值模拟方法。

背景技术

面板堆石坝防渗体的变形主要取决于堆石体变形，堆石体变形过大，将使面板产生裂缝，从而影响其防渗性能，甚至危及坝体稳定。因此，面板堆石坝工程安全问题归根结底是堆石体的变形问题。面板堆石坝变形包括瞬时变形和时效变形，前者由坝体自重及水荷载作用引起，后者是坝体蓄水后的流变变形。已建面板堆石坝变形监测结果表明，流变变形作为面板堆石坝变形的重要组成部分，是造成大坝长期变形和防渗体失效的主要原因。对于修建在深覆盖层上的高面板堆石坝，由于大规模可压缩地基的存在，并且深覆盖层地基和部分坝体在运行期处于饱水渗流状态，因此，在渗流作用影响下，大坝的流变变形问题更加突出。

目前国内面板堆石坝数值分析中堆石料采用的的本构模型一般可分为非线性弹性模型和弹塑性模型两大类,其中非线性弹性模型主要有DuncanE-v模型、DuncanE-B模型、清华K-G模型,弹塑性模型主要有“南水”双屈服面模型和喻殷宗泽双屈服面模。然而，有关渗流作用对覆盖层砂砾石和堆石料流变特性的影响，已有研究较少涉及。

发明内容

本发明的目的是提供一种考虑水力耦合效应的堆石坝流变数值模拟方法，揭示渗流作用下面板堆石坝流变机理，为面板堆石坝的建设提供理论依据和技术支撑。

本发明所采用的技术方案是，一种考虑水力耦合效应的堆石坝流变数值模拟方法，其特征在于，具体按照以下步骤实施：

步骤1、开展考虑渗流作用的堆石坝流变实验，对比分析不同渗流水压力作用下的流变实验结果，通过数值方法模拟堆石坝流变试验结果，对现有流变模型进行考虑渗流作用的适应性验证，选用指数衰减曲线作为流变模型，根据实验结果修正该流变模型，基于隐式积分迭代方法，获得流动法则方程、硬化方程和一致性方程，建立散粒体、刚性防渗体、以及两者之间界面的本构数值积分统一算法模拟堆石坝坝体和地基非线性变形；

步骤2、采用新近提出的抛物Signorini型变分不等式方法模拟堆石坝和地基的非稳定渗流过程。为了克服渗流过程中潜在渗出面的溢出点在自由面迭代步中不稳定，采用带自适应罚Heaviside函数的抛物Signorini型变分不等式方法消除渗流溢出点的奇异性，并实现对非稳定渗流溢出点和自由面点的准确定位；

步骤3、通过步骤1得到的变形场和步骤2得到的渗流场，建立变形和渗流耦合控制方程，根据Kozeny-Carman模型，在小变形假设下，建立孔隙率的变化与体积应变关系，基于散度定理获得动量和质量守恒方程的等效积分弱形式，采用有限元方法对耦合方程等效积分弱形式进行空间离散，采用双场交叉迭代算法求解渗流—应力耦合过程，直到满足收敛条件；

步骤4、基于面向对象的概念，将步骤1得到的坝体-地基流变模型和堆石料弹塑性模型积分统一算法和步骤3中的渗流—变形耦合分析方法在ABAQUS软件平台上进行数值模拟，得到考虑水力耦合效应的堆石坝流变数值模拟结果。

本发明的特点还在于，

步骤1具体如下：

步骤1.1、采用指数衰减曲线模拟堆石体的流变特性，剪切流变和时间的关系可以表示为：

ε_vf＝a₁+b₁p,B_v＝a₂+b₂p (2)

ε_sf＝a₁+b₁q,B_s＝a₂+b₂q (4)

式中：ε_v为体积应变；ε_s为轴向应变；ε_vf为极限体积流变应变；ε_sf为极限轴向流变应变；σ表示应力状态；t为时间；B_v相当于t＝1是的体积流变量占总量的比值；B_s相当于t＝1是的剪切流变量占总量的比值；a₁、a₂、b₁、b₂四参数均为流变实验拟合值；p＝(σ₁+σ₂+σ₃)/3为八面体正应力；为广义剪应力；σ₁、σ₂、σ₃分别为大中小主应力；

进行有限元计算时，△t时段内粘弹塑性应变增量为：

式中：Δε_v(t,σ)和Δε_s(t,σ)分别为t时刻体积应变增量和剪切应变增量；

步骤1.2、采用南水双屈服面弹塑性模型模拟散粒体材料非线性特性，双屈服面由椭圆函数和幂函数组成，

f₁＝p²+r²q² (7)

式中，r表示屈服面参数，对于石料一般取2，p＝(σ₁+σ₂+σ₃)/3，为八面体正应力；为八面体剪应力；弹塑性应力应变关系的表达式为：

其中：[D]为弹性矩阵；A₁、A₂为塑性系数，A₁、A₂表示如下：

E_t＝E_i(1-R_fR_s)² (12)

式中；r、s表示屈服面参数，一般取2；p＝(σ₁+σ₂+σ₃)/3，为八面体正应力；为八面体剪应力；σ₁、σ₂、σ₃分别为大中小主应力。E_t表示切线弹性模量；E_i表示t时刻的弹性模量；K表示回弹体积变形模量；G表示回弹剪切模量；R_f表示破坏比；R_s表示t时刻的应力水平；E_ur—表示回弹模量；Pa表示标准大气压；μ表示弹性泊松比，可取0.3；v_t表示切线泊松比；C_d表示σ₃等于一个标准大气压时的最大收缩应变；n_d表示收缩体应变随σ₃增加而增加的幂指数；R_d表示发生最大收缩时的(σ₁-σ₃)与偏应力的渐近值(σ₁-σ₃)_ult之比；

对表达式求逆，可得增量形式的弹塑性应力应变关系：

{Δσ}＝[D]_cp{Δε} (17)

式中：[D]_ep为弹塑性矩阵

其中：S_x＝σ_x-P，S_y＝σ_y-P，S_z＝σ_z-P (23)

M₁＝K_p+4G/3，M₂＝K_p-2G/3 (24)

上式：S_x、S_y、S_z为x,y,z方向应力水平；p＝(σ₁+σ₂+σ₃)/3，为八面体正应力；为八面体剪应力；σ₁、σ₂、σ₃分别为大中小主应力。K_p、K_ur、P、M₁、M₂、A、B、C、D、Q为试验参数；G表示回弹剪切模量；s、r为表示屈服面参数，一般取2；A₁、A₂为塑性系数；

弹塑性矩阵中的剪切项系数d_ij＝d_ji不为零，i＝1,2,3，j＝4,5,6，而且d₁₁≠d₂₂≠d₃₃,

结合步骤1.1，得到堆石体内任一点的总应变为

{ε}＝{ε^e}+{ε^p}+{ε_v}+{ε_s} (27)

其增量形式：

{Δε}＝{Δε^e}+{Δε^p}+{Δε_v}+{Δε_s} (28)

得到应力应变增量形式：

{Δσ}＝[D]({Δε^e}+{Δε^p}+{Δε_v}+{Δε_s}) (29)

式中：{ε}为总应变；{ε^e}为瞬时弹性变形；{ε^p}为瞬时塑性变形；{ε_v}为体积流变；{ε_s}为剪切流变。{Δσ}为应力增量；{Δε^e}弹性变形增量；{Δε^p}为塑性变形增量；{Δε_v}体积流变变形增量；{Δε_s}剪切流变变形增量。[D]为弹性矩阵；

通过流动法则方程、硬化方程和一致性方程，得到在复杂应力状态下堆石体流变应力应变增量形式为：

式中：[D_ep]为弹塑性矩阵；[D_n]为弹性矩阵；{△σ^t}为t时刻堆石料流变引起的应力松弛列阵，{△σ_n}为堆石荷载作用n天的应力增量矩阵，{△ε_n}为堆石荷载作用n天的应变增量矩阵；a₁₁、a₁₂、a₂₁、a₂₂是屈服面对应的矩阵系数；{b_n}、{c_n}屈服面f₁、f₂对应的弹性矩阵；q₀是初始剪应力，p₀是初始正应力；为屈服面f₂对应的体积流变；/>公为屈服面f₁对应的剪切流变；式(24)-(30)是为了简化公式(23)存在，μ为泊松比；[D]为弹性矩阵；

步骤1.3、把接触单元类比为薄层单元，薄层单元的本构矩阵为

式中：z为接触面法向；D₁、D₁、D₁、k_zz、k_yz、k_zx均为各向异性弹性本构中弹性矩阵的弹性系数；E、μ分别为弹性模量和泊松比；k_zz0为接触面初始法向刚度；可取D₁；v为法向相对变形；t为接触面厚度；k_zx0和k_yz0为初始切向刚度，可取D₃；c为凝聚力；f为摩擦系数；α为摩擦剪切移角；r为失效比；τ为接触面切应力，γ为接触面切应变；

步骤2具体如下：

步骤2.1、采用自适应惩罚Heaviside函数的变分不等式方法确定饱和渗流过程自由面和移出点，假定潜在出渗面的边界条件为Signorini型互充条件，根据变分不等式的要求，将整个渗流区域定义为一个新边界问题，达西定律重新定义为：

φ＝z+p_ω/ρ_ωg (47)

其中：t为时间；v为渗流速度；k为渗透系数矩阵；φ为总水头；z为垂直位置坐标；P_w为孔隙水压力；H(φ-z)是Heaviside罚函数；ρ_w为水的密度；g为重力加速度；

初始条件和边界条件由下式控制：

式中：u₀为整个区域Ω上的初始位移矢量；为位移边界Γ_u上的预设的位移矢量；为牵引边界上Γ_t上的预设牵引矢量；

步骤2.2、对于渗流过程，Dirichlet和Neumann边界条件为

式中：为水头边界Γ_φ上的预设水头；/>为流量边界Γ_q上的预设流量，对于不透水边界，/>

渗流自由面满足Signorini类型的补充边界条件为

潜在出渗边界Γ_s满足以上边界条件，在出渗点满足φ＝0及q_n＝0；

经有限元离散后，离散型提法的迭代格式表述为：在有限维试探向量空间中求一向量/>使得对/>都有

(ψ-φ^k+1)^TKφ^k+1≥(ψ-φ^k+1)^Tq^k (51)

式中：K为整体渗透劲度矩阵；k_e为单元渗透劲度矩阵；k为自由面的迭代步；n为有限元网格节点总数；B为有限元模型的几何矩阵；H_λ是罚Heaviside函数，自适应罚Heaviside函数的定义如下：

式中：λ₁和λ₂是域每一单元相关联的两个参数，λ₁定义为单元内最低积分点与最低节点的垂直距离；λ₂定义为单元内最高积分点与最高节点的垂直距离；ζ是为了放大λ₁和λ₂两个参数的值而引入的计算参数；

步骤2.3、非稳定渗流过程数值计算过程中，判断渗流自由面准确定位计算的收敛条件为

式中：ε₁和ε₂为指定的容许误差，取ε₁＝ε₂＝10^-3；φ为总水头；i为迭代次数；等于上游面边界上的已知水头值。

通过确定的以上边界条件和收敛条件，得到堆石体和地基的非稳定渗流过程。

步骤3具体如下：

步骤3.1、构建的变形和渗流耦合控制方程如下：

耦合变形和渗流过程受连续介质力学的动量守恒定律支配，所以变形和渗流耦合过程的控制方程描述如下：

式中：D是切向弹性模量张量；u是位移矢量；α是Biot固结系数；δ是Kronecker三角矢量；f是体积力矢量；ε_v是体积应变；S_w是与水压缩性相关的存储量；t为时间；φ＝z+p_w/ρ_wg为总水头；z为垂直位置坐标；p_w为孔隙水压力；H(φ-z)是Heaviside罚函数，φ≧z时，H＝1；φ<z时，H＝0；ρ_w为水的密度；v为渗流速度；

步骤3.2、根据Kozeny-Carman模型，在小变形假设下，孔隙率的变化与体积应变相关，及：

式中：k和k₀分别为当前和初始的渗透系数，n和n₀分别为当前和初始的孔隙率。ε_v为堆石体体积应变；当堆石料颗粒不可变形时β₁取1。

采用双场交叉迭代算法求解渗流—应力耦合过程，直到满足收敛条件。其中收敛准则以前后两次迭代的变形和水头满足指定容许误差作为标准。

|φⁱ⁺¹-φⁱ|≤ε₁,|εⁱ⁺¹-εⁱ|≤ε₂ (59)

式中：φ为总水头；i为迭代次数；ε为应力变形；ε₁和ε₂为指定的容许误差，取ε₁＝ε₂＝10^-3。

步骤4具体如下：

将本项目坝体-地基流变模型和堆石料弹塑性模型积分统一算法和渗流-变形耦合分析方法在ABAQUS软件平台上进行计算，建立考虑渗流作用的深厚覆盖层高面板堆石坝流变计算模型。通过建立的水力耦合模型模拟渗流作用下面板堆石坝的流变变形，得到最终应力变形量。

本发明的有益效果是，一种考虑水力耦合效应的堆石坝流变数值模拟方法，前人研究中，针对渗流作用对材料流变影响机理的研究不够深入，此发明考虑了深覆盖层面板堆石坝的流变受坝体堆石料和地基砂砾石料流变的共同影响，而渗流作用直接影响浸润线以下地基及部分坝体的流变特性。所以基于考虑渗流作用的堆石料和地基砂砾石流变试验、深覆盖层高面板堆石坝流变变形数值仿真，深入揭示渗流作用下深覆盖层高面板堆石坝流变机理。

附图说明

图1是本发明一种考虑水力耦合效应的堆石流变数值模拟方法中点增量法的示意图。

图2是本发明坝体-地基流变模型本构数值积分算法的流程图；

图3是本发明双场交叉迭代过程的流程图；

图4是本发明地基沉降时间曲线数值和实测结果对比图；

图5是本发明地基和坝体孔隙水压力和渗漏速度变化过程图；

图6是本发明考虑水力耦合效应的堆石流变数值模拟方法流程图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明一种考虑水力耦合效应的堆石坝流变数值模拟方法，结合附图1～3，流程图如图6所示，具体按照以下步骤实施：

步骤1、开展考虑渗流作用的堆石坝流变实验，对比分析不同渗流水压力作用下的流变实验结果，通过数值方法模拟堆石坝流变试验结果，对现有流变模型进行考虑渗流作用的适应性验证，选用指数衰减曲线作为流变模型，根据实验结果修正该流变模型，基于隐式积分迭代方法，获得流动法则方程、硬化方程和一致性方程，建立散粒体、刚性防渗体、以及两者之间界面的本构数值积分统一算法模拟堆石坝坝体和地基非线性变形；

步骤1具体如下：

步骤1.1、采用指数衰减曲线模拟堆石体的流变特性，剪切流变和时间的关系可以表示为：

ε_vf＝a₁+b₁p,B_v＝a₂+b₂p (2)

ε_sf＝a₁+b₁q,B_s＝a₂+b₂q (4)

式中：ε_v为体积应变；ε_s为轴向应变；ε_vf为极限体积流变应变；ε_sf为极限轴向流变应变；σ表示应力状态；t为时间；B_v相当于t＝1是的体积流变量占总量的比值；B_s相当于t＝1是的剪切流变量占总量的比值；a₁、a₂、b₁、b₂四参数均为流变实验拟合值；p＝(σ₁+σ₂+σ₃)/3为八面体正应力；为广义剪应力；σ₁、σ₂、σ₃分别为大中小主应力；

因为体积流变和剪切流变具有相同的衰减规律，进行有限元计算时，△t时段内粘弹塑性应变增量为：

式中：Δε_v(t,σ)和Δε_s(t,σ)分别为t时刻体积应变增量和剪切应变增量；

步骤1.2、采用南水双屈服面弹塑性模型模拟散粒体材料非线性特性，双屈服面由椭圆函数和幂函数组成，

f₁＝p²+r²q² (7)

式中，r表示屈服面参数，对于石料一般取2，p＝(σ₁+σ₂+σ₃)/3，为八面体正应力；为八面体剪应力；当采用正交流动时，采用中点增量法求解非线性变形过程，解决堆石体渗透特性和切线变形模量等材料参数对应力或变形的依赖性，弹塑性应力应变关系的表达式为：

其中：[D]为弹性矩阵；A₁、A₂为塑性系数，A₁、A₂表示如下：

E_t＝E_i(1-R_fR_s)² (12)

式中；r、s表示屈服面参数，一般取2；p＝(σ₁+σ₂+σ₃)/3，为八面体正应力；为八面体剪应力；σ₁、σ₂、σ₃分别为大中小主应力。E_t表示切线弹性模量；E_i表示t时刻的弹性模量；K表示回弹体积变形模量；G表示回弹剪切模量；R_f表示破坏比；R_s表示t时刻的应力水平；E_ur—表示回弹模量；Pa表示标准大气压；μ表示弹性泊松比，可取0.3；v_t表示切线泊松比；C_d表示σ₃等于一个标准大气压时的最大收缩应变；n_d表示收缩体应变随σ₃增加而增加的幂指数；R_d表示发生最大收缩时的(σ₁-σ₃)与偏应力的渐近值(σ₁-σ₃)_ult之比；

对表达式求逆，可得增量形式的弹塑性应力应变关系：

{Δσ}＝[D]_ep{Δε} (17)

式中：[D]_ep为弹塑性矩阵

/>

其中：S_x＝σ_x-P，S_y＝σ_y-P，S_z＝σ_z-P (23)

M₁＝K_p+4G/3，M₂＝K_p-2G/3 (24)

上式：S_x、S_y、S_z为x,y,z方向应力水平；p＝(σ₁+σ₂+σ₃)/3，为八面体正应力；为八面体剪应力；σ₁、σ₂、σ₃分别为大中小主应力。K_p、K_ur、P、M₁、M₂、A、B、C、D、Q为试验参数；G表示回弹剪切模量；s、r为表示屈服面参数，一般取2；A₁、A₂为塑性系数；

弹塑性矩阵中的剪切项系数d_ij＝d_ji不为零，i＝1,2,3，j＝4,5,6，因而可以考虑堆石料的剪胀特性，而且d₁₁≠d₂₂≠d₃₃,从而可以考虑应力引起的各向异性；

结合步骤1.1，得到堆石体内任一点的总应变为

{ε}＝{ε^e}+{ε^p}+{ε_v}+{ε_s} (27)其增量形式：

{Δε}＝{Δε^e}+{Δε^p}+{Δε_v}+{Δε_s} (28)得到应力应变增量形式：

{Δσ}＝[D]({Δε^e}+{Δε^p}+{Δε_v}+{Δε_s}) (29)

式中：{ε}为总应变；{ε^e}为瞬时弹性变形；{ε^p}为瞬时塑性变形；{ε_v}为体积流变；{ε_s}为剪切流变。{Δσ}为应力增量；{Δε^e}弹性变形增量；{Δε^p}为塑性变形增量；{Δε_v}体积流变变形增量；{Δε_s}剪切流变变形增量。[D]为弹性矩阵；

通过流动法则方程、硬化方程和一致性方程，得到在复杂应力状态下堆石体流变应力应变增量形式为：

/>

式中：[D_ep]为弹塑性矩阵；[D_n]为弹性矩阵；{△σ^t}为t时刻堆石料流变引起的应力松弛列阵，{△σ_n}为堆石荷载作用n天的应力增量矩阵，{△ε_n}为堆石荷载作用n天的应变增量矩阵；a₁₁、a₁₂、a₂₁、a₂₂是屈服面对应的矩阵系数；{b_n}、{c_n}屈服面f₁、f₂对应的弹性矩阵；q₀是初始剪应力，p₀是初始正应力；为屈服面f₂对应的体积流变；/>公为屈服面f₁对应的剪切流变；式(24)-(30)是为了简化公式(23)存在，μ为泊松比；[D]为弹性矩阵；

步骤1.3、把接触单元类比为薄层单元，从而有效地改善了接触面附近应力分布的连续性。薄层单元的本构矩阵为

式中：z为接触面法向；D₁、D₁、D₁、k_zz、k_yz、k_zx均为各向异性弹性本构中弹性矩阵的弹性系数；E、μ分别为弹性模量和泊松比；k_zz0为接触面初始法向刚度；可取D₁；v为法向相对变形；t为接触面厚度；k_zx0和k_yz0为初始切向刚度，可取D₃；c为凝聚力；f为摩擦系数；α为摩擦剪切移角；r为失效比；τ为接触面切应力，γ为接触面切应变；

步骤2、采用新近提出的抛物Signorini型变分不等式方法模拟堆石坝和地基的非稳定渗流过程。为了克服渗流过程中潜在渗出面的溢出点在自由面迭代步中不稳定，采用带自适应罚Heaviside函数的抛物Signorini型变分不等式方法消除渗流溢出点的奇异性，并实现对非稳定渗流溢出点和自由面点的准确定位；

步骤2具体如下：

步骤2.1、采用自适应惩罚Heaviside函数的变分不等式方法确定饱和渗流过程自由面和移出点，假定潜在出渗面的边界条件为Signorini型互充条件，根据变分不等式的要求，将整个渗流区域定义为一个新边界问题，达西定律重新定义为：

φ＝z+p_ω/ρ_ωg (47)

其中：t为时间；v为渗流速度；k为渗透系数矩阵；φ为总水头；z为垂直位置坐标；P_w为孔隙水压力；H(φ-z)是Heaviside罚函数；ρ_w为水的密度；g为重力加速度；

初始条件和边界条件由下式控制：

式中：u₀为整个区域Ω上的初始位移矢量；为位移边界Γ_u上的预设的位移矢量；为牵引边界上Γ_t上的预设牵引矢量；

步骤2.2、对于渗流过程，Dirichlet和Neumann边界条件为

式中：为水头边界Γ_φ上的预设水头；/>为流量边界Γ_q上的预设流量，对于不透水边界，/>

渗流自由面满足Signorini类型的补充边界条件为

潜在出渗边界Γ_s满足以上边界条件，在出渗点满足φ＝0及q_n＝0；

经有限元离散后，离散型提法的迭代格式表述为：在有限维试探向量空间中求一向量/>使得对/>都有

(ψ-φ^k+1)^TKφ^k+1≥(ψ-φ^k+1)^Tq^k (51)

/>

式中：K为整体渗透劲度矩阵；k_e为单元渗透劲度矩阵；k为自由面的迭代步；n为有限元网格节点总数；B为有限元模型的几何矩阵；H_λ是为了消除自由面迭代过程中可能出现的数值不稳定性及网格依赖性而引入的罚Heaviside函数，在数值计算过程中，根据问题的收敛条件，逐步放宽罚Heaviside函数的定义，直至获得问题的解答。自适应罚Heaviside函数的定义如下：

式中：λ₁和λ₂是域每一单元相关联的两个参数，λ₁定义为单元内最低积分点与最低节点的垂直距离；λ₂定义为单元内最高积分点与最高节点的垂直距离；ζ是为了放大λ₁和λ₂两个参数的值而引入的计算参数，ζ建议在1～10之间取值。程序开始时，ζ值取1；计算过程中，当程序经过一定步数的迭代难以收敛时，根据问题的收敛条件和单元网格尺寸，以0.5～1的增量逐步增大ζ值；

步骤2.3、非稳定渗流过程数值计算过程中，判断渗流自由面准确定位计算的收敛条件为

式中：ε₁和ε₂为指定的容许误差，取ε₁＝ε₂＝10^-3；φ为总水头；i为迭代次数；等于上游面边界上的已知水头值。

通过确定的以上边界条件和收敛条件，得到堆石体和地基的非稳定渗流过程。

步骤3、基于连续介质力学和动量守恒定律，在Biot理论框架下，通过步骤1得到的变形场和步骤2得到的渗流场，建立变形和渗流耦合控制方程，根据Kozeny-Carman模型，在小变形假设下，建立孔隙率的变化与体积应变关系，基于散度定理获得动量和质量守恒方程的等效积分弱形式，采用有限元方法对耦合方程等效积分弱形式进行空间离散，采用双场交叉迭代算法求解渗流—应力耦合过程，直到满足收敛条件；

步骤3具体如下：

步骤3.1、构建的变形和渗流耦合控制方程如下：

耦合变形和渗流过程受连续介质力学的动量守恒定律支配，所以基于连续介质力学和动量守恒定律，在Biot理论框架下，变形和渗流耦合过程的控制方程描述如下：

式中：D是切向弹性模量张量；u是位移矢量；α是Biot固结系数；δ是Kronecker三角矢量；f是体积力矢量；ε_v是体积应变；S_w是与水压缩性相关的存储量；t为时间；φ＝z+p_w/ρ_wg为总水头；z为垂直位置坐标；p_w为孔隙水压力；H(φ-z)是Heaviside罚函数，φ≧z时，H＝1；φ<z时，H＝0；ρ_w为水的密度；v为渗流速度；

步骤3.2、堆石体的渗透特性与应力或变形的依赖关系是研究渗流-应力耦合过程的关键。堆石和地基的变形和孔隙率的改变将造成渗透性的变化，为了反应堆石体的渗透特性与应力或变形的依赖关系，根据Kozeny-Carman模型，在小变形假设下，孔隙率的变化与体积应变相关，及：

式中：k和k₀分别为当前和初始的渗透系数，n和n₀分别为当前和初始的孔隙率。ε_v为堆石体体积应变；当堆石料颗粒不可变形时β₁取1。

采用双场交叉迭代算法求解渗流—应力耦合过程，直到满足收敛条件。

其中收敛准则以前后两次迭代的变形和水头满足指定容许误差作为标准。

式中：φ为总水头；i为迭代次数；ε为应力变形；ε₁和ε₂为指定的容许误差，取ε₁＝ε₂＝10^-3。

步骤4、基于面向对象的概念，将步骤1得到的坝体-地基流变模型和堆石料弹塑性模型积分统一算法和步骤3中的渗流—变形耦合分析方法在ABAQUS软件平台上进行数值模拟，得到考虑水力耦合效应的堆石坝流变数值模拟结果。

步骤4具体如下：

将本项目坝体-地基流变模型和堆石料弹塑性模型积分统一算法和渗流-变形耦合分析方法在ABAQUS软件平台上进行计算，建立考虑渗流作用的深厚覆盖层高面板堆石坝流变计算模型。通过建立的水力耦合模型模拟渗流作用下面板堆石坝的流变变形，得到最终应力变形量。

为了验证本发明面板堆石坝水力耦合模型的合理性和准确性，以苗家坝和九甸峡面板堆石坝工程为研究实例，对其开展数值仿真计算。将应力变形计算成果与实际监测资料开展对比分析，对比结果见图4、图5。

从图4、图5中可以看出，采用本发明进行计算的结果与实际工程结果基本保持一致，在施工期的结果与实际工程较接近。蓄水后，随着水荷载的增大，模拟计算结果有所偏差。综上，可以看出考虑水力耦合作用的面板堆石坝变形比不考虑水力耦合作用的变形更符合实际工程。蓄水过程中，实测和数值计算水压力均增加，孔隙水压力在蓄水完成时达到峰值，之后缓慢减小并趋于稳定。采用本发明进行数值计算获得的孔隙水压力与实测结果吻合较好。渗漏速度变化过程与孔隙水压力变化过程基本相似。

本发明考虑了渗流对面板堆石坝流变的影响，面板堆石坝的流变受坝体堆石料流变的影响，而渗流作用直接影响浸润线以下地基及部分坝体的流变特性。基于考虑渗流作用的堆石料流变试验、面板堆石坝流变变形数值仿真及数据挖掘结果，揭示渗流作用下面板堆石坝流变机理。为面板堆石坝建设提供理论依据和技术支撑。

Claims (4)

1.一种考虑水力耦合效应的堆石坝流变数值模拟方法，其特征在于，具体按照以下步骤实施：

步骤1、开展考虑渗流作用的堆石坝流变实验，对比分析不同渗流水压力作用下的流变实验结果，通过数值方法模拟堆石坝流变试验结果，对现有流变模型进行考虑渗流作用的适应性验证，选用指数衰减曲线作为流变模型，根据实验结果修正该流变模型，基于隐式积分迭代方法，获得流动法则方程、硬化方程和一致性方程，建立散粒体、刚性防渗体、以及两者之间界面的本构数值积分统一算法模拟堆石坝坝体和地基非线性变形；

所述步骤1具体如下：

步骤1.1、采用指数衰减曲线模拟堆石体的流变特性，剪切流变和时间的关系可以表示为：

ε_vf＝a₁+b₁p,B_v＝a₂+b₂p (2)

ε_sf＝a₁+b₁q,B_s＝a₂+b₂q (4)

式中：ε_v为体积应变；ε_s为轴向应变；ε_vf为极限体积流变应变；ε_sf为极限轴向流变应变；σ表示应力状态；t为时间；B_v相当于t＝1是的体积流变量占总量的比值；B_s相当于t＝1是的剪切流变量占总量的比值；a₁、a₂、b₁、b₂四参数均为流变实验拟合值；p＝(σ₁+σ₂+σ₃)/3为八面体正应力；为广义剪应力；σ₁、σ₂、σ₃分别为大中小主应力；

进行有限元计算时，△t时段内粘弹塑性应变增量为：

式中：Δε_v(t,σ)和Δε_s(t,σ)分别为t时刻体积应变增量和剪切应变增量；

步骤1.2、采用南水双屈服面弹塑性模型模拟散粒体材料非线性特性，双屈服面由椭圆函数和幂函数组成，

f₁＝p²+r²q² (7)

式中，r表示屈服面参数，对于石料取2，p＝(σ₁+σ₂+σ₃)/3，为八面体正应力；为八面体剪应力；弹塑性应力应变关系的表达式为：

其中：[D]为弹性矩阵；A₁、A₂为塑性系数，A₁、A₂表示如下：

E_t＝E_i(1-R_fR_s)² (12)

式中；r、s表示屈服面参数，值取2；p＝(σ₁+σ₂+σ₃)/3，为八面体正应力；为八面体剪应力；σ₁、σ₂、σ₃分别为大中小主应力，E_t表示切线弹性模量；E_i表示t时刻的弹性模量；K表示回弹体积变形模量；G表示回弹剪切模量；R_f表示破坏比；R_s表示t时刻的应力水平；E_ur—表示回弹模量；Pa表示标准大气压；μ表示弹性泊松比，可取0.3；v_t表示切线泊松比；C_d表示σ₃等于一个标准大气压时的最大收缩应变；n_d表示收缩体应变随σ₃增加而增加的幂指数；R_d表示发生最大收缩时的(σ₁-σ₃)与偏应力的渐近值(σ₁-σ₃)_ult之比；

对表达式求逆，可得增量形式的弹塑性应力应变关系：

{Δσ}＝[D]_ep{Δε} (17)

式中：[D]_ep为弹塑性矩阵

其中：S_x＝σ_x-P，S_y＝σ_y-P，S_z＝σ_z-P (23)

M₁＝K_p+4G/3，M₂＝K_p-2G/3 (24)

上式：S_x、S_y、S_z为x,y,z方向应力水平；p＝(σ₁+σ₂+σ₃)/3，为八面体正应力；为八面体剪应力；σ₁、σ₂、σ₃分别为大中小主应力，K_p、K_ur、P、M₁、M₂、A、B、C、D、Q为试验参数；G表示回弹剪切模量；s、r为表示屈服面参数；A₁、A₂为塑性系数；

弹塑性矩阵中的剪切项系数d_ij＝d_ji不为零，i＝1,2,3，j＝4,5,6，而且d₁₁≠d₂₂≠d₃₃,

结合步骤1.1，得到堆石体内任一点的总应变为{ε}＝{ε^e}+{ε^p}+{ε_v}+{ε_s} (27)

其增量形式：

{Δε}＝{Δε^e}+{Δε^p}+{Δε_v}+{Δε_s} (28)

得到应力应变增量形式：

{Δσ}＝[D]({Δε^e}+{Δε^p}+{Δε_v}+{Δε_s}) (29)式中：{ε}为总应变；{ε^e}为瞬时弹性变形；{ε^p}为瞬时塑性变形；{ε_v}为体积流变；{ε_s}为剪切流变，{Δσ}为应力增量；{Δε^e}弹性变形增量；{Δε^p}为塑性变形增量；{Δε_v}体积流变变形增量；{Δε_s}剪切流变变形增量，[D]为弹性矩阵；

通过流动法则方程、硬化方程和一致性方程，得到在复杂应力状态下堆石体流变应力应变增量形式为：

式中：[D_ep]为弹塑性矩阵；[D_n]为弹性矩阵；{△σ^t}为t时刻堆石料流变引起的应力松弛列阵，{△σ_n}为堆石荷载作用n天的应力增量矩阵，{△ε_n}为堆石荷载作用n天的应变增量矩阵；a₁₁、a₁₂、a₂₁、a₂₂是屈服面对应的矩阵系数；{b_n}、{c_n}屈服面f₁、f₂对应的弹性矩阵；q₀是初始剪应力，p₀是初始正应力；为屈服面f₂对应的体积流变；/>公为屈服面f₁对应的剪切流变；式(24)-(30)是为了简化公式(23)存在，μ为泊松比；[D]为弹性矩阵；

步骤1.3、把接触单元类比为薄层单元，薄层单元的本构矩阵为

式中：z为接触面法向；D₁、D₁、D₁、k_zz、k_yz、k_zx均为各向异性弹性本构中弹性矩阵的弹性系数；E、μ分别为弹性模量和泊松比；k_zz0为接触面初始法向刚度；可取D₁；v为法向相对变形；t为接触面厚度；k_zx0和k_yz0为初始切向刚度，可取D₃；c为凝聚力；f为摩擦系数；α为摩擦剪切移角；r为失效比；τ为接触面切应力，γ为接触面切应变；

步骤2、采用新近提出的抛物Signorini型变分不等式方法模拟堆石坝和地基的非稳定渗流过程，采用带自适应罚Heaviside函数的抛物Signorini型变分不等式方法消除渗流溢出点的奇异性，并实现对非稳定渗流溢出点和自由面点的准确定位；

步骤3、通过步骤1得到的变形场和步骤2得到的渗流场，建立变形和渗流耦合控制方程，根据Kozeny-Carman模型，在小变形假设下，建立孔隙率的变化与体积应变关系，基于散度定理获得动量和质量守恒方程的等效积分弱形式，采用有限元方法对耦合方程等效积分弱形式进行空间离散，采用双场交叉迭代算法求解渗流—应力耦合过程，直到满足收敛条件；

步骤4、基于面向对象的概念，将步骤1得到的坝体-地基流变模型和堆石料弹塑性模型积分统一算法和步骤3中的渗流—变形耦合分析方法在ABAQUS软件平台上进行数值模拟，得到考虑水力耦合效应的堆石坝流变数值模拟结果。

2.根据权利要求1所述的一种考虑水力耦合效应的堆石坝流变数值模拟方法，其特征在于，所述步骤2具体如下：

步骤2.1、采用自适应惩罚Heaviside函数的变分不等式方法确定饱和渗流过程自由面和移出点，假定潜在出渗面的边界条件为Signorini型互充条件，根据变分不等式的要求，将整个渗流区域定义为一个新边界问题，达西定律重新定义为：

φ＝z+p_ω/ρ_ωg (47)

其中：t为时间；v为渗流速度；k为渗透系数矩阵；φ为总水头；z为垂直位置坐标；P_w为孔隙水压力；H(φ-z)是Heaviside罚函数；ρ_w为水的密度；g为重力加速度；

初始条件和边界条件由下式控制：

式中：u₀为整个区域Ω上的初始位移矢量；为位移边界Γ_u上的预设的位移矢量；/>为牵引边界上Γ_t上的预设牵引矢量；

步骤2.2、对于渗流过程，Dirichlet和Neumann边界条件为

式中：为水头边界Γ_φ上的预设水头；/>为流量边界Γ_q上的预设流量，对于不透水边界，/>

渗流自由面满足Signorini类型的补充边界条件为

潜在出渗边界Γ_s满足以上边界条件，在出渗点满足φ＝0及q_n＝0；

经有限元离散后，离散型提法的迭代格式表述为：在有限维试探向量空间中求一向量/>使得对/>都有

(ψ-φ^k+1)^TKφ^k+1≥(ψ-φ^k+1)^Tq^k (51)

式中：K为整体渗透劲度矩阵；k_e为单元渗透劲度矩阵；k为自由面的迭代步；n为有限元网格节点总数；B为有限元模型的几何矩阵；H_λ是罚Heaviside函数，自适应罚Heaviside函数的定义如下：

式中：λ₁和λ₂是域每一单元相关联的两个参数，λ₁定义为单元内最低积分点与最低节点的垂直距离；λ₂定义为单元内最高积分点与最高节点的垂直距离；ζ是为了放大λ₁和λ₂两个参数的值而引入的计算参数；

步骤2.3、非稳定渗流过程数值计算过程中，判断渗流自由面准确定位计算的收敛条件为

||φⁱ⁺¹-φⁱ||₁＜ε₁||φⁱ||₁且||φⁱ⁺¹-φⁱ||_∞＜ε₂||φⁱ||_∞ (56)

式中：ε₁和ε₂为指定的容许误差，取ε₁＝ε₂＝10^-3；φ为总水头；i为迭代次数；||φⁱ||_∞等于上游面边界上的已知水头值，通过确定的以上边界条件和收敛条件，得到堆石体和地基的非稳定渗流过程。

3.根据权利要求2所述的一种考虑水力耦合效应的堆石坝流变数值模拟方法，其特征在于，所述步骤3具体如下：

步骤3.1、构建的变形和渗流耦合控制方程如下：

耦合变形和渗流过程受连续介质力学的动量守恒定律支配，所以变形和渗流耦合过程的控制方程描述如下：

式中：D是切向弹性模量张量；u是位移矢量；α是Biot固结系数；δ是Kronecker三角矢量；f是体积力矢量；ε_v是体积应变；S_w是与水压缩性相关的存储量；t为时间；φ＝z+p_w/ρ_wg为总水头；z为垂直位置坐标；p_w为孔隙水压力；H(φ-z)是Heaviside罚函数，φ≧z时，H＝1；φ<z时，H＝0；ρ_w为水的密度；v为渗流速度；

步骤3.2、根据Kozeny-Carman模型，在小变形假设下，孔隙率的变化与体积应变相关，及：

式中：k和k₀分别为当前和初始的渗透系数，n和n₀分别为当前和初始的孔隙率，ε_v为堆石体体积应变；当堆石料颗粒不可变形时β₁取1，

采用双场交叉迭代算法求解渗流—应力耦合过程，直到满足收敛条件，其中收敛准则以前后两次迭代的变形和水头满足指定容许误差作为标准，

|φⁱ⁺¹-φⁱ|≤ε₁,|εⁱ⁺¹-εⁱ|≤ε₂ (59)

式中：φ为总水头；i为迭代次数；ε为应力变形；ε₁和ε₂为指定的容许误差，取ε₁＝ε₂＝10^-3。

4.根据权利要求3所述的一种考虑水力耦合效应的堆石坝流变数值模拟方法，其特征在于，所述步骤4具体如下：

将坝体-地基流变模型和堆石料弹塑性模型积分统一算法和渗流-变形耦合分析方法在ABAQUS软件平台上进行计算，建立考虑渗流作用的深厚覆盖层高面板堆石坝流变计算模型，通过建立的水力耦合模型模拟渗流作用下面板堆石坝的流变变形，得到最终应力变形量。