地震触发海底滑坡物理模型及数学模型试验分析

李佳林冯秀丽刘 隆

(1.中国海洋大学海洋地球科学学院,山东 青岛 266100;2.自然资源部第二海洋研究所,浙江 杭州 310012)

南海北部陆坡区,属于张裂性大陆边缘,并处于环太平洋火山地震带,该区域构造运动剧烈,地震频发[1-3]。近几十年来,为开发海上油气资源,在南海地区搭建了大量的海上工作平台,同时也面临着各种海洋地质灾害的威胁,其中之一便是海底滑坡[4]。海底滑坡通常规模巨大,运动速度快[5-6],且触发因素复杂,而地震正是触发海底滑坡最普遍也最具破坏性的因素之一[7]。研究海底滑坡的地震响应及动力学特征,对于保障人类海上正常生产施工作业具有重要意义。

海底滑坡具有骤发性和不可预测性,且海洋观测技术限制,因此实地监测海底滑坡困难重重[8],所以对海底滑坡的研究多为室内分析。Coleman和Garrison[9]率先通过海底钻孔与室内分析确定了海底滑坡现象的存在依据。Milne[10]调查了1616—1861年间300余次海底地震作用,并率先指出地震可引发水下滑坡。地震作用产生的超孔隙水压力,使得土体有效应力降低,土体强度降低,是导致海底滑坡的重要因素[11]。

在海底滑坡室内分析方面,主要有物理模型试验和数值模拟稳定性评价两种方法。Schwarz[12]和Mohrig等[13-14]利用水槽模型并采用倾倒法研究了海底滑坡的运动形态和速度;Ilstad等[15-17]和Elverhoi等[18]则基于水槽试验分析了黏粒含量对滑坡的影响。海底斜坡稳定性评价主要有极限平衡法、数值分析法、极限分析上限法三种分析方法。

本文采用水槽模型试验方法,依据南海北部陆坡工程环境构建模拟陆坡模型,探讨引发陆坡失稳的临界地震条件等,并基于ANSYS有限元强度折减法,分析模拟陆坡在不同地震加载下的稳定系数变化。综合两部分试验结果,分析地震对于海底滑坡触发的影响作用。

1 水槽模型试验分析

1.1 试验设备及模型构建

本次水槽模型试验采用自制加工玻璃水槽和水平振动台进行(图1)。水槽长1.4 m,宽0.4 m,高1.0 m。槽边使用不锈钢材料进行焊接包裹加固,以保证水槽振动条件下的牢固性,后侧不锈钢板下方留有出水口,方便每组试验结束后排水。振动台台面尺寸为1 m×1 m,振动频率范围为1～10 Hz,额定水平加速度为4 m/s2。

图1 试验水槽及振动台Fig.1 Test tank and shaking table

本次试验用土为东营黄河口周边水库的砂质粉砂土样,平均粒径约4.4Φ,黏粒质量分数约为3%,密度约为1 840 kg/m3,土体力学性质方面,该土样黏聚力约为3.5 k Pa,内摩擦角约为25°。相较于南海北部陆坡土体,该土样整体粒度较粗,砂土含量较高,黏性土含量较少,土体黏聚力偏低,内摩擦角略大[19-23]。为了使模拟陆坡土体性质与真实海底一致,更接近南海北部陆坡的真实情况,选用马来西亚高岭土与原试验土样按不同比例搅拌混合。由于高岭土属于黏土,平均粒径较细,且饱和时黏性较高,混入高岭土可以降低试验土体的平均粒径,增加黏粒含量,土体黏聚力相应增高,并且当黏粒含量较高时,土体液化趋势增强,更易发生“滑水现象”,其内摩擦角也相应降低。因此,混入高岭土使得原试验土样土体性质与南海北部陆坡土体性质更加接近,使得模拟试验更有意义。

南海北部各段陆坡平均坡度大致在2°～5°,部分区域整体坡度小于1°[7],因此,模拟陆坡的坡度铺设也尽量与海底情况一致。针对不同高岭土含量及斜坡坡度设置大量预实验,加载不同强度的地震并观察现象,在保证试验顺利进行且有效试验组数充足的情况下最终确定本次试验土样中高岭土混入量为18%,斜坡坡度为5°。另外,海洋地球物理勘探资料显示,在南海北部陆坡区域的海底表层有一层厚度约1～20 m 的软弱结构层[24]。软弱层通常抗剪强度较低,更易受扰动而失稳,并且会呈现一个良好的滑动面。为了使模拟滑坡现象更加明显,也更好地与实际海底情况相贴合,配置高岭土混入量70%的土样作为软弱层。这样本次试验确定模拟陆坡坡度为5°,并共用到3组土样:原粉砂土样(土样Ⅰ)、含18%高岭土土样(土样Ⅱ)和含70%高岭土土样(土样Ⅲ)。3种土样粒度参数及力学参数见表1。

