基于强度折减法的海底边坡三维稳定性分析方法

摘要

本发明涉及一种基于强度折减法的海底边坡三维稳定性分析方法。该方法包括S1.选定待分析区域；S2.对待分析区域进行前处理；S3.将导出的excel文件，通过ABAQUS数值接口程序生成ABAQUS导入文件.inp；S4.构建自重应力场；S5.进行稳定性分析；S6.对数值计算结果进行分析。本发明通过ABAQUS数值接口程序，使得不仅模型能够导入ABAQUS中，并且能够正确地进行运算并得出相应的研究结果。

说明书

技术领域

本发明涉及一种边坡稳定性分析方法，特别是一种基于强度折减法的海底边坡三维稳定性分析方法。

背景技术

随着我国近海海域工程开发，海底潜在地质灾害尤其是海底边坡失稳日益成为威胁近海工程的重要因素。海底滑坡是指海底斜坡上未固结得到松软沉积物或有软弱结构面的岩石，在外力作用下沿斜坡中的软弱结构面发生滑动的现象。海底边坡稳定性分析方法的研究目前仍然处于发展阶段。随着有限元的产生和发展，应用严格的应力应变方法的分析结构的稳定和变形成为可能，逐渐应用到岩土工程边坡稳定分析当中。

目前大多数的ABAQUS有限元边坡稳定性研究都是基于传统的二维及三维有限元模型得出的结果，已经不能满足研究及工程精度的需要，需要一种可以具体体现海底地貌的有限元建模来完成海底边坡稳定前期的数值计算工作。

GOCAD(Geological Object Computer Aided Design)地质建模软件是法国Nancy大学开发的主要应用于地质领域的三维可视化建模软件。GOCAD软件具有强大的三维建模、可视化、地质解译和分析的功能。既可以进行表面建模，以可以进行实体建模；既可以设计空间几何对象，也可以表现空间属性分布。

ABAQUS是一套功能强大的工程模拟的有限元软件，其解决问题的范围从相对简单的线性分析到许多复杂的非线性问题。ABAQUS包括一个丰富的、可模拟任意几何形状的单元库。ABAQUS除了能解决大量结构(应力/位移))题，还可以模拟其他工程领域的许多问题，例如热传导、质量扩散、热电耦合分析、声学分析、岩土力学分析(流体渗透/应力耦合分析)及压电介质分析。

发明内容

本发明所要解决的技术问题是提供一种能够满足工程精度需要、能实现GOCAD和ABAQUS之间的接口的基于强度折减法的海底边坡三维稳定性分析方法。

本发明解决上述技术问题所采用的技术方案如下：

基于强度折减法的海底边坡三维稳定性分析方法，包括如下步骤：

S1.选定待分析区域；

S2.对待分析区域进行前处理：

S21.将待分析区域的数据导入GOCAD中形成三维地貌图；

S22.将形成的三维地貌图导入Surfer中，并经过插值处理后形成海底表面图；

S23.将海底表面图导入GOCAD中，通过由点-线-面-体的顺序形成网格，组成不规则的上表面条件下的三维海底地形图；

S24.将GOCAD中的三维海底地形图导出为excel文件；

S3.将导出的excel文件，通过ABAQUS数值接口程序生成ABAQUS导入文件.inp，数值接口程序包括如下过程：

S31.将excel文件导入到数值接口程序中；

S32.对excel文件中的生成单元网格各个节点的三维坐标参数，即p0至p7的三维坐标参数，按照下表的方式进行修改：

S33.从数值接口程序中导出为ABAQUS导入文件.inp；

S4.构建自重应力场：

S41.编写导入文件.inp，将边界条件为限制前后两面的x方向位移，左右两面的y方向唯一，底面限制其xy两面的位移；

S42.编写完成后，将导入文件.inp导入ABAQUS中生成模型；

S43.执行ABAQUS中的相应模块，将自重应力场中的单元中心点的六个应力分量F11、F22、F33、F12、F13、F23输出到外部文件.txt并导出；

S44.删除外部文件.txt中数据行以外的文件，并在每个单元的数据前面增加模型名，格式为“模型名-1.序号”，并每一行将同一个单元所属的六个应力分量用逗号分隔；

S5.进行稳定性分析：

