一种高位远程滑坡后破坏工程防护效果的预测方法与流程

本发明涉及高位滑坡技术领域，尤其涉及一种高位远程滑坡后破坏工程防护效果的预测方法。

背景技术：

目前，受地球内外动力响应作用，我国西部山区区域地质环境复杂，人类工程活动强烈，是特大型高位滑坡的易发区。高位滑坡从高陡斜坡位置剪出凌空坠落，经过势动能转化，具有滑体运动速度极快、冲击能量巨大和远程成灾的灾害特征，容易对周围建筑物造成毁灭性的摧毁，对人类的生命和财产安全造成极大的危害。高效的工程防护手段可以对该类型滑坡的后破坏过程进行人为干预，从而减轻滑坡给人类带来的惨痛伤害。而滑坡的冲击破坏力和工程防护强度之间的平衡一直是地质灾害研究领域定量评估的难点问题，为达到高效经济的工程防护效果，数值仿真的反演和预测定量评估技术显得尤为重要。

现有技术针对高位滑坡后破坏运动过程对防护工程的冲击动力学研究几乎没有，仅从两个方面略有体现：

(1)现有防护技术大多聚焦于滑坡失稳之前的静力学工程防护，采用单纯有限元方法进行锚杆、抗滑桩和挡土墙等技术工程防护；

(2)军工业武器防护技术及装备的应用。

上述现有技术方案缺少防护工程对滑坡后破坏动力学响应的数值仿真研究技术，缺少对于动力学响应的滑坡防护强度和结构设计的定量化支撑。

技术实现要素：

本发明的目的是提供一种高位远程滑坡后破坏工程防护效果的预测方法，该方法能完整且多次反复的获得高位滑坡后破坏的动力学特征，并定量化分析滑坡后破坏过程对防护工程的应力和应变定量的分布情况。

本发明的目的是通过以下技术方案实现的：

一种高位远程滑坡后破坏工程防护效果的预测方法，所述方法包括：

步骤1、根据滑坡数字高程数据建立滑坡三维地形，基于光滑粒子流三维sph3d高保真技术建立滑坡模型；

步骤2、采用有限元fem技术建立防护工程模型；

步骤3、基于所建立的滑坡模型，根据不同的滑坡类型选用不同的边界层摩擦模型和滑坡运动参数；

步骤4、基于所建立的防护工程模型，选取防护工程的结构设计及材料强度参数；

步骤5、通过数值仿真获得滑坡后破坏动力学特征及滑体与防护工程之间的相互作用，并由此得到相应的滑坡速度、冲击能量、防护工程位移及应力云图数值结果图；

步骤6、根据所得到的数值结果图，获得滑坡后破坏运动速度和冲击能量，以及防护工程损毁特征及应力分布情况。

由上述本发明提供的技术方案可以看出，上述方法能完整且多次反复的获得高位滑坡后破坏的动力学特征，可以定量化分析滑坡后破坏过程对防护工程的应力和应变定量的分布情况，为滑坡后破坏动力学和防护工程设计提供研究及试验平台。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。

图1为本发明实施例提供的高位远程滑坡后破坏工程防护效果的预测方法流程示意图；

图2为本发明实施例所建立的滑坡模型的示意图；

图3为本发明实施例所提供的滑坡后破坏运动与结构体相互作用过程中一个时刻的结果示意图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。

本发明实施例所述方法采用光滑粒子流(sph)模拟高位滑坡的后破坏运动过程，采用有限元(fem)技术仿真防护工程模型，在两种技术相结合的情况下进行滑坡冲击防护工程的反演和对比分析，确定高效且经济的滑坡后破坏工程防护手段。下面将结合附图对本发明实施例作进一步地详细描述，如图1所示为本发明实施例提供的高位远程滑坡后破坏工程防护效果的预测方法流程示意图，所述方法包括：

步骤1、根据滑坡数字高程数据建立滑坡三维地形，基于光滑粒子流三维sph3d高保真技术建立滑坡模型；

在该步骤中，首先基于离散节点和网格单元组成的滑坡数字高程数据，建立滑坡三维地形存储数组；

获得包含有滑体位置的数字高程数据，根据滑坡与滑体的基底网格和节点的一致性对应关系，在所建立的滑坡三维地形存储数组中找到相对应的高程数据，并在该高程数据的基础上进行高度叠加，即为滑体在该位置处的最大高度；