表1 三组土样粒度及力学参数Table 1 Grain size parameters of three soil samples

利用3组土样铺设模拟陆坡(图2),其中,含18%高岭土的土样Ⅱ作为陆坡主要滑动部分;白色部分为软弱层土样Ⅲ,提供滑动面的作用;基底为原试验土样Ⅰ,该部分对试验影响较小。陆坡最高点土层厚度为0.3 m,最低点土层厚底为0.2 m,实际坡顶距坡底垂直距离为0.1 m。斜坡倾角约为5°,中间软弱层厚度为0.02 m。

图2 模拟陆坡铺设Fig.2 The laying of simulated slope

水槽右侧留出0.4 m 的空间,利用挡板与土样间隔开,作为滑坡的沉积区,用来观察比较每组试验滑坡体运移量,反映滑坡的规模与剧烈程度。

向陆坡铺好的水槽中注水。通过预实验发现,水深超过0.5 m 时,在加载较大的振动参数时会造成水体溢出,水深小于0.5 m 时,振动带来的波浪效应对斜坡土体破坏干扰明显。因此确定注水深度为0.5 m,实际坡顶距水面0.2 m。本次试验为相对理想模型,即不考虑波浪等水动力作用的影响,以及水槽空间限制而产生的边壁效应,重点讨论地震对斜坡稳定性的影响。

1.2 模拟陆坡振动分析

1.2.1 滑坡触发临界地震条件

根据振动台工作参数,设置不同振动频率与加速度的组合。其中频率分别为1～5 Hz,加速度分别为0.1g～0.4g(g为重力加速度),除去工作台性能限制无法实现的组合外,共确定16组不同振动参数的实验。分别按照上述方法铺设陆坡,并加载水平方向每组振动,时长2 min,记录试验现象见表2。

表2 振动试验现象Table 2 Vibration test phenomenon

当振动频率＜3 Hz时(试验1～试验5),即使加速度增加至控制台的最大值,斜坡表面无明显液化流动或是块体位移现象,只有表层少量土体颗粒进入上层水体悬浮,水体略微产生浑浊(图3);而当控制振动频率≥3 Hz且频率固定时,发现加速度过小时仍然出现类似情况,对应于试验6 和实验10。当加速度≥1.96 m·s-2(即0.2g)时,才逐渐在斜坡表面产生滑动现象。随着振动强度的增加,斜坡破坏更加严重,沿斜坡向下滑动规模增大,更多的土体沿滑动面滑向右侧沉积区。由如图4所示的振动频率、加速度最大组(试验16),斜坡底部的滑坡形态可发现表面浊流、碎屑流向右卷动位移现象十分明显,右侧挡板处几乎堆积满滑塌下来的试验土体。

图3 较小振动参数下斜坡表面情况Fig.3 Slope surface conditions with small vibration parameters

图4 试验16斜坡底部滑坡形态Fig.4 Slope bottom landslide morphology of test 16

因此,确定本次地震触发滑坡物理模型试验在斜坡坡角5°、土样饱和密度1.85 g·cm-3、黏粒含量约5%并且软弱层存在的情况下,滑坡发生的临界地震条件为:振动频率为3 Hz,振动加速度为1.96 m·s-2。

由于部分组别试验现象较为重合,后续数据分析结果也基本吻合。本文精选出5组具有代表性的试验组别,并为了方便理解进行重新编号(表3),并继续对试验测量数据进行对比分析。

表3 振动试验参数设置(简化)Table 3 Vibration test parameter setting(simplified)

1.2.2 振动前后土体强度变化

利用电子微型贯入仪在每组振动试验前后,分别于斜坡同一水平高度点进行贯入阻力测试。贯入阻力结果可以很好地对振动试验前后进行对比,模拟陆坡土体固结程度变化,将其反映为土体强度的变化。5组试验贯入阻力随深度变化值见图5。

每次试验贯入深度为20 cm,若中途即将超出贯入仪量程则停止贯入,保护压感探头。根据模拟陆坡土层铺设情况,从距斜坡最高点表层20 cm 垂直向下分别为:8 cm 厚上层土体(土样Ⅱ)、2 cm 厚软弱层土体(土样Ⅲ)和10 cm 厚下层土体(土样Ⅰ)。由图5可知,振动后(红色点线)斜坡土体各深度贯入阻力值相较于振动前(黑色点线)总体降低。部分点位情况发生反转,有可能是仪器测量误差导致,或是人为贯入速度每组存在差异而导致。8～10 cm 深处为软弱层土体,当贯入到该层土体后,阻力值有一个明显的持平甚至折减,反映了软弱层土体相较于周围土体强度性质更低,也是本试验引入软弱层这一条件的目的。

