CN102621009A - 模拟堆石体长期变形的试验方法 - Google Patents

Abstract

一种模拟堆石体长期变形的试验方法，其步骤包括：⑴根据室内三轴试验获取的堆石体等时应力应变曲线，建立强度劣化模型；⑵三维堆石颗粒集合的随机模拟；⑶利用堆石颗粒随机模型和强度劣化模型计算并绘制堆石颗粒流变的长期变形曲线。其优点是：本发明的模拟试验方法可以实时观察颗粒在模拟流变过程中的滑移、旋转和破碎，方便地提取各个组构量，通过试验数据绘制堆石体的长期变形曲线，可更加准确地求取非线性应力应变参数，为堆石坝的数值计算提供依据，并能对一些堆石坝工程的长期变形问题等的研究提供一定的参考依据和理论指导。

Description

模拟堆石体长期变形的试验方法

技术领域

本发明涉及水利水电工程和岩土工程技术领域，具体的说是一种模拟堆石体长期变形的试验方法。

背景技术

已建堆石坝的原型观测资料表明，在蓄水后大坝变形并未结束，在相当长的一段时间内还会继续增加。长期变形对堆石坝的安全运行有负面影响，例如过大的长期变形会导致混凝土面板破坏，从而影响其防渗性能，甚至危及坝体的安全。研究堆石体长期变形对于掌握混凝土面板堆石坝体的应力与变形规律、实行大坝与止水系统的运行监控和安全评估、科学指导坝工建设、保障工程质量安全等方面有着重要的指导意义。

在面板堆石坝的应力与变形分析中，选用合理的堆石体本构模型以及准确的模型参数是整个分析的关键，堆石体的参数一般由室内或现场试验获得，然而受试验条件、缩尺效应的限制和堆石材料自身性质的离散性，使得测定的力学特性参数与实际值存在一定的差异，由此计算的堆石坝应力与变形与实测值差别较大，如室内流变试验一般几个月就完成了，而坝体的后期变形往往持续数年。

在堆石体流变的室内试验中，由于试样尺寸及其最大粒径的限制，高接触应力-破碎和颗粒重新排列-应力释放、调整和转移的循环过程很快结束，并进入单纯的颗粒重新排列过程，同样由于试样体积有限，这种单纯颗粒重新排列过程引起的变形效应难以测出，所以堆石体的室内流变试验宏观上表现为变形迅速平稳。且受试验仪器的限制，室内试验难以对堆石细观组构进行动态观测。而新兴的模拟试验方法可以实时观察颗粒在加载过程中的滑移、旋转和破碎，方便地提取各个组构量，是研究堆石体长期变形的有效手段。

发明内容

本发明的目的是提供一种模拟堆石体长期变形的试验方法，该方法可模拟堆石料的长期变形这一复杂的、非线性演化问题，可以实时观察颗粒在加载过程中的滑移、旋转和破碎，方便地提取各个组构量，是研究堆石料长期变形的有效手段。

本发明模拟堆石体长期变形的试验方法，其步骤如下：

(1)根据室内三轴试验获取的堆石体等时应力应变曲线，建立强度劣化模型

(1.1)取样，对堆石体试样进行三轴试验，绘制等时应力应变曲线；

(1.1.1)试样安装；

(1.1.2)启动电动机，控制剪切应变速率，开始剪切；试样每产生一定的轴向应变，测计一次测力计读数和轴向变形值；

(1.1.3)当测力计读数出现峰值时，剪切继续进行，至超过5％的轴向应变为止；当测力计读数无峰值时，剪切进行到轴向应变为15％～20％；

(1.1.4)关闭电动机，拆除试样，整理数据，绘制等时应力应变曲线；

(1.2)根据等时应力应变曲线建立时间变化与堆石颗粒的强度、劣化程度、劣化速度之间的等量关系，得到以下公式：

$f\left(\mathrm{\τ}\right)=f_0(1-\frac{\mathrm{a\τ}}{\mathrm{\τ}+n})$

式中：f(τ)表示τ时刻的强度；f₀为瞬时强度；a为反映劣化程度的参数；n为反映强度劣化快慢的参数；

(2)三维堆石颗粒集合的随机模拟