再根据滑体数字高程数据确定滑体范围的边界；

对于处于滑体范围内部的节点，根据建模粒度的需要动态的加密或稀疏化生成相应数量的滑体粒子，生成的滑体粒子的初始位置根据其所在单元内的具体位置采用线性插值的方式获得，同时在三维高度方向上，根据建模粒度需要动态的生成相应层数的粒子，粒子的位置同样根据线性插值获得；

再根据赋予的滑体的物质属性，将属性参量作为滑体单元数组的一部分，实现滑坡基底数组和滑体数组的构建，如图2所示为本发明实施例所建立的滑坡模型的示意图。

步骤2、采用有限元fem技术建立防护工程模型；

在该步骤中，防护工程模型的建立采用几何建模软件建立三维几何模型，再导入网格划分软件，进行细致均匀的网格划分操作。

其中，几何建模软件采用solidworks软件完成，该软件是达索公司旗下的solidworks子公司负责开发的三维软件，无论多复杂的几何模型，都可分解成有限数量的几何构型特征，而每一种构成特征，都可用有限的参数完全约束。同时，基于特征的实体模型化系统，研究人员采用基于智能特性的功能去生产模型，如腔、壳、倒角及圆角，可以任意勾画草图，较易改变模型，这一功能特性给工程设计者提供了设计上的简便和灵活。

而对于模型网格生成，采用功能强大的cae应用软件包—hypermesh软件完成，该软件具有强大的网格划分处理功能，它使cae分析工程师可以在高度交互和可视化环境下进行仿真分析工作，与其他的有限元前后处理器比较，hypermesh的图形用户界面容易理解和学习，特别是它支持直接输入已有的三维cad几何模型如proe所建的几何模型，并且导入的效率和生成的模型质量都很高，可以大大减少很多重复性的工作，使得cae工程师能够投入更多精力和时间到分析计算工作上去。

步骤3、基于所建立的滑坡模型，根据不同的滑坡类型选用不同的边界层摩擦模型(流体和颗粒体具有一定差异性)和滑坡运动参数(可以通过实验方法得到)；

在该步骤中，所选取的滑坡运动参数包括两类：

一类为滑坡物质的物性参量，包括物质的密度、粘度、速度、初始温度、初始内能、压力比热容、散体密度、散体粒径、散体碰撞恢复系数；

另一类为滑坡所需采用的边界层摩擦模型，包括层流模型、湍流模型、塑性模型、宾汉姆模型、摩擦模型或voellmy模型，具体来说：

1、层流模型(laminar)：当流速很小时，流体分层流动，互不混合，称为层流。层流模型是典型的光滑流线型流动模型，是牛顿流动定律的延伸，雷诺系数相对较小。碎屑流层流运动中剪切力与流速呈正比关系，而与碎屑流的厚度成反比，具体公式为：

参数：μ为流体动态粘度，v为碎屑流运动速度，h为碎屑流体厚度。

2、湍流模型(turbulent)：当流速很大时，流线不再清楚可辨，层流被破坏，相邻流层间不但有滑动，还有混合，称为湍流，即紊流。湍流基底剪切阻力与碎屑流流速的平方成比例，可采用曼宁公式计算：

参数：ρ为碎屑流密度，v为碎屑流运动速度，h为碎屑流体厚度，n为曼宁粗糙系数，g为重力加速。

3、塑性模型(plastic)：适用于塑性剪切模型，塑性行为发生时，基底剪力强度是恒定的，具体公式为：

τ＝c

参数：c代表一个恒定剪切力，例如粘土不排水剪切力。

4、宾汉姆模型(bingham)：宾汉姆模型是塑性和粘性组合阻力方程，适用于弹性临界屈服强度和粘性强度之间的流体材料，基底阻力通过3次方程求得：

参数：μ为流体宾汉姆粘度，τyield为宾汉姆屈服强度。

5、摩擦模型(frictional)：适用于颗粒材料基底摩擦阻力的计算，基底剪切力与碎屑流有效应力是正比关系，孔隙水压力的存在起着关键作用，具体公式为：

τ＝(σz-u)(z＝b)tanφ＝σz'(z＝b)tanφ

参数：σz为碎屑流正应力，u为孔隙水压力，z为碎屑流厚度，σz'为有效应力。

