某大堤空洞位置与弹性波中纵横波关系的数值模拟研究

2013-12-23张振华袁梦玲

三峡大学学报(自然科学版) 2013年2期

张振华钟高周杨袁梦玲，3

（1.三峡大学土木与建筑学院，湖北宜昌 443002；2.水利部堤防安全与病害防治工程技术研究中心，郑州 450033；3.黄河上游水电开发有限责任公司大坝管理中心，西宁 810000）

堤防的主要功能是在某一范围能够抵御一定防洪标准的洪水侵害，是防洪体系中的一个重要组成部分，主要起着抵御洪水、保护人民生命财产安全的作用.我国现有的大堤普遍存在着堤身填筑质量不高、堤基渗透变形、沉降变形和堤岸失稳等问题，一旦发生险情将严重威胁到人民的生命及财产安全，只有及时地发现堤防工程的缺陷，才能防患于未然［1］.而堤防工程堤身内部的缺陷很难直接被发现，且只有通过较为先进的勘探方法才能快速准确的判断堤身内部的缺陷类型和空间位置.由于历史原因，黄河下游堤防的堤身、堤基及近堤范围内存在多种隐患，1950年6月，黄委制定了《堤防大检查实施办法》，据此河务部门广泛发动沿河群众，普遍开展堤身隐患探查，通过探查共查出隐患7 605处，其中獾狐洞穴559处，战沟残缺175 处，坑穴50 处，井穴、暗洞、红薯窖198处.1952年，山东堤防探查发现碉堡66处，军沟206条，防空洞101处，其中齐河县发现一条军沟长20多m，宽2m 多，洞内有尸体和铁锅等物.1953年在河南武陟厂探查出一组獾洞，位于堤顶下3～9m，洞穴横三竖四，上下交错，全长达300m.大量隐患的发现，加深了对消灭隐患重要性和长期性的认识，从而将普查隐患工作定为制度［2］.

我国堤坝隐患检测始于20世纪80年代，近几十年来，堤防检测方法由人工探视法到破损法，再发展到无损检测法.在堤身无损检测方面使用较为普遍的方法有地震勘探法、电法、水声勘探法、探地雷达法、弹性波测试法、层析成像法、放射性测量法和综合测井法，各种方法各有优缺点，其中地震勘探法在很多物探领域得到了广泛的应用，一般地，是利用地震波所产生的多种弹性波中的一种波来检测堤身内部的缺陷［3－4］.

对于这些方法，主要是通过实验的方法来研究，然而，不管是做现场实验还是做室内实验，都会受到一定条件的限制，且要耗费较大的财力和物力；而数值模拟方法也是研究的手段之一，其特点是可在室内利用软件在计算机上重复地进行，可节省大量的人力和财力.因此，本文基于Flac3D 方法，对大堤堤身在有空洞条件下的动力响应进行计算，研究黄河大堤模型在动力响应下横波和纵波的传播特征，获得有空洞的堤身模型在无阻尼条件下堤身空洞位置与弹性波中的纵波与横波的关系差异，为堤防工程空洞缺陷的检测提供参考.

1 纵波与横波在堤防中传播基本特征

当敲打地面时，地面会振动，这些振动向远处传播就形成了波动（见图1）.在这些波动中包括体波和面波，常见的体波有纵波和横波，都是从震源呈放射状向外在介质中传播，其中纵波的介质质点的振动方向与传播方向一致，而横波的介质质点的振动方向与传播方向垂直.由于纵波的传播速度最快，发生振动后会首先到达观测点，所以又叫P 波（Primary Wave），而横波速度较慢，在纵波之后到达，所以又叫S波（Secondary Wave）［5］.横波在传播过程中，媒质内部会发生剪切变形，固体有切变弹性，可以传播横波，液体和气体都没有切变弹性，不能传播横波.纵波在传播过程中，介质中的各部分时而受压，时而受拉，每时每刻都在发生体积改变；而当固体体积改变时，其内部会产生弹性恢复力，液体和气体也会产生弹性恢复力，所以纵波在这3种状态的介质中都可传播.

图1 堤防地震波法检测示意图

2 动力计算

本文应用Flac3D 软件中提供的动力计算模块对模型进行动力计算，对模型输入地震波荷载，模拟研究有空洞的某段黄河大堤堤身在地震荷载作用下的纵波和横波的传播特征.

2.1 数值网格模型及材料参数

以黄河大堤的某段堤身为研究对象，模拟大堤常见缺陷中的空洞.根据堤防基本断面设计标准确定黄河大堤横截面尺寸［2］，堤顶宽取12m，堤身高取15 m，堤身长取100m，大堤背水坡度为1∶1.5，临水坡度为1∶3.采用六面体八节点单元对研究的堤段进行三维数值网格划分，数值网格模型共划分单元数356 800个，节点数383 307个，数值网格模型具体如图2所示.

图2 黄河大堤数值网格模型

黄河大堤堤身中的空洞尺寸小到一两米，大到数10m［2］，本文拟在堤顶表面以下深度h 为3m、6m、9 m 和12m 的位置分布设置边长为1m 的立方体的空洞，进行大堤空洞位置与横波及纵波传播的对应关系研究，空洞的具体位置如图3所示.

图3 堤身模型Ⅰ－Ⅰ剖面图

堤身材料设为均质土体，土体密度取为2 200kg／m3，弹性模量取为180MPa，泊松比为0.3，重力加速度为10m／s2.

