基于Abaqus的边坡可靠度计算

基于Abaqus的边坡可靠度计算

一、利用Abaqus建立边坡模型并进行强度折减计算

Abaqus作为强大的有限元通用软件，可以对许多岩土工程领域进行数值模拟。当然，边坡工程也是其中之一。
我国自然灾害频发，传统的极限平衡法不能得到边坡工程的应力应变数据，而有限元方法可以将土体的本构模型考虑在内，求出土体的应力应变并通过Abaqus强大的后处理功能绘制出云图。故有限元强度折减法近年来被越来越多技术人员熟知，对于较为复杂的边坡模型更加适用。
Abaqus软件分为前处理—计算—后处理三大模块。前处理是指利用Abaqus模仿实际物体的参数、形状、约束等因素。计算是指对模型进行有限元计算。后处理是指对有限元计算结果进行分析、处理、出图等操作。
本章将详细讲述这三部分过程，为可靠度计算提供模型基础。

1.在CAE中建立边坡模型

Abaqus早期是通过命令建立模型，并没有友好的用户交互界面，后期为了方便更多人使用，开发出了CAE模块，win10系统通过系统搜索“cae”即可调出。

当然，也可以试着用命令行控制模型建立，此种方法虽然门槛较高，但是更加容易帮助我们了解有限元软件的运行原理。同样在win搜索中输入“abaqus command”即可

CAE(Computer Aided Engineering)是指计算机辅助工程，顾名思义就是通过计算机对工程问题进行模拟，以达到预测、复现的效果，这样可以减少预算，避免损失。

进入cae界面后，首先要进行模型轮廓的建立，本文使用简单的二维模型进行强度折减计算，并且在后续的章节使用通用的模型进行分析。

Create part—更改模型名称—选择"2D Planar"—Continue

通过以上操作确定模型为二维Deformable模型

用Create lines画出边坡轮廓线（暂时不管尺寸）

用Add Dimension对每条边进行赋值

赋值完成后如图所示

通过上述操作后点击鼠标中键，完成轮廓绘制

至此便完成了边坡轮廓的建立
此算例为费康老师《Abaqus在岩土工程中的应用》强度折减法算例

2.将CAD模型导入CAE

当然，有时我们需要把已有的较为复杂的CAD图导入CAE中，具体操作参考我的另一篇文章
CAD图形导入Abaqus2020方法

3.建模过程

边坡有限元强度折减法建模需要用到七个模块(Module)分别是Part-Property-Assembly-Step-Load-Mesh-Job,接下来按这个顺序操作

1. Part
   在绘制完边坡轮廓后，要建立点集，在后面关键字修改时会用到
   点集是abaqus是十分重要的一个概念，在软件学习过程中也需要养成建立点集的好习惯，方便对不同区域进行操作。
   操作流程：Tools—Set—Create—Continue—选择边坡整体，点击Done或者鼠标中键

2. Property
   Property分为创建材料—创建截面属性—赋予截面属性，下面按照这个顺序进行操作
   创建材料：Create Material
   需要对以下几个材料属性进行赋值
   密度：General—Density
   Elastic（弹性）：Mechanical—Elasticity—Elastic
   Mohr Coulomb Plasticity(摩尔库伦属性)：Mechanical—Plasticity—Mohr Coulomb Plasticity（分为粘聚力和内摩擦角）
   下图是对这些参数赋值后的截图
   
   从左下和右下两张图中可以看出，Abaqus是通过场变量控制折减系数的，关于折减系数的生成，可以利用excel完成
   
   从图中可知，除了第7行中没有公式，其他行将红色方框中的数改为Di(i为行)，之后通过改变A7和B7中的数值即可完成自动折减，当然，此过程还可以通过matlab或Python完成，今后会不断完善。
   
   创建截面属性：Create Section—Solid—Homogeneous—Continue—OK
   赋予截面属性：Assign Section—选择边坡轮廓—点击鼠标中键—OK
   

