技术领域
[0001] 本
发明涉及计算土
力学、颗粒
流变学与防震减灾技术领域，尤其是涉及一种强震触发的岩土介质渐进破坏、类固-液相变行为的模拟方法。
背景技术
[0002] 强震触发的大型滑坡往往存在分阶段的表现形式：滑坡高速启动前的有限
变形阶段和失稳破坏后的极大变形阶段(以应变幅值的大小为参量)。在有限变形状态下，岩土介质表现为
固体力学特性，将开始进入大变形阶段。而在极大变形状态下，岩土介质发生黏性流动，处于大变形的力学状态，具有不同于固体的本构行为特征。边坡岩土材料从有限变形到极大变形的变化过程，涉及到物理状态由类固态到类液态的转变，灾变演化过程复杂，国际上缺少相关研究。
发明内容
[0003] 本发明的目的就是为了克服上述
现有技术存在的
缺陷而提供一种岩土介质渐进破坏、类固-液相变行为的模拟方法。
[0004] 本发明的目的可以通过以下技术方案来实现：
[0005] 一种岩土介质渐进破坏、类固-液相变行为的模拟方法，包括以下步骤：
[0006] 1)选择滑坡的研究区域，获取研究区域内滑坡体尺寸、峰值与残余强度参数和
密度，以及该区域内有效
地震台站的
地震波时程数据，即地震波
加速度时程，并将研究区域离散为多个粒子；
[0007] 2)根据滑坡体尺寸确定粒子的初始间距Δp；
[0008] 3)根据连续性方程计算粒子的密度ρ，计算式为：
[0009]
[0010] 其中，ρi为目标粒子i的密度，ρj为邻近粒子j的密度，mj为邻近粒子j的
质量，t为时间， 为目标粒子i在α方向上的速度， 为邻近粒子j在α方向上的速度，Wij为邻近粒子j对目标粒子i产生影响的核函数， 为目标粒子i在α方向上的坐标，N为目标粒子i支持域中的粒子总数；
[0011] 4)在SPH方法中引入人工
应力项 以防止相邻粒子的距离太小时出现
张力，其中，斥力项fn的计算公式为：
[0012]
[0013] 其中，n为指数因子， 为目标粒子i的人工应力张量在α和β方向上的分量，为邻近粒子j的人工应力张量在α和β方向上的分量，W(rij)为目标粒子i和邻近粒子j间距离rij的核函数，W(Δp)为初始间距的核函数；
[0014] 5)在动量方程上增加人工应力项，并根据修正后的动量方程获取粒子的速度，具体为：
[0015]
[0016] 其中， 为目标粒子i在α方向上的速度， 为目标粒子i在α和β方向上的应力，为邻近粒子j在α和β方向上的应力，δαβ为狄拉克函数在α和β方向上的分量，Fi为目标粒子i的重力， 为目标粒子i在β方向上的坐标；
[0017] 6)根据滑坡体分阶段的运动表现形式建立统一的本构模型，并计算粒子的应力状态；
[0018] 7)将滑坡模型的底部边界设置为无滑移边界，竖直边界设置为自由滑移边界，获取岩土介质在自重应力下的初始
应力分布；
[0019] 8)将底部边界设置为经过地震
波速度修正的无滑移边界并施加地震波，计算地震作用下的应力分布，则有：
[0020] VAB＝VA-VB＝ξ(VA-Vseismic)
[0021] VB＝(1-ξ)VA+ξVseismic
[0022] 其中，VAB为边坡粒子A相对于边界粒子B的相对速度，ξ为防止数值出现奇异性的修正系数，Vseismic为地震波速度，VB为边界粒子B的速度，VA为边坡粒子A的速度；
[0023] 9)将滑坡模型的密度场、速度场、位移场和应力场通过跳蛙法进行更新；
[0024] 10)判断是否达到预先设置的时间步数，若是，则停止模拟，若否，则继续模拟计算。
[0025] 所述的步骤4)中，指数因子n的取值为2.55，用以保证人工应力的影响范围只局限在最邻近粒子间。
[0026] 所述的步骤6)中，所述的统一本构模型由Drucker-Prager本构模型和
牛顿
流体本构模型组成；
[0027] 计算小变形阶段的固体力学特性时采用Drucker-Prager本构模型，具体为：
[0028]
[0029] 其中， 为应力率张量在α和β方向上的分量， 为偏应变率张量在α和β方向上的分量， 为应变率张量的球体部分，K为体积模量，sαβ为剪应力张量在α和β方向上的分量，G为
剪切模量，Ψ为剪胀
角，为塑性应变率算子，J2为应力张量的第二不变量；
[0030] 计算大变形阶段的
流体力学特性时，采用
牛顿流体本构模型，具体为：
[0031]
[0032] 其中，τ为剪应力，η为
粘度，为剪应变率。
[0033] 所述的步骤8)中，在滑坡模型竖直边界外侧产生自由场边界粒子以防止地震波发生反射，具体为：
[0034]
[0035]
[0036] 其中， 和 分别为边坡粒子A沿x和y轴的应力， 和 分别为边坡粒子A沿x和y轴的速度， 和 分别为边界粒子B沿x、y轴的速度， 和 分别为边界粒子B沿x、y轴的应力，Cp和Cs为地震波的压缩波波速和剪切波波速。
[0037] 所述的步骤9)中，为了缩小计算误差，时间间隔Δt取值满足Courant-Friedrichs-Lewy条件，即：
[0038]
[0039] 其中，Ccour为柯朗系数，h为光滑长度，c为声速。
[0040] 与现有技术相比，本发明具有以下优点：
[0041] 1)本发明首次将SPH方法引入强震触发的岩土介质渐进破坏、类固-液相变行为的数值仿真。改进的SPH方法突破了传统数值方法无法桥接岩土介质滑动与流动不同运动模式等局限性。将地震效应、地质条件、
地貌条件、统一本构模型和分阶段运动控制方程等进行有机结合，捕捉滑动面贯通过程，进而准确预测滑动距离，从而合理评估地震触发滑坡-碎屑流致灾危险性，为工程防震减灾规划与设计提供科学依据，进一步提升我国地震地质灾害防灾减灾
水平。
[0042] 2)针对现有SPH方法只适用于刻画静力条件下流滑灾害的现状，本发明从边界条件对源程序进行改进，使其能再现强震动力条件下岩土介质渐进破坏与类固-液相变行为。地震波以速度的形式施加在底部边界，创建自由场边界防止地震波发生反射。
[0043] 3)本发明创新性地建立统一本构模型来描述岩土介质从小变形到大变形阶段的运动特征。强震触发的大型滑坡往往存在分阶段的表现形式：滑坡高速启动前的有限变形阶段和失稳破坏后的极大变形阶段。本发明通过联合Drucker-Prager本构模型与牛顿流体本构模型构建一个新的统一本构模型，全面描述滑坡启动、高速流动、堆积过程的运动特性，将岩土材料的破坏过程和破坏后流动行为两个环节统一在一个本构方程中。
附图说明
[0044] 图1为唐家山几何尺寸。
[0045] 图2为汶川地震波时程曲线。
[0046] 图3为t＝5s唐家山等效塑性应变
云图。
[0047] 图4为t＝10s唐家山等效塑性应变云图。
[0048] 图5为t＝15s唐家山等效塑性应变云图。
[0049] 图6为t＝16s时唐家山滑坡运动构型。
[0050] 图7为t＝24s时唐家山滑坡运动构型。
[0051] 图8为t＝32s时唐家山滑坡运动构型。
[0052] 图9为t＝40s时唐家山滑坡运动构型。
[0053] 图10为t＝50s时唐家山滑坡堆积构型。
[0054] 图11为滑坡过程中监测点的速度变化。
具体实施方式
[0056] 本发明的目的是提出一种强震触发的岩土介质渐进破坏、类固-液相变行为的数值计算方法，借鉴颗粒流变学，构建考虑类固-液相态转变的计算模型，再现强震作用下滑坡-碎屑流全过程的运动规律，从而合理评估地震触发滑坡-碎屑流致灾危险性。
[0057] 如图12所示，为了实现以上目的，本发明采用的技术方案如下：
[0058] 1)选择滑坡的研究区域，获取研究区域内滑坡体尺寸、峰值与残余强度参数和密度，以及该区域内有效地震台站的地震波时程数据，即地震波加速度时程，并将研究区域离散为多个粒子；
[0059] 2)根据滑坡体尺寸确定合适的粒子初始间距Δp。
[0060] 3)根据连续性方程计算粒子的密度，计算式为
[0061]
[0062] 其中，ρi为目标粒子i的密度，ρj为邻近粒子j的密度，mj为邻近粒子j的质量，t为时间， 为目标粒子i在α方向上的速度， 为邻近粒子j在α方向上的速度，Wij为邻近粒子j对目标粒子i产生影响的核函数， 为目标粒子i在α方向上的坐标，N为目标粒子i支持域中的粒子总数；
[0063] 4)在SPH方法中引入人工应力项 以防止相邻粒子的距离太小时出现张力，其中，斥力项fn的计算公式为：
[0064]
[0065] 式中，n为指数因子，指数因子n的取值为2.55，用以保证人工应力的影响范围只局限在最邻近粒子间。 为目标粒子i的人工应力张量在α和β方向上的分量， 为邻近粒子j的人工应力张量在α和β方向上的分量，W(rij)为目标粒子i和邻近粒子j间距离的核函数，W(Δp)为初始间距的核函数；
[0066] 5)在动量方程上增加人工应力，并根据修正后的动量方程获取粒子的速度，具体为：
[0067]
[0068] 其中， 为目标粒子i在α方向上的速度， 为目标粒子i在α和β方向上的应力，为邻近粒子j在α和β方向上的应力，δαβ为狄拉克函数在α和β方向上的分量，Fi为目标粒子i的重力， 为目标粒子i在β方向上的坐标；
[0069] 6)根据滑坡体分阶段的运动表现形式建立统一的本构模型，并计算粒子的应力状态，统一本构模型由Drucker-Prager本构模型和牛顿流体本构模型组成；
[0070] 计算小变形阶段的固体力学特性时采用Drucker-Prager本构模型，具体为：
[0071]
[0072] 其中， 为应力率张量在α和β方向上的分量， 为偏应变率张量在α和β方向上的分量， 为应变率张量的球体部分，K为体积模量，sαβ为剪应力张量在α和β方向上的分量，G为剪切模量，Ψ为剪胀角，为塑性应变率算子，J2为应力张量的第二不变量；
[0073] 计算大变形阶段的流体力学特性时，采用牛顿流体本构模型，具体为：
[0074]
[0075] 其中，τ为剪应力，η为粘度，为剪应变率；
[0076] 7)将滑坡模型的底部边界设置为无滑移边界，竖直边界设置为自由滑移边界，获取岩土介质在自重应力下的初始应力分布；
[0077] 8)将底部边界设置为经过地震波速度修正的无滑移边界并施加地震波，计算地震作用下的应力分布，则有：
[0078] VAB＝VA-VB＝ξ(VA-Vseismic) (6)
[0079] VB＝(1-ξ)VA+ξVseismic(7)
[0080] 其中，VAB为边坡粒子A相对于边界粒子B的相对速度，ξ为防止数值出现奇异性的修正系数，Vseismic为地震波速度，VB为边界粒子B的速度，VA为边坡粒子A的速度，在滑坡模型竖直边界外侧产生自由场边界粒子以防止地震波发生反射，具体为：
[0081]
[0082]
[0083] 其中， 和 分别为边坡粒子A沿x和y轴的应力， 和 分别为边坡粒子A沿x和y轴的速度， 和 分别为边界粒子B沿x、y轴的速度， 和 分别为边界粒子B沿x、y轴的应力，Cp和Cs为地震波的压缩波波速和剪切波波速；
[0084] 9)将滑坡模型的密度场、速度场、位移场和应力场通过跳蛙法进行更新，为了缩小计算误差，时间间隔Δt取值满足Courant-Friedrichs-Lewy条件，即：
[0085]
[0086] 其中，Ccour是柯朗系数，h代表光滑长度，c代表声速；
[0087] 10)判断是否达到预先设置的时间步数，若是，则停止模拟，若否，则继续模拟计算。
[0089] 以汶川地震触发的唐家山滑坡启动与崩滑全过程为例
[0090] 地震前，唐家山斜坡尺寸如图1所示，坡度约为39°，经过长时间的地质
风化作用，斜坡分为强风化岩与弱风化岩部分。2008年5月12日，唐家山在汶川地震作用下发生失稳破坏，滑坡体物质(主要为强风化岩)解体后以较大的速度沿地形运动，形成较大体积方量的碎屑流，导致84人遇难。
[0091] 采用本发明的数值仿真方法，对汶川地震触发的唐家山滑坡启动与崩滑全过程灾害进行仿真模拟。具体过程如下：
[0092] 1)汶川地震前，唐家山边坡的研究区域选择为水平距离1712米，高度为872米，其尺寸如图1，计算参数如表1所示。选择研究区域内的有效地震动台站，获取该台站的地震波时程数据，如图2。
[0093] 表1唐家山计算参数
[0094]
[0095] 2)将研究区域离散成14752个粒子，其中边界粒子762个，粒子初始间距为8米，粒子的初始速度为0，其他计算参数按照表2设置。
[0096] 表2 SPH模拟中所用参数
[0097]
[0098] 3)根据连续性方程计算粒子的密度ρ。
[0099] 4)根据公式(2)计算人工应力项，防止相邻粒子的距离太小时出现张力。
[0100] 5)在动量方程上增加人工应力项，根据修正后的动量方程(3)计算粒子的速度。
[0101] 6)针对滑坡分阶段的运动表现形式，联合Drucker-Prager本构模型和牛顿流体本构模型，建立统一的本构模型，根据公式(4)和(5)，计算粒子的应力状态。
[0102] 7)将滑坡模型的底部边界设置为无滑移边界，竖直边界设置为自由滑移边界，获取唐家山边坡在自重应力下的初始应力分布。
[0103] 8)将底部边界设置为经过地震波速度修正的无滑移边界并施加地震波，根据公式(6)和(7)计算地震作用下的应力分布。
[0104] 在边坡模型竖直边界外侧产生自由场边界粒子以防止地震波发生反射，参见公式(8)和(9)。
[0105] 9)将滑坡模型的密度场、速度场、位移场、应力场等通过跳蛙法(Leapfrog Method)进行更新。为了缩小计算误差，根据公式(10)，设置时间间隔Δt。
[0106] 10)判断有无达到预先设置的时间步数，若有，则停止计算，若没有，则继续计算。等效塑性变形计算结果如图3-5，可见，滑移面在第15s时贯通。随后，滑坡启动，图6-10展示了滑坡体从第16s到第50s时的崩滑运动、堆积过程的构型。滑坡从启动到崩滑全过程中三个监测点的速度时程见图11。
[0107] 本发明借鉴颗粒流变学的思想，联合Drucker-Prager本构模型与牛顿流体本构模型构建了一个新的统一本构模型，实现地震触发的岩土介质从小变形到大变形阶段的跨相态描述，重现滑坡-碎屑流全过程的运动规律，从而准确评估地震触发滑坡-碎屑流致灾危险性，为工程防震减灾规划与设计提供科学依据，进一步提升我国地震地质灾害防灾减灾水平。