(2.1)建立堆石颗粒随机模型

(2.1.1)按照颗粒级配曲线生成不等边椭球的半径

a＝r₁+(r₂-r₁)rand1

b＝r₁+(r₂-r₁)rand2

c＝r₁+(r₂-r₁)rand3

式中：a，b，c不等边椭球的三个极半轴的长度，r₁为某一粒径组的粒径下限，r₂为相应粒径组的粒径上限，rand1、rand2、rand3为区间[0，1]均匀分布上的独立随机数；

(2.1.2)在[n_min，n_max]之间均匀分布颗粒的顶点数：

n＝n_min+(n_max-n_min)rand

式中：rand为区间[0，1]均匀分布上的随机数，n_min、n_max根据要模拟的颗粒的实际形状随机取值；若实际颗粒棱角较多，则n_min、n_max都相应取大些，若棱角较少则n_min、n_max取小些；若堆石体中不同颗粒的棱角数目差不多，则n_min、n_max的值相差较小，若不同颗粒的棱角数目相差较大，则n_min、n_max的值相差较大；

(2.1.3)采用球坐标确定凸多面体的顶点：

z_i＝z₀+ccosθ_i

式中：x_i，y_i，z_i不等边椭球上第i个点的坐标；x₀，y₀，z₀是不等边椭球的球心坐标；θ_i是球坐标系中的天顶，在[0，π]区间内均匀分布；

是球坐标系中的方位角，在[0，2π]区间内均匀分布；

(2.1.4)在不等边椭球上随机布完点后，遍历每个顶点p_i，寻找与顶点p_i距离最近的顶点p_i1，然后再在剩下的顶点中找一个点p_i2，使得其余所有顶点均在p_i、p_i1、p_i2这三个点确定的三角形平面的同一侧；遍历完所有顶点后，删除具有相同节点的面；

(2.1.5)采用判别准则判别相互侵入关系

根据判别准则，如果两个凸多面体是相互分离的，进入步骤(2.1.6)；如果两个凸多面体相互侵入，回到步骤(2.1.1)重新生成凸多面体；

本发明三维颗粒侵入判别准则如下：当凸多面体A₁A₂A₃…A_n所有顶点在凸多面体B₁B₂B₃…B_n任意一个面B₁B₂B₃的一侧，而凸多面体B₁B₂B₃…B_n的形心在另一侧，那么可以断定这2个凸多面体是相互分离的；此种方法简单明了，避免了前述研究中判别方法中出现的不必要的麻烦；

由于所有颗粒均在不等边椭球面上内接而成，为了进一步加快颗粒侵入的判断，当两颗粒球心的距离大于两颗粒最大极半轴之和时，颗粒是互不侵入的；

(2.1.6)存储经上述步骤生成的全部凸多面体的几何信息；

(2.2)插入界面单元发生颗粒破碎

在颗粒内部的所有细观单元之间插入界面单元；当界面单元的刚度降低到0时，界面单元完全失效，显示颗粒发生破碎；

(3)利用堆石颗粒随机模型和强度劣化模型计算并绘制堆石颗粒流变的长期变形曲线

(3.1)将整个流变过程平均划分为N个时间步长；

(3.2)根据强度劣化模型计算堆石颗粒从0时刻开始每经过一个时间步长后新的强度；

(3.3)根据所得强度对堆石颗粒随机模型进行平衡迭代计算，求取并保存不同强度对应的静力平衡状态的变形值；

(3.4)将得到的变形值用平滑的曲线串联，得到整个流变过程的长期变形曲线。

所述整个模拟过程所需要的时间，为本领域技术人员对于不同岩性岩样的预估流变所需时间，可以输入后根据结果进行调整，不影响模拟结果；一般地，输入的时间不小于长期变形曲线趋平的临界点时间。

本发明模拟堆石体长期变形的试验方法的优点是：本发明的模拟试验方法可以实时观察颗粒在模拟流变过程中的滑移、旋转和破碎，方便地提取各个组构量，经与室内试验对比，符合试验结果，通过模拟试验数据绘制堆石体的长期变形曲线，可更加准确地求取非线性应力应变参数，如长科院幂函数流变模型参数c、d、η、m、c_α、d_α、c_β、d_β、λ_v等，为堆石坝的数值计算提供依据，并能对一些堆石坝工程的长期变形问题等的研究提供一定的参考依据和理论指导。