3. Assembly
   任何一个结构都可以视为一个（Instance），它由很多个部件（Part）构成，使用Assembly（装配）功能可以把多个部件组装起来，创建为一个实体（Instance），并在整个坐标系中为这些实体定位，形成一个完整的装配件。即使只有一个Part，也需要进行装配，可以把部件理解为流水线上的配件，Assembly为加工的机器，经过加工后的配件才能进行计算。
   操作流程：Create Instance—OK

4. Step
   分析步分为"Load"和"Reduce"两步，代表加载和折减，这两步的名称按自己习惯定义，但是之后提取ODB文件时，需要与分析步名称一致才可成功提取。
   操作流程1：Create Step— 修改分析步名称（“Load”）—Initial(默认)—Static,General(默认)—Continue—Other—勾选Unsymmetric选项（摩尔库伦失效模型必须使用非对称存储方式）—OK
   操作流程2：Create Step— 修改分析步名称（“Reduce”）—Initial(默认)—Static,General(默认)—Continue—Other—勾选Unsymmetric选项（摩尔库伦失效模型必须使用非对称存储方式）—OK
   控制场变量FV的输出：Field Output Requests Manager—双击第一行—选择State/Field/User/Time中的FV—OK

5. Load
   对边坡模型施加重力，有两种方式 1.施加重力（需要定义密度） 2.施加场力 两种方法均可，下面以施加重力为例说明
   操作流程：Create Load—Step（选为刚才定义的第二个分析步）—Category（点选“Mechanical”）—Types for Selected Step（点选Gravity）—Continue—点击Region后面的小箭头，扩选所需加载重力的区域（全部扩选）—点击鼠标中间确认—Component1（填“0”）—Component（填“-9.8”）表示Y方向的负方向施加9.8的重力加速度—OK
   成功施加重力后应显示下图界面
   
   对边坡进行边界条件控制
   操作流程：Create Boundary Condition—Step（选择第一个分析步）—Category（Mechanical）—Types for Selected Step—Diplacement/Rotation（控制位移和旋转）—Continue—同时按住键盘上的"Shift"点选左右两个边界—点击鼠标中键—选择U1且位移为0（表示仅有水平方向约束）—OK—再以同样的方式约束下边界，但注意下边界约束U1,U2方向均为0
   
   至此就完成了模型load模块的操作

6. Mesh
   对模型网格进行划分
   操作流程：首先点选Object中的Part见下图，表示对part进行网格操作
   
   区域划分
   **操作流程：Partition Face Sketch：将模型划分为规则的三块（方便单元划分） **
   

网格类型
操作流程：点击Element Type—扩选模型—下拉Family框，选择—Plane Strain（平面应变问题)—OK
网格种类
操作流程：点击Mesh Controls—Element Shape（选择“Quad”）—Technique（“Structure’”）—OK
布种
操作流程：点击Seed Part—默认点击Apply—如果种子较疏松，则通过改变Approximate global size：（改为较小的值）—OK
创建网格
操作流程：点击Mesh Part—选择全部区域—点击鼠标中键

创建完网格后要修改关键字，将关键字插入，具体操作见动态图

第一个关键字
*initial conditions, type=field, variable=1
part-1-1. set-1, 0.5
其中set-1为一开始创建的点集，表示强度折减的区域为整个边坡，0.5表示强度折减从0.5开始进行
第二个关键字
*field, VARIABLE=1
part-1-1. set-1,3
其中3表示强度折减最多折减到3
注意“set-1”必须与之前整个边坡的点集名称一致！！！
7. Job
操作流程：Create Job—Continue—OK—Job Manager—Submit
通过以上7步就完成了边坡强度折减法的建模和计算过程

4.后处理

计算完成后可以任意的将计算结果进行后处理分析，比如绘制出U1-FV1曲线，观察位移和折减系数的变化关系，通过上述建模后并利用费康老师书中的绘制方法会弹出错误信息
"X Value is not monotonic and interpolation is undefined"
经过尝试后发现，是因为我修改了折减系数的步长，如果将材料中内摩擦角和粘聚力按照费康老师书中设置，就可以绘制出曲线。

