人工黄土坡面重力侵蚀数值模拟

2020-08-21王春振赵华荣

中国农村水利水电 2020年8期

金鑫，刘喆，王春振，宋颖，赵华荣，王迪，李霞

(1.桂林理工大学广西环境污染控制理论与技术重点实验室，广西桂林 541004；2.桂林理工大学岩溶地区水污染控制与用水安全保障协同创新中心，广西桂林 541004

0 引言

当地表土石在地震、降雨、径流和海浪、风、冻融、人工采掘等任何一种应力作用时，在重力作用下，而产生迁移和堆积的一种侵蚀过程称为重力侵蚀[1]。重力侵蚀是我国黄土高原地区土壤侵蚀的三种主要类型之一[2-5]，但由于影响重力侵蚀的因素较多，且其发生具有随机性和突发性，因而较难准确的测定产沙量[6]。相较于对水力侵蚀过程机理的认识，对重力侵蚀的研究还较为薄弱，侵蚀机理尚未完全明确[7]。近年来，随着研究技术和方法的不断应用，对重力侵蚀的认识逐渐加深。例如，谭万沛等[8]探究了不同的坡度条件对重力侵蚀的影响，指出重力侵蚀的发生形式和侵蚀产沙量在不同坡度条件下存在较大差异。Gabet[9]发现重力侵蚀的出现次数在不同下垫面植被类型条件下有较大波动，且重力侵蚀的发生次数在下垫面植物根系具有较深的深度及较高的强度时，而出现了明显减少。郑书彦[10,11]等对重力侵蚀的发生机理分别从滑坡等土壤侵蚀角度进行了探讨。朱同新[12]等研究指出在不同的空间位置处和不同的时间阶段内，重力侵蚀发生强度和类型是不同的。此外，王光谦等[13]研究了沟坡冲蚀条件下重力侵蚀的力学机制，并通过力学方法建立了坡沟重力侵蚀的概化力学模型，李容全[14]通过14C法计算了在重力侵蚀下沟谷后退速度，叶浩[15]通过GPS研究了砒砂岩地区重力侵蚀引起的沟缘线后退速度等。朱同新等[16-18]等对黄土高原重力侵蚀发生现象进行了观察和统计分析，杨吉山[19]通过对桥沟小流域的研究指出不同的沟道发育阶段，重力侵蚀的特点也不尽相同。在大多数学者的研究中，重力侵蚀主要是在不同的地形、地貌、植被覆盖等条件下进行研究。但从沟坡系统力学稳定性的角度来看，由于沟坡系统重力侵蚀具有随机性和多发性，基础观测资料相对较少，相关研究相对薄弱[20]。因此，本文通过Rhino和FLAC 3D数字模拟软件的二次开发，建立人工黄土坡面侵蚀概化模型，并对黄土坡面重力侵蚀破坏区进行研究，为今后的水土流失相关研究提供参考。

1 实验材料与方法

1.1 实验概况

本实验用土为陕西榆林岔巴沟黄土，属于第四纪时期形成的土壤堆积物，含多量钙质或黄土结核，在干燥时较坚硬，被水流润湿后，更容易剥落和遭受侵蚀，持水性一般，土体垂直节理发育，平均容重为1.15 g/cm3。本实验用土地区选自岔巴沟黄土侵蚀典型代表地区作为实验用土采集地，采用机械挖掘表层土样后，将土样装车运回。土样采集回来放置于库房将土壤经风干后，过1 cm的筛网，剔除天然杂草、石砾以及较大的土壤结块。为保证土槽的透水性及透水状况接近于天然坡面，先在土槽底部铺上10 cm细沙，并盖上透水纱布。根据土槽的体积及所需容重进行分层填土，以10 cm为一层并压实，每填完一层把表面刮毛以减少各层间的差异性。完成一组实验后，将表层10 cm土壤刨松，然后填入新土进行下一组实验。

每次正式降雨实验前一天，选用30 mm/h降雨强度进行前期预降雨，降雨至坡面有积水出现，并用塑料布覆盖并静置24 h。前期降雨的目的一是保持坡面下垫面有相同的前期土壤含水量，二是减少下垫面条件的空间变异性。

实验中降雨系统采用全自动不锈钢模拟降雨器，该降雨装置配有旋转下喷式喷头，该设备采用潜水泵提供动力，其扬程为50 m，降雨高度为6 m，雨强可调节范围为10～200 mm/h，雨滴直径为0.1～6 mm。实验采取两个坡度(20°、25°)，一个雨强(60 mm/h)。每场次降雨持续时间为90 min，每组实验分为4个场次进行连续性降雨，因此使坡面细沟得到充分的发育，便于观测在重力侵蚀作用下坡面细沟不同发育阶段的形态特征。

1.2 人工黄土坡面重力侵蚀数字分析原理及建模方法

1.2.1 FLAC 3D数值模拟基本原理

FLAC 3D(Fast Lagrangian Analysis of Continua)是由Itasca公司研发推出的连续介质力学分析软件，内置有静力、动力、蠕变、渗流、温度5种计算模式，各种模式间可以相互耦合，以模拟计算复杂的工程力学行为。作为有限差分软件，相对于其他有限元软件，在算法上有以下优点：FLAC 3D采用混合离散法[21]来模拟材料的塑性破坏和塑性流动，较通常采用的离散集成法更为准确合理；FLAC 3D采用显示差分法求解微分方程，可以方便地求得应力增加、不平衡力并跟踪系统的演化过程[22]。

在FLAC 3D中采用混合离散方法，区域将被离散为常应变多面体单元[23]，在运算过程中，每个多面体又进一步离散为以该多面体顶点为顶点的常应变四面体，所有变量在四面体上计算，取多面体单元内四面体应力、应变的加权平均值作为多面体单元的应力、应变值。如图1所示，对任意一个四面体，节点编号为1～4，第n面表示与节点n相对的面，设其内任一点的速率分量为vi，则可由高斯公式得

(1)

式中：V为四面体的体积；S为四面体的外表面；nj为外表面的单位法向向量分量。

对于常应变单元，vi为线性分布；nj在每个面上为常量，由式(2)可得：

(2)

式中：l为节点l的变量；(l)为面l的变量。

图1 四面体Fig.1 Tetrahedron

在FLAC 3D中计算的对象以节点为主，将力与质量集中在节点上，通过运动方程进行求解，公式如下。

(3)

FLAC3D由速率求解某一步时长的单元应变增量，公式如下。

(4)

对于静态问题，在不平衡力中加入非黏性阻尼，使得系统振动逐渐衰减至达到平衡状态。此时式3变为：

(5)

其二，常见的书面体育文本主要包括体育新闻报道、体育教学资料、体育营销文案、体育学术论文等，所涉体育项目众多（罗永洲，2012）。尽管如此，体育类文本也呈现出一定共有特征。体育文本往往牵扯到国际赛事、国际组织、运动员姓名、比赛规则、技术要领等方面，具有专业性强，专有名词、专业术语、缩略语、合成词繁多等特征。因此，体育文本的翻译过程需要处理大量的特殊词汇。鉴于体育文本的这一特点，相对于其他类型的翻译而言，体育类翻译对译员的专业知识和词汇的积累要求更高。

(6)

式中：∝为阻尼系数，默认值为0.8。

由上可知FLAC 3D的计算过程，如图2所示。

图2 FLAC 3D计算流程Fig.2 Calculation process of FLAC 3D

1.2.2 Rhino三维建模软件介绍

Rhino是美国Robert Mcneel公司开发的一款基于NURBS为主体结构的专业三维建模软件，采用自由曲面的建模技术和特征实体的操作模式[24]。相较于POLYGON(多边形)，NURBS可用极少的控制点创建具有高精度的复杂曲线并建立曲面模型，可精确的实现自由造型产品模型制作，同时支持ACIS、Parasolid、3ds、dwg以及点云资料等多种格式文件的输入输出。在Rhino建模之前为确保模型精度，在模型核心参数公差值的选择上，采取Absolute tolerance默认的绝对公差值，默认Absolute tolerance值为0.01 mm。

1.2.3 复杂地形的FLAC 3D建模方法

虽然FLAC 3D专为岩土工程力学开发，内置有丰富的弹、塑性材料本构模型。但是在建立复杂计算模型时，仍是以命令驱动模式为主要输入模式，这种方式对于模型的建立过于繁杂，因此本研究采用3D造型软件Rhino结合FLAC 3D内置的fish语言在初始单元模型基础上编写了前处理程序，实现了对复杂多层地形建模的二次开发[25]。

首先使用Rhino读入原始侵蚀坡面点云数据，采用MeshPatch指令进行网格补丁，建立三角网格，对网格进行修剪、调整网格坐标原点位置、调整模型单位为米，生成STL文件类型的地表模型[图3(a)]。在FLAC 3D中通过fish语言程序提取STL文件的地表高程信息，然后使用generate zone brick生成六面体实体网格，采用全局坐标系定义坐标轴x、y、z方向网格单元数目及相邻单元尺寸大小比率，然后调用generate topography geomset命令填充地表模型与新建六面体实体网格之间的间隙，建立接触面。采用Rhino软件较好地解决了在FLAC 3D中前期复杂侵蚀地形建模困难的问题，成功实现建模。仿真实例表明，所建立的三维模型能够真实反映侵蚀坡面地形，模拟效果良好。

1.3 人工黄土坡面重力侵蚀的有限元数值模拟

1.3.1 人工黄土坡面概化模型及有限元计算模型

图3(a)为采用Rhino软件绘制的侵蚀坡面地形图，图3(b)为采用FLAC 3D建立的侵蚀坡面概化模型。该模型底部和四周设为固定约束边界，坡面设为自由边界[25]。该模型土层设置为单一土层，平均厚度为1 m，长度设为4 m，宽度设为1.2 m。坡面坡度为20°、25°，20°模型有限元网格共有节点445 511 个，单元40 万个；25°模型有限元网格共有节点445 511 个，单元40 万个。

图3 20°及25°人工黄土坡面地形图及概化模型Fig.3 Topographic maps and generalized models of 20° and 25° artificial loess slopes

1.3.2 黄土物理力学指标

根据边坡工程经验、现场资料分析、现场及室内岩土物理力学试验和有限元计算的需要[25]，模型土层材料物理力学参数的具体特征取值见表1。

表1 黄土坡面性质参数Tab.1 Properties of loess slope

2 结果与分析

2.1 黄土坡面重力侵蚀发育情况

分别在坡面0.9、1.6、2、3 m处设置监测点，观测坡面上中下3个部位处位移变化情况。模型运行到24 000 步时，计算最大不平衡力低于系统设置的默认值(10-5)，模型计算停止。从图4中可以看出，各点的位移变化趋势相同，随着步长的增加，各点位移增大，然后降低至某一值后保持稳定。观察图4位移曲线变化规律可划分为不同的侵蚀发育时期，将开始至位移最高点阶段划分为侵蚀发育期，从最高点降至稳定值阶段划分为侵蚀发育成熟期，从稳定值以后划分为稳定期，由图4可见，坡面下部侵蚀发育完成时间要快于中上部侵蚀发育，这与实验过程观察结果一致。且从图4中可知，在坡面上部0.9 m处整体位移最大，坡面下部3 m处整体位移最小，当坡度增大时，中下部位监测点位移均有一定程度增大，其中1.6、2、3 m处增加幅度分别为10.7%、55.4%、33.4%。上部监测点位移减小，减小幅度为4.5%。分析原因为，坡度增加导致坡面径流速度增大、坡面承雨面积减小，从而在一定程度上加剧了坡面中下部侵蚀发育程度，减缓了坡面上部的侵蚀发育。

图4 监测点位移情况Fig.4 Displacement of the monitoring point

2.2 黄土坡面位移分布规律

在最大不平衡力低于系统默认值时，系统达到平衡状态，得出模型3个方向的位移大小和坡面整体位移、竖直方向位移和水平方向位移分布图(图5、图6、图7、图8)。从整体位移云图来看，位移最大部分出现在坡顶及坡面上部；竖直方向位移云图与总体位移云图分布规律基本一致，且竖直方向最大位移数值也出现在坡顶及坡面上部，由此可知坡面位移是以向下侵蚀滑落为主。由图7见，最大水平位移出现在坡面下部沟坡地带，而坡面其他位置的水平位移基本为零，这表明相对于其他位置，下部出现水平位移的沟坡处易在水平方向发生侵蚀坍塌或滑坡，而朝沟底方向滑动。由图8可知，纵向最大位移出现在坡面中下部，这与侵蚀实验观测结果一致，侵蚀实验中坡面中下部侵蚀发育最为剧烈。

图5 20°及25°坡面整体位移图Fig.5 Overall slope displacement of 20° and 25°

图6 20°及25°坡面竖直方向位移图Fig.6 Vertical displacement of 20° and 25° slope

图7 20°及25°坡面水平方向位移图Fig.7 Horizontal displacements of 20°and 25°slopes

图8 20°及25°坡面纵向位移图Fig.8 Longitudinal displacements of 20°and 25°slopes

2.3 黄土坡面应力场分布规律

FLAC 3D分析计算出坡面模型在平衡状态时所受到的应力大小及分布规律。如图9、图10中所示(FLAC 3D中定义为以拉应力为正，压应力为负)，应力分布基本一致，从坡面应力图分布来看，未出现拉应力，因此坡面受力基本以压应力为主，即为坡面受到破坏时主要以压剪形式产生的破坏模式为主。从图9可以看出，边坡内部，最大主应力(即第一主应力)方向与水平面夹角逐渐减小，且随着坡度的增加土体内部受到的压应力也随之增大，表明深层土体主要受铅直方向的压应力，表现为受压屈服。从图10分析可知，由于坡度增大，导致应力集中，土体内部主应力及最大受压屈服区体积增加，因此，在一定程度上坡度的增加，加大了土体的受力，从而导致加剧了重力侵蚀的发生。

图9 坡面第一应力分布图Fig.9 Distribution of the first stress on the slope

图10 坡面第三应力分布图Fig.10 Distribution of the third stress on the slope

2.4 黄土坡面塑性区分布规律

图11为模型计算达到平衡状态时坡面塑性状态指示图，图中主要获取剪切屈服区域(shear)和张拉屈服区(tension)[26]。每个屈服区函数赋有两种状态：now和past，其中now表示该单元在本次计算步数中正处于屈服面上，past表示曾处于屈服面上，但现处于弹性范围，因此只有处于now状态单元才对模型的破坏起作用[27]。由图11可见，张拉屈服区主要分布在坡面顶部及坡面中上部位，说明在坡顶及坡面中上部位发生张拉破坏的可能性较大，张拉破坏较为严重；剪切屈服区域分布较为广泛，其主要分布在坡面上下部沟坡区域及土体内部，说明这些区域易发生水平剪切变形，且当坡度增加时，坡面上部位剪切屈服区减少，张拉屈服区有所增加。

图11 坡面塑性状态分布Fig.11 Distribution of slope plastic state

从图11中可知，以20°坡面为例，根据第一场次模型平衡状态时的塑性区分布，对应第二场次连续性降雨结束后侵蚀地形，可以看出坡面中下部出现明显侵蚀细沟且随着降雨进行，细沟开始溯源侵蚀，侵蚀地形符合第一场次模型平衡状态时的塑性区分布规律。同时对三次连续性降雨坡面进行受力分析，发现本场次侵蚀地形基本符合上一场模型平衡状态时的塑性区分布情况，可以说明FLAC3D在一定程度上能够对侵蚀地形发育起到预测的作用。

通过编程FISH语言，调用命令流对20°、25°两个模型塑性屈服区体积进行计算，结果如表2所示。可以看出，在黄土坡面上塑性屈服区主要以剪切屈服区为主，其中在20°坡面上，剪切屈服区占总屈服区体积的87.9%，张拉屈服区占总屈服区体积的12.1%；在25°坡面上，剪切屈服区占总屈服区体积的85.4%，张拉屈服区占总屈服区体积的14.6%，说明坡面主要以剪切破坏模式为主。当坡度增加至25°时，剪切屈服区体积增加了0.271 m3，张拉屈服区体积增加了0.166 m3，总屈服区体积增加0.437 m3，相较于20°坡度下剪切屈服区、张拉屈服区和总屈服区体积分别增加了7.6%、34.08%、10.8%。说明随坡度增加，张拉屈服区破坏体积占总屈服区破坏体积增大，但剪切塑性屈服区的体积占全部屈服区体积的大部分，因此重力侵蚀破坏仍是以剪切破坏模式为主。