S51.将转换后的重力场文件导回原ABAQUS模型之中，在模型顶部生成surface(表面)；

S52.定义分析步长，并对粘聚力和内摩擦角进行折减运算,并输入到输 入ABAQUS程序之中，公式如下：

$C_m=\frac{C}{F_m}---\left(1\right)$

其中，C、分别为粘聚力和内摩擦角，C_m、分别为折减后的粘聚力和内摩擦角，F_m为强度折减系数，数值范围为0.5-10；

S53.在ABAQUS中，分析波浪载荷在t＝0、t＝0.25T、t＝0.5T和t＝0.75T时刻，分别在海床表面施加竖直方向不同时刻的波浪压力，根据微幅波理论，水体的动压力可以由式(3)表示：

$P_z=\rho g\frac{H\cosh \left[k(z+h)\right]*cos(kx-\sigma t)}{2\cos h\left(kh\right)}-\rho gz---\left(3\right)$

当z＝-h时，海底面压力为：

$P_0(x,t)=\frac{\rho gh}{2\cos \left(kh\right)}cos(kx-\sigma t)+\rho gh---\left(4\right)$

其中，z为波面相对高度，H为波高，σ为波动圆频率，ρ为液体密度，k为波数，h为静水面至海底的垂直距离；

S54.利用ABAQUS对模型整体进行计算，输出变量为安全系数F_s与海底面选取点的位移；

S6.对数值计算结果进行分析：

S61.缩减积分单元位移等值线图分析：根据在t＝0、t＝0.25T、t＝0.5T和t＝0.75T四个时刻计算终止时，ABAQUS中分别形成的缩减积分单元位移等值线图，若出现从多个滑动面向单圆滑动面演变，再由单一滑动面向多个滑动面分离的趋势，则进行S62；

S62.通过安全系数判断海底边坡稳定性：

S621.将安全系数F_S与海底面选取点的位移作为输出变量，利用ABAQUS提供的Combine函数，得到不同时刻的安全系数F_S随位移的变化关系图；