附图说明

图1为本发明与室内试验级配曲线对比图。

图2为堆石颗粒随机模型示意图。

图3为界面单元失效颗粒破碎示意图。

图4为强度劣化模型原理示意图。

图5为堆石颗粒强度随时间步长劣化流程图。

图6为本发明与室内试验的岩石强度与时间关系曲线对比图。

图7为本发明与室内试验的轴向应变曲线对比图。

图8为本发明不同颗粒强度下的轴向应变曲线图。

图9为本发明与室内试验的体积应变曲线对比图。

图10为本发明不同颗粒强度下的体积应变曲线图。

具体实施方式

下面结合附图，对本发明进行进一步说明：根据图1-10所示，一种模拟堆石体长期变形的试验方法，其步骤如下：

(1)根据室内三轴试验获取的堆石体等时应力应变曲线，建立强度劣化模型

(1.1)取样，对堆石体试样进行三轴试验，绘制等时应力应变曲线；

(1.1.1)试样安装；

(1.1.2)启动电动机，控制剪切应变速率，开始剪切；试样每产生一定的轴向应变，测计一次测力计读数和轴向变形值；

(1.1.3)当测力计读数出现峰值时，剪切继续进行，至超过5％的轴向应变为止；当测力计读数无峰值时，剪切进行到轴向应变为15％～20％；

(1.1.4)关闭电动机，拆除试样，整理数据，绘制等时应力应变曲线；

(1.2)根据步骤(1.1)不同岩性试样的流变试验得到等时应力应变曲线，确定各时刻的强度值，采用双曲线型式进行拟合，本实施例中的试验数据来自文献：许洪发，软岩强度和弹模的时间效应研究；李良权等，基于流变试验的向家坝砂岩长期强度评价；沈明荣等，红砂岩长期强度特性的试验研究；得到公式：

$f\left(\mathrm{\τ}\right)=f_0(1-\frac{\mathrm{a\τ}}{\mathrm{\τ}+n})$

式中：f(τ)表示τ时刻的强度；f₀为瞬时强度；a＝(f₀-f_∞)/f₀反映了劣化程度，f_∞为长期强度，a越大，劣化越显著；n为反映强度劣化快慢的参数，n越小，劣化越快。根据不同岩性试样的蠕变试验确定各时刻的强度值，采用此强度劣化模型进行拟合，曲线拟合的很好，相关系数均在0.99以上，参见图6。

模型的强度参数有f_t、c、

目前对岩石抗拉长期强度研究较少，在此假设岩石在任意时刻的压拉强度比保持不变，即f_t(τ)＝f_c(τ)/CT，下标c和t分别表示抗压强度和抗拉强度，CT为压拉强度比。大量的试验资料表明，岩土介质的长期抗剪强度低于瞬时强度值，这里主要指粘聚力c和内摩擦角

且相比于粘聚力的降低程度，内摩擦角的变化要小得多，为此本发明假定内摩擦角不随时间变化，强度的劣化体现在粘聚力和抗拉强度随时间减小。考虑强度劣化的界面单元破坏准则可表示为图4。

(2)三维堆石体颗粒集合的随机模拟

对于人工堆石材料，其颗粒的形状和空间分布具有显著的随机性，可通过抽样统计的方法获取堆石体颗粒的特征参数，如级配、颗粒尺寸、颗粒形状等。本发明运用随机模拟技术生成形状随机的堆石颗粒及其在空间中的分布，采用Fortran语言编写堆石体随机颗粒模型的生成程序SGG(StochasticGranule Generator)。

堆石体的堆石颗粒随机分布问题，与混凝土和土石混合体有所不同，颗粒没有水泥砂浆和土的包裹，而是依靠颗粒之间的相互接触来保持平衡，此外还要满足一定的空隙率，故处理起来比较困难。本发明在建立堆石体随机模型的时候，借鉴了混凝土粗骨料随机投放算法的思路，并提出了新的方法，提高了颗粒投放效率。