2.2 边界条件的设置

对图2中大堤模型左右2个堤身侧面及底面静力学边界为法向约束边界，堤身上表面及前后两个坡面为自由边界.对图3 数值模型中堤身的空洞采用Flac3D 软件中提供的空模型（null模型）进行模拟.对于Flac3D 软件中的静态边界条件可以用来吸收边界上的入射波，相当于在模型的法向和切向分别设置自由阻尼器，所提供的法向和切向黏性力可以分别由以下两式求出.

式中，vn和vs分别表示模型边界上法向和切向的速度分量，ρ是介质的密度，Cp和Cs分别表示P波和S波的波速.

对图2中大堤模型各侧面的纵波和横波传播边界设置为自由场运动边界（自由场边界提供了与无限场地相同的效果［6］）.在动力边界施加后，原先施加在这些边界上的静力边界条件将会自动被去掉，如果需要在静态边界上输入动力荷载，那么就只能输入应力时程，也可将加速度时程和速度时程通过转换公式（3）和（4）变成应力时程，然后再施加到静态边界上：

式中，σn和σs分别表示施加在静态边界上的法向应力和切向应力；系数“2”表示施加的能量中只有一半是作为动力向上传播的，而另一半则是向边界下部传播的；其它符号意义同式（1）和（2）.

2.3 阻尼的设置

Flac3D 有限差分软件中提供了3 种阻尼形式，一种是瑞利阻尼，一种是局部阻尼，还有一种是滞后阻尼.而实践证明，用瑞利阻尼计算所得到的加速度响应规律与实际中最相符，因此大堤堤身在纵波和横波传播过程中的阻尼设置采用瑞利阻尼.根据经验的做法，大堤堤身土体的阻尼比为0.5%［6］；经计算，大堤模型完成一个振荡周期需要的时间约为1.3s，据此可得大堤的自振频率为0.78Hz.

2.4 地震波的施加

Flac3D 方法可以通过在模型边界或内部节点施加动载荷来模拟材料受到外部或内部动力作用下的反应.程序允许的动力载荷输入可以是加速度时程，速度时程，应力（压力）时程，集中力时程［6］.本文所计算的模型材料类型是均质土体，可采用地震速度时程直接输入到模型顶部，并选择横波和纵波分别作为施加的地震波.

在数值模拟中，纯剪切力源就可以模拟横波在模型中的传播，因此在本文采用纯剪切力源来模拟横波（如图4（a））.横波产生的震源的施加在1－1剖面上，深度为堤顶表面以下5m，具体坐标为（0，50.5，10）.由于在均匀各向同性介质中，爆炸震源产生的是纯纵波，因此本文采用爆炸震源模拟纵波在堤身中的传播过程.爆炸震源的模拟是在网格单元的某些节点上加作用力，使它们的作用效果相当于球腔爆炸震源的对称径向作用力（如图4（b））.纵波产生的震源的施加在1－1剖面上，深度为堤顶表面下5m，具体坐标为（0，50.5，10）.

图4 横波和纵波震源示意图［7］

2.5 计算结果及分析

采用Flac3D 方法，利用动荷载频率为10Hz（其它大于10Hz频率的计算结果有类似的规律，但10 Hz栋荷载频率的规律最为明显，因此本文仅以动荷载频率为10Hz进行说明）纯剪切震源和爆炸震源分别模拟横波和纵波在大堤模型上传播，可获得位于1－1剖面上、水平位置为x＝40m、竖直位置位于堤顶下深度为3m、6m 和9m 处时的横波和纵波传播时堤顶表层质点的振动图像，具体见图5和图6.

图5 空洞位于3种不同深度时横波传播时域波形图

鉴于对每一种材料而言都对应着一个固定的横波或纵波波速，因此在同一种均质介质中，横波或纵波速是一个定值，也就是说质点首先出现波谷的位置所连成的斜线的斜率就是该列横波或纵波的波速.图5和图6中堤顶表面各测点第一次出现波谷连线的斜率即为横波河纵波的波速.

图6 空洞位于3种不同深度时纵波传播时域波形图

图5中横波波速发生变化的位置已经用红色圆圈标出.从图5中可知，图5（a）中x 坐标位置由42～46m 处波速发生了变化，图5（b）中x 坐标位置由43～48m 处波速发生了变化，图5（c）中x 坐标位置由45～50m 处波速发生了变化.由此可以判断频率为10Hz的纯剪切震源产生的横波在大堤中的传播特征为表层质点振动发生延迟的x 坐标范围均在42 m 左右，而本次计算设置的空洞的水平位置为x＝40 m.因此，可以通过表层质点振动发生延迟的x 坐标范围，初步判断出堤身空洞的水平位置范围.

由图6中纵波波速发生变化的位置已经用粉红色圆圈标出.从图6可以看出，图6（a）中x 坐标位置为42m 处波速发生了变化，图6（b）中x 坐标位置为44m 处波速发生了变化，图6（c）中x 坐标位置为46 m 处波速发生了变化.由此可判断频率为10Hz的爆炸震源产生的纵波在大堤中传播特征为随着深度在增加，表层质点振动发生延迟的x 坐标也逐渐推后.因此，结合纵波在有空洞大堤中的传播特征及其反射的规律（当纵波在堤身中传播遇到空洞时，其波速会立刻减小，并在空洞与堤身的分界面上发生反射，经反射的横波使表面质点的振动延迟），可以初步判断出空洞深度的范围.