S622.通过安全系数F_S随位移的变化关系图得到安全系数值与时刻的对应关系，若安全系数值随时刻变化呈现上下起伏的波浪趋势，取安全系数值的最小值，并利用《建筑边坡工程技术规范》(GB>

作为优选，S21中，待分析区域的数据为用cad软件生成并以.dxf格式保存的待分析区域的海底等深线水深数据。

作为优选，S53中，分析步长至少为一个周期T，周期T为波浪循环，完成一次正弦曲线的时间。

作为优选，S53中，波数k的计算方法为：

S531.在MATLAB环境中，设波数为未知数x，计算公式为：

f＝inline('x-(/(g*tanh(x*h)))')

S532.取f的范围是-100至100；

S533.查找f在-20只20之间的零点；

S534.经过多次迭代后，得到波数k的值。

本发明中的matlab函数：fzero，用于查找一元连续函数的零点；Ezplot，画图函数，它无需数据准备，直接画出函数图形；inline用来定义一个类的内联函数。

本发明中的DXF是Autodesk公司开发的用于AutoCAD与其它软件之间进行CAD数据交换的CAD数据文件格式。DXF是一种开放的矢量数据格式。

本发明同现有技术相比具有以下优点及效果：

1、由于本发明通过ABAQUS数值接口程序，使得不仅模型能够导入ABAQUS中，而且并能够正确的进行运算并得出相应的研究结果。

2、由于本发明通过Visual C#语言平台开发了GOCAD与ABAQUS软件之间的应用接口程序，实现了只需利用CAD等深线图，即可在有限元分析软件ABAQUS中生成复杂的海底地形及地貌模型。

3、由于本发明通过Visual C#语言平台开发了ABAQUS场输出文件转换程序，解决了大数据量条件下重力场数据格式转换，经过该文件转换后的场数据文件完全符合导回ABAQUS软件中进行下一步静力及动力荷载施加计算的要求。

4、由于本发明形成了独特的海底边坡三维稳定性分析方法及技术路线，利用有限元分析手段同时兼顾到了波浪动力荷载作为外部荷载对海底边坡应力状态的影响，提高了计算的精度，简化了海底稳定分析的流程。

附图说明

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

图1为本发明的技术路线示意图。

图2为本发明的模型转换控制技术路线示意图。

图3为本发明的安全系数值随时刻变化示意图。

图4为本发明的单元网格各个节点所在位置示意图。

具体实施方式

下面结合实施例对本发明做进一步的详细说明，以下实施例是对本发明的解释而本发明并不局限于以下实施例。

实施例1：

由于ABAQUS程序模型构造能力不强，尤其是处理复杂的海底表面更是力不从心。GOCAD程序可以满足这一要求。GOCAD可以快速通过CAD生成地形面，如果经过Surfer程序将地形用有限差分法进行克里格插值处理后，三维地貌将更为平滑。GOCAD所生成的网格与ABAQUS所要求网格形状一致(六面体网格)且可用Excel文件导出，这就为接口的实现创造了可能。

基于上述论述，稳定性分析技术路线为：1.将原始的CAD等高线图导入GOCAD中进行前处理，形成基本的等高线点位；2.将生成的等高线导入Surfer中利用克里格插值法进行数据插值，再导回GOCAD程序中形成海底表面；3.把处理后的等高线进行三维地质建模，导出建模网格数据；4.将网格数据导入应用程序接口进行数据处理，并生成新的ABAQUS计算输入文件.inp，导入ABAQUS中进行运算。其技术路线示意图如图1所示。

另外，现有的研究使GOCAD和ABAQUS之间的接口得以实现，但由于其网格解译顺序出现问题，导致可以建出模型，但是不能运算，程序报错。本发明基于Visual C#语言平台，实现只需导入经GOCAD和Surfer两个程序处理后导出的Excel格式的网格数据文件，自动解译数据，即可完成ABAQUS新网格以及附带的边界条件的生成。本发明不仅模型能够导入ABAQUS中，而且并能够正确的进行运算并得出相应的研究结果。该过程的模型转换控制技术路线图如图2所示。

本实施例的基于强度折减法的海底边坡三维稳定性分析方法，包括如下步骤：

S1.选定待分析区域，在本实施例中，待分析区域长1610m，宽890m，波长36.5m，周期T＝5s，波高H＝4m；

S2.对待分析区域进行前处理：

S21.将待分析区域的数据导入GOCAD中形成三维地貌图；

S22.将形成的三维地貌图导入Surfer中，并经过插值处理后形成海底表面图；

S23.将海底表面图导入GOCAD中，通过由点-线-面-体的顺序形成网格组成不规则的上表面条件下的三维海底地形图；

S24.将GOCAD中的三维海底地形图导出为excel文件；

S3.将导出的excel文件，通过ABAQUS数值接口程序生成ABAQUS导入文件.inp，数值接口程序包括如下过程：

S31.将excel文件导入到数值接口程序中；

S32.对excel文件中的生成单元网格各个节点的三维坐标参数，即p0至p7(如图4所示)，按照下表的方式进行修改：

S33.从数值接口程序中导出为ABAQUS导入文件.inp；

S4.构建自重应力场：

S41.编写导入文件.inp，将边界条件为限制前后两面的x方向位移，左右两面的y方向唯一；底面限制其xy两面的位移；

S42.编写完成后，将导入文件.inp导入ABAQUS中生成模型；

S43.执行ABAQUS中的相应模块，将自重应力场中的单元中心点的六个应力分量F11、F22、F33、F12、F13、F23输出到外部文件.txt并导出；

S44.删除外部文件.txt中数据行以外的文件，并在每个单元的数据前面增加模型名，格式为“模型名-1.序号”，并每一行将同一个单元所属的六个应力分量用逗号分隔；

S5.进行稳定性分析：

S51.将转换后的重力场文件导回原ABAQUS模型之中，在模型顶部生成surface；

S52.定义分析步长，并对粘聚力和内摩擦角进行折减运算,并输入到输入ABAQUS程序之中，公式如下：

$C_m=\frac{C}{F_m}---\left(1\right)$

其中，C、分别为粘聚力和内摩擦角，Cm、分别为折减后的粘聚力和内摩擦角，Fm为强度折减系数，数值范围为0.5-10；

S53.在ABAQUS中，分析波浪载荷在t＝0、t＝0.25T、t＝0.5T和t＝0.75T时刻，分别在海床表面施加竖直方向不同时刻的波浪压力，根据微幅波理论，水体的动压力可以由式(3)表示：

$P_z=\rho g\frac{H\cosh \left[k(z+h)\right]*cos(kx-\sigma t)}{2\cos h\left(kh\right)}-\rho gz---\left(3\right)$

当z＝-h时，海底面压力为：

$P_0(x,t)=\frac{\rho gh}{2\cos \left(kh\right)}cos(kx-\sigma t)+\rho gh---\left(4\right)$

其中，z为波面相对高度，H为波高，σ为波动圆频率，ρ为液体密度，k为波数，h为静水面至海底的垂直距离；

S54.利用ABAQUS对模型整体进行计算，输出变量为安全系数F_S与海底面选取点的位移；

S6.对数值计算结果进行分析：

S61.缩减积分单元位移等值线图分析：根据在t＝0、t＝0.25T、t＝0.5T和t＝0.75T四个时刻计算终止时，ABAQUS中分别形成的缩减积分单元位移等值线图，若出现从多个滑动面向单圆滑动面演变，再由单一滑动面向多个滑动面分离的趋势，则进行S62；

S62.通过安全系数判断海底边坡稳定性：

S621.将安全系数F_S与海底面选取点的位移作为输出变量，利用ABAQUS提供的Combine函数，得到不同时刻的安全系数F_S随位移的变化关系图；

S622.通过安全系数F_S随位移的变化关系图得到安全系数值与时刻的对应关系，若安全系数值随时刻变化呈现上下起伏的波浪趋势，取安全系数值的最小值，并利用《建筑边坡工程技术规范》(GB>

上述S21中，待分析区域的数据为用cad软件生成并以.dxf格式保存的待分析区域的海底等深线水深数据。

上述S53中，分析步长至少为一个周期T，周期T为波浪循环，完成一次正 弦曲线的时间。

上述S53中，波数k的计算方法为：

S531.在MATLAB环境中，设波数为未知数x，计算公式为：

f＝inline('x-(/(g*tanh(x*h)))')

S532.取f的范围是-100至100；

S533.查找f在-20至20之间的零点；

S534.经过多次迭代后，得到波数k的值。

波数k的计算方法具体为：

上述S43中，为了将自重应力场的场数据提取出来，需要首先进入Visualization后处理模块，打开相应的计算结果文件。再将自重应力场中的单元中心点的六个应力分量输出到外部文件.txt中导出，然后将场数据导入有静水压力作用下的模型中才能进行进一步的折减分析。但是，导出的场数据不能被ABAQUS/CAE直接读取，必须通过处理后通过ABAQUS/command模块读取。

上述S44中，将外部文件.txt中删除数据行以外的文件并每个单元的数据前面增加模型名，格式为“模型名-1.序号”，并每一行将同一个单元所属的不同应力方向的应力数据，即六个应力分量F11、F12、F13、F12、F13、F23用逗号分隔，如原导出数据为：

通过场文件转换程序将文件转换为：

模型名-1.1,F11，F12，F13，F12，F13，F23

模型名-1.2,F11，F12，F13，F12，F13，F23

模型名-1.3,F11，F12，F13，F12，F13，F23

上述S622中，《建筑边坡工程技术规范》(GB 50330-2013)规定：5.3.1 除校核工况外,边坡稳定性状态分为稳定、基本稳定、欠稳定和不稳定四种状态,可根据边坡稳定性系数按表5.3.1确定。

表5.3.1 边坡稳定性状态划分

注：F_st——边坡稳定安全系数。

5.3.2边坡稳定安全系数Fst应按表5.3.2确定,当边坡稳定性系数小于边坡稳定安全系数时应对边坡进行处理。

表5.3.3 边坡稳定安全系数F_m

注：1.地震工况时，安全系数仅适用于塌滑区内无重要建(构)筑物的边坡；

2.对地质条件很负责或破坏后果极严重的边坡工程，其稳定安全系数应适当提高。

上述F_m为强度折减系数，数值范围从0.5-10之间选取，本实施例中采用的是

此外，需要说明的是，本说明书中所描述的具体实施例，其零、部件的形状、所取名称等可以不同。凡依本发明专利构思所述的构造、特征及原理所做的等效或简单变化，均包括于本发明专利的保护范围内。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，只要不偏离本发明的结构或者超越本权利要求书所定义的范围，均应属于本发明的保护范围。