(2.1)建立堆石颗粒随机模型

随机数的产生是随机模拟技术的基础，本发明采用Mersenne Twister算法生成[0，1]区间均匀分布的随机数。考虑到堆石颗粒的形状及其在空间中的分布是随机均匀的，因此采用均匀分布模型模拟堆石体。影响堆石体随机颗粒模型的因素有：颗粒级配曲线(粒径)、颗粒的形状、位置坐标等。

根据颗粒级配曲线确定每组粒径的上下限，按粒径从大到小生成每组粒径区间内的颗粒直到颗粒含量满足级配要求。实际颗粒集合体是不会出现相交和侵入现象的，故还要判断颗粒间的相互侵入关系。

实际的堆石体颗粒是通过人工爆破、破碎而成，主要形态是凸多面体。为了能反映堆石颗粒的不规则形态，在不等边椭球内生成随机凸多面体颗粒，参见图2。生成方法为在某不等边椭球上随机布点，颗粒的顶点数随机取值，然后按照给定的算法连点形成凸多面体，具体方法如下。

(2.1.1)按照颗粒级配曲线生成不等边椭球的半径

a＝r₁+(r₂-r₁)rand1

b＝r₁+(r₂-r₁)rand2(1)

c＝r₁+(r₂-r₁)rand3

式中：a，b，c不等边椭球的三个极半轴的长度，r₁为某一粒径组的粒径下限，r₂为相应粒径组的粒径上限，rand1、rand2、rand3为区间[0，1]均匀分布上的独立随机数。

(2.1.2)为了保证生成颗粒的形状具有足够的随机性，颗粒的顶点数在[n_min，n_max]之间均匀分布：

n＝n_min+(n_max-n_min)rand    (2)

式中rand为区间[0，1]均匀分布上的随机数，n_min、n_max可根据颗粒的实际形状取值，本实施例中n_min取8，n_max取16。

(2.1.3)采用球坐标确定凸多面体的顶点：

z_i＝z₀+ccosθ_i

式中：x_i，y_i，z_i不等边椭球上第i个点的坐标；x₀，y₀，z₀是不等边椭球的球心坐标；θ_i是球坐标系中的天顶，在[0，π]区间内均匀分布；

是球坐标系中的方位角，在[0，2π]区间内均匀分布。

(2.1.4)在不等边椭球上随机布完点后，遍历每个顶点p_i，寻找与顶点p_i距离最近的顶点p_i1，然后再在剩下的顶点中找一个点p_i2，使得其余所有顶点均在p_i、p_i1、p_i2这三个点确定的三角形平面的同一侧。遍历完所有顶点后，删除具有相同节点的面。

(2.1.5)颗粒侵入判别在三维随机颗粒模型生成中至关重要，判别准则的好坏直接影响程序实现的难易与算法的效率和生成颗粒的质量。根据侵入关系的判别准则，如果两个凸多面体是相互分离的，进入步骤(2.1.6)；如果两个凸多面体相互侵入，回到步骤(2.1.1)重新生成凸多面体。

本发明三维颗粒侵入判别方法如下：当凸多面体A₁A₂A₃…A_n所有顶点在凸多面体B₁B₂B₃…B_n任意一个面B₁B₂B₃的一侧，而凸多面体B₁B₂B₃…B_n的形心在另一侧，那么可以断定这2个凸多面体是相互分离的。此种方法简单明了，避免了前述研究中判别方法中出现的不必要的麻烦。

由于所有颗粒均在不等边椭球面上内接而成，为了进一步加快颗粒侵入的判断，当两颗粒球心的距离大于两颗粒最大极半轴之和时，颗粒是互不侵入的。

(2.1.6)存储经上述步骤生成的全部凸多面体的几何信息；

(2.2)插入界面单元发生颗粒破碎

颗粒内部划分若干个细观单元，允许颗粒发生变形，颗粒的滑移、旋转、颗粒间的接触采用变形体离散元模拟，采用动态松弛法求解运动方程，使用显式时步步进的中心差分方法更新节点的位移以及单元应变。