图5 贯入阻力随深度变化值Fig.5 Penetration resistance varies with depth

同样利用微型十字板剪切仪,于试验前后分别在斜坡同一水平高度进行剪切试验,反映振动前后土体抗剪强度的变化。5组试验结果见图6。

图6 振动前后土体抗剪强度变化Fig.6 Variation of shear strength of soil before and after vibration

从试验结果可知,振动后各组土体抗剪强度明显降低。相较于前2组试验,当振动参数增大,后2组试验表现为振动后土体抗剪强度的减小更加明显。表明振动强度与土体强度破坏程度呈正相关。

综上所述,地震作用会使得土体强度降低、抗剪强度下降,从而引发滑坡,并且地震强度与土体结构破坏的程度呈正相关关系。

1.2.3 试验滑坡过程

在试验过程中,随着振动的开始,会在水体中扬起“粉尘”,并在模拟陆坡表面形成浑浊流体(图7)。该流体规模不大,但是速度很快,夹带着第一层土体和部分被破坏的软弱结构层,迅速向模拟陆坡前部运动,陆坡前部的高度升高,随后越过挡板,在挡板右侧区域形成大面积浑浊并沉积(图8)。上层土体变形情况不易观察,但通过白色的软弱层可直观地发现,原本相对平整均匀的软弱层变得凹凸不平,厚度极不均匀。说明振动作用打破了原本土体的稳定结构,发生了表层滑动及内部扭曲变形。

图7 陆坡表面浑浊流体Fig.7 Turbid fluid on the continental slope surface

图8 斜坡底部滑坡形态Fig.8 Slope bottom landslide morphology

海底表层土体经常处于弱重力条件下,工程性质极不稳定。所以振动初期,处于弱重力状态下的土体颗粒首先进入上层水体悬浮,水体表面形成浑浊。随着振动进行,土体内部有效应力下降,导致孔隙水压力上升,进而发生土体液化流动,表现为图中大规模的浑浊流体向下运动情况。

试验后,排出所有水体,观察发现陆坡表层存在起伏,有明显的流体流动形成的痕迹。在挡板位置处发现分叉的,形如树杈的沉积,如图9所示。这些沉积都有一个共同的特点,形成的切割较明显且切割得很深,这就表明该浑浊流体在运动时是高速向下的,且内部性质存在差异性,分为高密度区与低密度区,因此会在挡板处形成树杈状的沉积。

图9 挡板处树杈状沉积Fig.9 Tree branch deposit at baffle

2 基于ANSYS强度折减法的模拟陆坡稳定性分析

海底陆坡土体的抗剪强度(τf)满足以下黏性土库伦公式:

式中:c为黏聚力;σ为剪切破坏面上的法向压应力;φ为内摩擦角。

陆坡失稳过程必然伴随着土体抗剪强度的降低,强度折减法的理论即引入一个折减系数Fs,对抗剪强度进行折减并带入有限元软件分析计算,直到软件显示陆坡失稳,此时的Fs视为该陆坡的稳定系数。折减后的黏聚力c′和内摩擦角φ′满足以下公式:

2.1 有限元模型建立

参照上述试验中物理模型尺寸,在ANSYS软件中构建有限元模型。为了达到ANSYS最佳计算效果,将物理模型尺寸按照1∶100等比例放大,并自坡顶、坡脚分别向左、右两侧水平方向延伸,坡脚处垂直向下延伸。修正后的示意图如图10 所示,红框区域为物理模型试验中的模拟陆坡原型,放大后的尺寸长为100 m,高约10 m,坡角约为5°,为主要的陆坡稳定性分析区域。利用ANSYS有限元软件,参照该模型尺寸进行建模,采用Plane 82四边形平面应变单元划分网格,单元尺寸为2 m,共划分为925个单元(图11),并采用Drucker-Prager屈服准则,保证计算连续不中断。

图10 有限元模型尺寸(m)Fig.10 The dimensions of the finite element model(m)

图11 有限元模型网格划分(m)Fig.11 Mesh division of finite element model(m)

根据三轴试验测量并参考海底陆坡土体力学特性参数表[25],确定有限元模型土体力学参数:弹性模量20 MPa、泊松比为0.3、重度为18.5 kn·m-3、黏聚力为16.8 kPa和内摩擦角为21.1°。

2.2 地震波合成与加载

由于水槽试验使用的振动台无法导出加速度时程曲线,故利用振动频率及加速度并根据式(4)进行转换,得到人工合成的间谐波进行加载。