6、voellmy模型：该模型的运动形式是颗粒摩擦和湍流运动的组合形式。对于一些实例分析中，包括雪崩、岩质崩滑、流动性滑坡、滑坡碎屑流、泥石流等，该流变模型更适用于块体运动的研究分析，具体公式为：

参数：f为碎屑流与基底的摩擦系数(f＝tanφ,φ为摩擦角)，ξ为碎屑流体的湍流系数，该系数表征流速和表面空气作用的影响，ρ为碎屑流密度，v为碎屑流运动速度，σz为碎屑流正应力，g为重力加速，z为碎屑流厚度。

步骤4、基于所建立的防护工程模型，选取防护工程的结构设计及材料强度参数；

在该步骤中，对于防护工程模型中采用的高强度铅和钢材质结构，本构模型采用等提出的修正johnson-cook强度模型，同时防护结构的偏斜采用刚体动力学模型进行描述；

对于防护工程模型中采用水泥、混凝土搭建的结构，本构模型采用考虑高压强、高应变率以及含损伤累积的hjc本构模型进行描述。

具体模型包括：

1)金属材质本构模型

为准确描述防护工程中金属材料受到强冲击作用产生的屈服应力及损伤演化，这里采用等提出的修正johnson-cook强度模型，在该模型中将材料的屈服强度表示为等效应变、等效应变率、损伤变量以及温度的函数，具体表示为：

其中，a,b,c,n,m为材料常数；d为材料的损伤变量，d＝0表示材料未产生损伤，d＝1表示材料完全失效；温度t＝(t-t0)/(tm-t0)为无量纲量，t0为室温，tm为材料的熔点；r为材料的累积损伤塑性应变，r＝(1-d)p，p是累积塑性应变；是用户自定义的参考应变率；材料在出现宏观裂纹时，损伤变量临界值将小于1，失效准则描述为d＝dc≤1，dc为损伤变量临界值。损伤变量d是累积塑性应变p的函数，公式如下：

pd为损伤阈值，pf为断裂塑性应变，与材料的应力三轴度、应变率和温度相关。在johnson和cook提出的剪切损伤演化模型中，pf公式为：

其中，d1-d5为材料常数，σ＝σm/σeq为材料的应力三轴度，σm＝(σx+σy+σz)/3为材料的平均正应力，弹体材料的压强采用gruneisen状态方程进行计算，具体为：

其中，

式中，系数γ＝1.99，cs＝3940，ss＝1.489。

2)混凝土本构模型

混凝土损伤大多引入johnson-cook、taylor-chen-kuszmaul(tck)和riedel-thoma-hiermaier(rht)等本构模型，这些模型均未考虑混凝土的压缩累积损伤，无法精确描述混凝土在超高速冲击下出现裂纹扩展等细节。本实施例将同时考虑了压缩强度的压力相关性、应变率效应和损伤软化效应的hjc本构模型引入到处理防护工程的混凝土在高速冲击作用下的变形损伤问题。

hjc本构模型共包含强度模型、损伤模型和状态方程三部分，其归一化的屈服强度函数为：

其中，σ^*＝σ/f′c，σ为实际计算得到的材料的等效强度；p^*＝p/f′c为材料的归一化压力，其中，p为真实压力；f′c为材料的准静态单轴抗压强度；d为材料损伤变量(0≤d≤1.0)；为材料的等效应变率，为应变率，为材料的参考应变率；a,b,n,c及smax为材料常数。

损伤值d由材料的等效塑性应变和塑性体积应变引起的损伤累积叠加形成，公式为：

式中，δεp及δμp分别为在δt时间内由材料的等效塑性应变和塑性体积应变造成的损伤增量；为恒定压力p下的损伤塑性应变，公式为：

其中，d1、d2为常数。

混凝土在冲击载荷下的压缩变形可以划分弹性阶段、过渡阶段和压实阶段三个阶段，其中：

1)弹性阶段(0≤μ≤μc)

p＝keμ(11)

其中，弹性体积模量ke＝pc/μc，pc和μc分别为采用单轴强度抗压试验所得到的材料压碎压力和压碎体积应变，μ＝ρ/ρ0-1为材料的单元体积应变，ρ和ρ0分别表示单元的实时密度和初始密度。

2)过渡阶段(μc＜μ≤μpl)