由于颗粒破碎位置的不确定性，在颗粒内部的所有细观单元之间插入界面单元。通过界面单元的起裂、扩展和失效，实现开裂扩展的数值模拟。界面单元的法向、切向应力与张开、滑移变形之间的关系采用内聚力模型来描述：在加载的初始阶段，界面单元的应力与相对位移之间满足线性关系，随着加载的进行，界面单元的应力状态达到破坏准则后，界面单元的刚度逐渐下降，承载能力降低，当刚度降低到0时，界面单元完全失效，颗粒发生破碎，参见图3。采用带拉伸截断的摩尔库伦准则(Mohr-Coulomb)作为界面的破坏准则。当界面单元的法向应力达到抗拉强度后发生拉伸破坏，当界面单元的切向应力超过抗剪强度则发生剪切破坏，优先考虑拉伸破坏。

(3)利用堆石颗粒随机模型和强度劣化模型计算并绘制堆石颗粒流变的长期变形曲线

(3.1)将整个流变过程平均划分为N个时间步长ΔT，一般以若干小时或者若干天为一个步长；

(3.2)根据强度劣化模型计算堆石颗粒从0时刻开始每经过一个时间步长ΔT后新的强度；

(3.3)根据所得强度对堆石颗粒随机模型进行平衡迭代计算，耗时10^-5～10^-6s，求取并保存不同强度对应的静力平衡状态的变形值，包括轴向应变值和体积应变值；

(3.4)将得到的变形值用平滑的曲线串联，得到整个流变过程的长期变形曲线。

按照本发明的技术方案，进行了堆石料三轴流变模拟试验，试样级配采用双江口心墙堆石坝堆石料流变试验级配。试样尺寸φ300mm×600mm，最大粒径d_max＝60mm，孔隙率为30％，共生成8586个颗粒，采用二阶四面体网格离散为123343个实体单元，204491个界面单元。颗粒级配曲线参见图1，由于堆石体含有大量的小颗粒，若模拟完整的级配曲线，总的颗粒数目将会十分庞大，使计算效率大大降低，为此，本实施例的模拟试样将级配曲线进行截断处理，将粒径小于10mm的颗粒用粒径为10mm的颗粒进行等量替换。

表1为最终的细观参数。室内试验成果来自长科院所做的双江口堆石料三轴试验(参见文献：堆石料蠕变试验方法研究，左永振，程展林等，长江科学院院报)。

表1细观参数

模拟试验中围压1.6MPa，应力水平0.8，细观参数如表1所示，模拟试验中的应力路径与室内试验一致。根据双江口花岗岩堆石料母岩的长期强度试验，确定劣化模型参数a为0.145，n为4.26。模拟试验结果与长科院所做双江口堆石料的流变试验成果对比参见图7、图9。可以看出模拟试验的轴向流变与体积流变与室内试验成果规律上相似，数值上略有差别。这表明本发明抓住了堆石料流变的主要机理，适合模拟堆石料的长期变形这一复杂的、非线性演化问题。室内试验中初期流变变形略大于模拟试验结果，流变变形收敛较快。这是由于室内流变试验的试样中存在很多小颗粒，这部分颗粒受外界环境影响劣化较快，加速了流变变形的发展。

受外界环境影响，堆石料颗粒的劣化导致了宏观流变变形的出现。颗粒劣化的程度和劣化速率必然会影响流变变形的大小和收敛速度，为此本发明对颗粒强度劣化模型进行了参数敏感性分析。分析颗粒强度劣化程度影响时，固定反映劣化速率的参数n为5。分析强度劣化速率时，固定反映劣化程度的参数a为0.2。结果表明，随着颗粒强度劣化程度的增大，宏观流变变形量逐渐增加，强度劣化程度由0.1增大至0.4时，轴向流变变形从1.14％急剧增加至6.96％；颗粒强度的劣化速率对最终流变变形量和流变速率均有较大影响，颗粒劣化的越快，宏观流变变形量越大，流变变形收敛的越快。

同时，本发明采用模拟试验方法进行了不同母岩强度的堆石料流变试验。颗粒单轴抗压强度分别为60MPa、120MPa，其余参数见表1，强度劣化模型参数a为0.15，n为3。模拟试验结果表明，母岩强度对堆石料流变特性影响较大，在其他参数相同的情况下，颗粒强度越低，堆石料的轴向和体积流变变形均越大，颗粒强度由60MPa增大至120MPa时，轴向流变变形从2.33％减小至0.69％，参见图8、图10。这与工程实践中软岩堆石料的后期变形较大这一事实相符合。