二、Abaqus强度折减法INP文件解读

在了解INP文件之前，首先要注意，INP并不是一种语言，而是按照Abaqus书写要求组织成的一个文档，其中包括模型的全部信息（部分子程序没有），在CAE中对模型进行计算，等同于调用Abaqus引擎对INP文件进行计算。所以，如果想要传输模型信息，INP当然是首选。

本文针对强度折减法INP文件进行讲解，避免陈述太多与强度折减法无关的讲解。

首先INP文件包括模型数据和历史数据两部分组成：1.模型数据 2.历史数据
1.模型数据的作用:定义一个有限元模型.包括单元,节点,单元性质,定义材料等等有关说明模型自身的数据.模型数据可被组织到零件中(零件可以被组装成一个有意义的模型)。
2.历史数据的定义是模型发生了什么----事情的进展,模型响应的荷载,历史被分成一系列的时步层序.每一步就是一个响应。
不论哪种数据，他们都由关键字行、数据行和注释行组成。关键字行一般用*开头，数据行必须紧跟关键字行，注释行以**开头

1.关键字

强度折减INP文件一般有下面这些关键字行
*Heading ：写作规范，INP文件必须以此为开头
*Part ：定义部件名称
*Node：节点信息，包括所有节点的编号和坐标
*Element：单元信息，包括所有单元的编号、单元类型，单元对应节点编号
*Nset：节点点集
*Elset：单元点集
*Solid Section：实体截面定义
*End Part：结束部件部分
*Assembly：装配
*Instance
*End Instance
*Nset：装配后的点集名称
*Elset：装配后的单元铭恒
*End Assembly：结束装配
*Material：定义材料
*Density：密度
*Mohr Coulomb：定义内摩擦角
*Mohr Coulomb Hardening：定义粘聚力
*Boundary：定义边界条件
*Step：定义分析步1，一开始为load分析步
*Static：
*Dload：定义重力加速度
*Output, field：定义场输出变量
*Node Output：定义节点输出
*Element Output,：定义单元输出
*End Step：结束分析步1
*Step：分析步2，一般为reduce分析步
*Static
*Node Output：节点输出
*Element Output：单元输出
*End Step：结束分析步2

2.数据

1.一个数据行包括空格在内不能超过256个字符。

2.所有的数据条目之间必须用“，”格开。

3.一行中必须包括指定说明的数据条目的数字。

4.每行结尾的空数据域可以省略。

5.浮点数最多可以占用20个字符。

6.整数可以是10个

7.字符串可以是80个

8.延续行可以被用到特定的情况。

3.强度折减法关键字

Abaqus是通过控制关键字实现强度折减法的，关键字中包含场变量的起始值和终止值，还包括参与强度折减单元的点集，这在之前建模时都有所提及，下面再细致的讲解一下

第一个关键字
*initial conditions, type=field, variable=1
part-1-1. set-1, 0.5
其中part1-1.set-1就是模型在装配后的点集编号，可以从cae界面左侧模型树中查询，Model-1—Assembly—Sets，只要关键字中的点集和装配后总体模型的点集名称一致，则关键字设置正确。“0.5”是指场变量的起始值。

** 第二个关键字**
field, VARIABLE=1
*part-1-1. set-1,3
其中"3"是指强度折减的终止值

从关键字可知，强度折减法通过变换场变量的值，将内摩擦角和粘聚力进行折减，首先从0.5折减，也就是将强度乘以二倍，如果边坡计算收敛，则继续增加折减系数，通过不断的折减最终以计算不收敛为结束准则。

三、抽样方法