在该阶段中，混凝土将发生破碎和断裂，加载的状态方程为：

p＝pc+kp(μ-μc)(12)

式中，塑性体积模量kp＝(pl-pc)/(μpl-μc)；pl为压实应力；压实体积应变μpl＝ρg/ρ0-1，对应于pl处的体积应变；ρg为颗粒密度。

在该阶段，卸载的状态方程为：

p＝[(1-f)ke+f·kp]·μ(13)

插值因子f＝(μmax-μc)/(μpl-μc)，kp和μpl如上述定义，μmax为卸载前混凝土达到的最大体积应变。

3)压实阶段(μ≥μpl)

在该阶段，混凝土完全被压碎，加载时，状态方程为：

其中，为修正的体积应变；k1,k2,k3为混凝土材料常数；卸载时，状态方程为：

p＝k1μ(15)

举例来说，如下表1所示为hjc混凝土强度模型参数列表，如下表2所示为hjc混凝土靶板状态方程参数的列表：

表1hjc混凝土强度模型参数

表2hjc混凝土靶板状态方程参数

步骤5、通过数值仿真获得滑坡后破坏动力学特征及滑体与防护工程之间的相互作用，并由此得到相应的滑坡速度、冲击能量、防护工程位移及应力云图数值结果图；

在所述步骤5数值仿真的过程中，对于滑坡后破坏滑体的运动过程采用sph光滑粒子法进行模拟，对于防护工程受滑体冲击过程采用fem有限元方法与sph粒子法相耦合的方法进行模拟。

其中，sph-fem耦合算法主要包括sph-fem固结算法、sph-fem接触算法、sph-fem转换算法，自主开发了lagrangian框架下的三维sph和fem耦合计算程序——csfp3d。根据求解问题的需要，csfp3d既可以单独使用sph或者fem进行数值计算，也可以进行sph-fem耦合计算，其中sph-fem耦合算法是该程序的关键。加入合适的本构模型和状态方程，可以模拟滑体冲击、结构破坏和流固耦合等使用单一方法计算困难的流体与大变形问题。

具体实现中，上述csfp3d程序是基于linux平台开发的，其中主程序和sph部分是c++编写，fem部分是fortran90编写，sph-fem耦合计算采用c++和fortran混合编程技术。主程序采用面向对象的设计方法，便于进行扩展和移植。csfp3d程序的单元类型、单元积分方式、材料模型、sph光滑核函数、sph临近粒子搜索方式和输出方式如下表3所示：

表3csfp3d程序主要内容

在实际应用中，应按照问题的具体形式以及对计算精度和效率的要求，在不同类型间进行选择。csfp3d程序框架主要有三部分组成：前处理、求解器和后处理。其中前处理依托软件truegrid，通过将数据文件进行格式转换，生成csfp3d程序所需要的输入文件。在满足cfl稳定性条件下，通过sph粒子间距、最小单元尺寸以及材料属性等计算得到sph-fem耦合计算的时间步长。主程序将该时间步长传递给sph计算模块、fem计算模块、sph-fem耦合计算模块分别进行计算，并且在每个时间步内在各个模块之间进行sph粒子和有限单元的数据传输，实现sph和fem的耦合计算。

通过上述数值仿真之后，计算结果后处理依托软件tecplot，按照程序控制信息提供的数据输出方式，输出所有场变量，生成相关动画，并按照程序控制信息提供的粒子/节点编号和变量类型编号，生成相关变量的时间历程曲线，获得相应的滑坡速度、冲击能量、防护工程位移及应力云图数值结果图，如图3所示为本发明实施例所提供的滑坡后破坏运动与结构体相互作用过程中一个时刻的结果示意图。

步骤6、根据所得到的数值结果图，获得滑坡后破坏运动速度和冲击能量，以及防护工程损毁特征及应力分布情况。

值得注意的是，本发明实施例中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

综上所述，本发明实施例所述方法具有如下优点：

1、能完整且多次反复的获得高位滑坡后破坏的动力学特征；

2、可以定量化分析滑坡后破坏过程对防护工程的应力和应变定量的分布情况；

3、可通过对比数值试验比选不同防护设计、抗冲击材料强度和布置位置的优越性，选取最优方案进行设计施工；

4、为滑坡后破坏动力学和防护工程设计提供研究及试验平台。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。