Claims (3)

1.一种模拟堆石体长期变形的试验方法，其特征在于包括如下步骤：

(1)根据室内三轴试验获取的堆石体等时应力应变曲线，建立强度劣化模型

(1.1)取样，对堆石体试样进行三轴试验，绘制等时应力应变曲线；

(1.2)根据等时应力应变曲线建立时间变化与堆石颗粒的强度、劣化程度、劣化速度之间的等量关系，得到以下公式：

$f\left(\mathrm{\τ}\right)=f_0(1-\frac{\mathrm{a\τ}}{\mathrm{\τ}+n})$

式中：f(τ)表示τ时刻的强度；f₀为瞬时强度；a为反映劣化程度的参数；n为反映强度劣化快慢的参数；

(2)三维堆石颗粒集合的随机模拟

(2.1)建立堆石颗粒随机模型

(2.1.1)按照颗粒级配曲线生成不等边椭球的半径

a＝r₁+(r₂-r₁)rand1

b＝r₁+(r₂-r₁)rand2

c＝r₁+(r₂-r₁)rand3

式中：a，b，c不等边椭球的三个极半轴的长度，r₁为某一粒径组的粒径下限，r₂为相应粒径组的粒径上限，rand1、rand2、rand3为区间[0，1]均匀分布上的独立随机数；

(2.1.2)在[n_min，n_max]之间均匀分布颗粒的顶点数：

n＝n_min+(n_max-n_min)rand

式中：rand为区间[0，1]均匀分布上的随机数，n_min、n_max根据颗粒的实际形状随机取值；

(2.1.3)采用球坐标确定凸多面体的顶点：

z_i＝z₀+ccosθ_i

式中：x_i，y_i，z_i不等边椭球上第i个点的坐标；x₀，y₀，z₀是不等边椭球的球心坐标；θ_i是球坐标系中的天顶，在[0，π]区间内均匀分布；

是球坐标系中的方位角，在[0，2π]区间内均匀分布；

(2.1.4)在不等边椭球上随机布完点后，遍历每个顶点p_i，寻找与顶点p_i距离最近的顶点p_i1，然后再在剩下的顶点中找一个点p_i2，使得其余所有顶点均在p_i、p_i1、p_i2这三个点确定的三角形平面的同一侧；遍历完所有顶点后，删除具有相同节点的面；

(2.1.5)根据侵入关系的判别准则，如果两个凸多面体是相互分离的，进入步骤(2.1.6)；如果两个凸多面体相互侵入，回到步骤(2.1.1)重新生成凸多面体；

(2.1.6)存储经上述步骤生成的全部凸多面体的几何信息；

(2.2)插入界面单元模拟颗粒破碎

在颗粒内部的所有细观单元之间插入界面单元；当界面单元的刚度降低到0时，界面单元完全失效，显示颗粒发生破碎；

(3)利用堆石颗粒随机模型和强度劣化模型计算并绘制堆石颗粒流变的长期变形曲线

(3.1)将整个流变过程平均划分为N个时间步长；

(3.2)根据强度劣化模型计算堆石颗粒从0时刻开始每经过一个时间步长后新的强度；

(3.3)根据所得强度对堆石颗粒随机模型进行平衡迭代计算，求取并保存不同强度对应的静力平衡状态的变形值；

(3.4)将得到的变形值用平滑的曲线串联，得到整个流变过程的长期变形曲线。

2.根据权利要求1所述的模拟堆石体长期变形的试验方法，其特征在于：所述侵入关系的判别准则为当凸多面体A所有顶点在凸多面体B任意一个面的一侧，而凸多面体B的形心在该面的另一侧，则这两个凸多面体是互不侵入的。

3.根据权利要求1所述的模拟堆石体长期变形的试验方法，其特征在于：所述侵入关系的判别准则为当两颗粒球心的距离大于两颗粒最大极半轴之和时，则这两个凸多面体是互不侵入的。

Images