由于可靠度的计算和蒙特卡罗方法密不可分，而蒙特卡洛法又和抽样密不可分，故本章会大致介绍一下抽样方法的相关知识。
蒙特卡洛法：蒙特卡洛(Monte Carlo)方法首先生成随机变量的样本，然后将随机变量的样本作为输入获得功能函数的样本，再统计是小样本的数量从而估算失效概率。

1.随机数的生成

在说明抽样方法前，首先要介绍随机数的产生，因为想要得到随机抽样样本的前提是可以生成较为“随机”的随机数，比方说，用什么方法才能从0-1之间抽取一个数，怎样保证抽样过程是真的随机？
随机数又分为真随机数和伪随机数，真随机数要求产生的数列之间没有相关关系， 在给定初值的情况下数列不具有重现性，但是伪随机数则恰好相反，伪随机数在给定初值的情况下会生成具有重复性的随机数，这就要求计算机在生成随机数时应该具有足够长的周期，使得在数列重复之前能产生足够多的随机数。
随机数一般都从0-1均匀分布随机数的生辰开始，然后通过逆变换法，将0-1均匀分布转化为其他分布类型。
具体相关知识可以参考祝玉学老师的《边坡可靠性分析》一书和张璐璐老师的《岩土工程可靠度理论》一书。在此不加赘述。

2.拉丁超立方抽样

常用的抽样方法有重要抽样方法和拉丁超立方抽样方法。本节主要讲解拉丁超立方抽样方法。具体内容可以从参考我的另一篇文章图解拉丁超立方抽样。抽样方法之所以被研究，是因为他们相较普通简单的抽样方法可以在达到相同计算精度的前提下抽取较少的样本，提高计算的效率。可以理解为，通过各种抽样方法抽取的样本能够更好的代表某个随机变量的统计特征。
拉丁超立方抽样可以在给定抽样函数的情况下，采取更有策略的抽样方法，获得更能代表抽样函数的随机数，从而统计减小误差。

四、参数相关性

上一章对抽样方法进行了讨论，但是通常情况下，功能函数包含多种随机变量，也就是受到许多不同因素的影响，而这些因素同时又具有相关性。比如边坡的安全系数可以受到粘聚力和内摩擦角的影响，但是实验结果表明，内摩擦角和粘聚力通常呈现负相关关系，也就是土体的内摩擦角大则粘聚力就有可能较小。这种相关关系对计算结果也有一定的影响，应该加以讨论。

1.相关系数矩阵

首先介绍相关系数的概念：相关系数是用来形容两组数据之间相关性的度量，相关系数取值区间在[-1,1]之间，如果相关系数大于零，则称两组数据正相关，反之称两者负相关，当相关系数等于0时，则称两组数据线性无关，注意这里所说的是线性无关，但有可能非线性相关，也就是拥有非线性的相关关系。假设有两个随机变量X、Y，通过抽样后得到了他们的样本向量X=(x1,x2,x3,…,xn)共n个元素，同理Y=(y1,y2,y3,…,yn)共n个元素。通过以下公式可以得到两组数据的相关系数。

现在假设有四个随机变量则，四个随机变量之中，每两个随机变量就可以通过上式计算出他们的相关系数，则四个随机变量可生成相关系数矩阵，如下图：

其中，随机变量和其本身的相关系数等于1，且随机变量1和2的相关系数等于随机变量2和随机变量1的相关系数相同。故相关系数矩阵是对角线为1的对称矩阵。
综上可知，相关系数矩阵就是用来形容一些随机变量之间相关性的矩阵。

计算相关系数的相关代码（python代码）如下：

import numpy
Y=[1,2,3,4,5,5]
X=[5,4,3,2,1,0]
t=numpy.corrcoef(X,Y)
print(t)

结果

ans =
[[ 1.         -0.98198051]
 [-0.98198051  1.        ]]

2.相关系数矩阵的乔列斯基分解(Cholesky)

什么是乔利厄斯基分解？
Cholesky 分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解。它要求矩阵的所有特征值必须大于零，故分解的下三角的对角元也是大于零的。Cholesky分解法又称平