式中:a为加速度;f为频率;a(t)为按每组加速度、频率进行变形后的加速度时程正弦函数曲线。同样选用表2所示的16组频率、加速度组合,具体转换函数见表4,这样便能将ANSYS边坡动力分析的加载条件与物理模型试验中振动台的加载条件对应。

表4 加速度时程曲线公式Table 4 Formula of acceleration curve with time

试验前需先进行模型模态分析,并利用Rayleigh阻尼法确定模型阻尼。最终根据材料特性选择阻尼比为0.05,并通过计算得到质量阻尼系数α=0.267,刚度阻尼系数β=0.009。

2.3 模拟陆坡稳定性分析

分别对有限元模拟陆坡进行静态以及16组地震加载动态折减分析,每组地震加载时间为10 s。以有限元不收敛或模拟陆坡实现塑性应变贯通为判断失稳标准,最终得到各组试验陆坡稳定系数见表5。

表5 各组试验陆坡稳定系数Table 5 Stability coefficient of each group of test slope

如图12所示,模拟陆坡在静力条件下稳定系数可高达4.90,但随着地震加载,稳定系数发生骤减。当加速度达到1.96 m·s-2,频率达到2 Hz时,稳定系数Fs＜1,表明此时陆坡将在地震作用下发生失稳。所以本次有限元模型分析中陆坡失稳的临界地震条件为:振动频率为2 Hz,振动加速度为1.96 m·s-2。

图12 几组试验斜坡失稳表现Fig.12 Slope instability in several groups of tests

选取斜坡表面3个特征点A、B和C(图11),分别位于坡顶、坡中和坡底,提取并分析每组试验3个点位的水平位移、速度、加速度随时间的变化情况。结果如图13所示,三点的位移、速度、加速度值均为点B＞点A＞点C,即坡中最大、坡顶次之、坡底最小。局部点位的加速度值可以明显大过水平加载的地震波加速度值,且在斜坡中部,地震响应最为明显,位移变形程度最大。在试验7的加速度曲线中,4～5 s间发生了明显的“削尖”现象,此现象可大致排除操作失误及外界因素干扰原因,推测为特征点较为靠近斜坡表面从而引起软件计算失误所导致。

图13 三点水平位移、速度、加速度随时间变化曲线(试验7)Fig.13 Curves of horizontal displacement,velocity and acceleration over time at three points(Test 7)

3 对比分析

本次试验有限元模型的模型构建及地震加载基本与物理模型试验相对应,保证了试验结果对比的有效性。综合分析两部分试验结果可知,本次试验构建的模拟陆坡模型在自然条件下处于稳定状态,但在一定强度的地震荷载下会发生失稳,陆坡失稳的临界地震条件为振动频率2～3 Hz、振动加速度1.96 m·s-2。

将2次试验中每组陆坡失稳发生的时间进行对比。选用表4中的5组试验组别,相对应的将物理模型试验中将陆坡表面开始发生破坏的时间作为失稳发生时间,将有限元模型试验中初始状态下有限元计算不收敛时的振动加载时间作为失稳发生时间。结果如图14所示。

图14 两次试验失稳发生时间Fig.14 The occurrence time of instability in two tests

由于试验方式不同,2组试验失稳时间存在跨量级的差距。但整体趋势基本一致,反映为加速度、频率各自变化对斜坡稳定性的影响趋势基本一致。并且对比可发现,相比于频率,振动加速度的增大带来的影响更为显著。通过表6中各组稳定系数对比也可发现这一规律。因此,加速度相比于频率是影响斜坡稳定性更重要的地震参数。

4 结 论

本文采用水槽模型试验方法,依据南海北部陆坡工程环境构建模拟陆坡模型,探讨引发陆坡失稳的临界地震条件等。基于ANSYS有限元强度折减法,分析模拟陆坡在不同地震加载下的稳定系数变化,并对物理模型试验结果进行验证。综合两部分试验结果,得到如下结论:

1)本次试验构建的模拟陆坡,在自然条件下处于稳定状态,随着模拟地震的加载会逐渐发生失稳,失稳发生的临界地震条件为:振动频率2～3 Hz,振动加速度1.96 m·s-2。在这一过程中,加速度的作用占主导地位。

2)振动作用使得原本处于水下弱重力状态的土体有效应力降低,孔隙水压力累积并发生土体液化,滑动力增加,进而导致滑坡发生。试验测得振动过后陆坡土体强度降低,抗剪强度下降,是为陆坡失稳滑动的原因所在。

3)地震波的荷载在斜坡中部响应最为突出,加速度瞬时特征值可超过地震波加速度值,位移变形程度也最大,坡顶次之,坡脚则最小。