文章目录

- 基于Abaqus的边坡可靠度计算
- - 一、利用Abaqus建立边坡模型并进行强度折减计算
  - - 1.在CAE中建立边坡模型
    - 2.将CAD模型导入CAE
    - 3.建模过程
    - 4.后处理
  - 二、Abaqus强度折减法INP文件解读
  - - 1.关键字
    - 2.数据
    - 3.强度折减法关键字
  - 三、抽样方法
  - - 1.随机数的生成
    - 2.拉丁超立方抽样
  - 四、参数相关性
  - - 1.相关系数矩阵
    - 2.相关系数矩阵的乔列斯基分解(Cholesky)
  - 五、通过Python批量生成INP文件
  - - 1.Python字符串处理基础知识
    - 2.关于强度折减法土体随机参数的生成
    - 3.服从对数正态分布参数的生成
    - 5.利用for循环批量生成INP文档
  - 六、批量计算INP文件
  - 七、批量提取ODB文件并进行可靠度计算
  - - 1.利用py文件提取odb中的场变量
    - 2.批量提取ODB结果文档
  - 八、可靠度结果后处理
  - - 1.CDF、PDF曲线绘制
    - 2.失效概率计算
后记

基于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,接下来按这个顺序操作

Part
在绘制完边坡轮廓后，要建立点集，在后面关键字修改时会用到
点集是abaqus是十分重要的一个概念，在软件学习过程中也需要养成建立点集的好习惯，方便对不同区域进行操作。
操作流程：Tools—Set—Create—Continue—选择边坡整体，点击Done或者鼠标中键
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
Assembly
任何一个结构都可以视为一个（Instance），它由很多个部件（Part）构成，使用Assembly（装配）功能可以把多个部件组装起来，创建为一个实体（Instance），并在整个坐标系中为这些实体定位，形成一个完整的装配件。即使只有一个Part，也需要进行装配，可以把部件理解为流水线上的配件，Assembly为加工的机器，经过加工后的配件才能进行计算。
操作流程：Create Instance—OK
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
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模块的操作
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分解法又称平方根法，是当A为实对称正定矩阵时，LU三角分解法的变形。(百度百科)
如图：L即为下三阶矩阵，LT是L的转置矩阵。
代码如下：

import numpy as np
co_mat1=[[1,0.5],[0.5,1]]
down_mat1=np.linalg.cholesky(co_mat1)
print(down_mat1)

结果：

[[1.        0.       ][0.5       0.8660254]]

关于为什么进行乔利厄斯基分解？
在抽样分析中，乔列斯基分解主要用来将两组独立样本转化为具有一定相关性的样本。

如何应用乔利厄斯基分解
将经过乔列斯基分解后的下三角矩阵左乘样本点矩阵就可以得到相关的样本点矩阵
假设[X,Y]=[[x1,x2,x3,x4,…,xn],[y1,y2,y3,y4,…,yn]],则样本矩阵为2行n列的矩阵，L为下三角矩阵，其形状为两行两列，用L左乘样本矩阵L[X,Y]=[X1,Y1]依然得到一个新的两行n列的矩阵。此矩阵中的参数具有一定的相关性。如图
此图为相关系数等于0.9时变换前后的散点图，蓝色点为独立比标准正态分布的样本点，红色点为经过相关性处理的样本点。

以matlab代码为例：

nsamples=input('请输入试验次数');
xr=[1,0.9;0.9,1] ;      %xr为相关系数矩阵
L=chol(xr,"lower");     %获得下三角矩阵
x=LHSNoCorr1(nsamples); %利用拉丁超立方抽样获得的样本点
y=LHSNoCorr1(nsamples); %利用拉丁超立方抽样获得的样本点
figure(1)
plot(x,y,".b")          %蓝色点为独立数据散点图
A=[x;y];                %获得原始样本点矩阵
B=L*A;                  %左乘下三角矩阵
x1=B(1,:);              %获得相关变化后的样本点矩阵
y1=B(2,:);              %获得相关变化后的样本点矩阵
hold on
plot(x1,y1,".r")        %红色点为相关数据散点图%拉丁超立方抽样函数，可扩展维度
function newz=LHSNoCorr1(M)
z=[];
for i=1:Mz(1,i)=norminv(rand(1,1)/M+(i-1)/M);  %rand(1,1)表示生成一行一列个01均匀分布随机数
end
[temp1,O1]=sort(rand(M,1));    %O1表示一个随机序列，他是对随机生成的M个正态分布进行排序后的数列
newz(1,:)=z(1,O1)';            %利用刚才生成的随机序列来对z进行乱序排列，这样去除了数据之间的相关性
end

五、通过Python批量生成INP文件

通过对INP文件格式的简单理解，可以看出，INP文件将有限元模型的不同步骤用关键字隔开定义，但是进行可靠度计算时，主要是对土体的抗剪强度参数进行变换，通过对土体参数赋新值，生成新的模型进行有限元计算。故只需要对INP文件中的材料部分进行批量修改即可，通过这种方式，批量生成INP计算文件，进行边坡可靠度计算。

1.Python字符串处理基础知识

批量生成INP计算文件包括几个步骤

把INP文件中始终不变的内容以字符串的形式提取到缓存中，比如边界条件、单元编号、单元节点对应关系、节点坐标等信息。
利用python生成新的材料属性并转化为字符串。
将前两步生成的字符串按照INP文件的书写格式合并起来，并存文txt文档，将txt文档修改为inp文档。
通过以上三步就生成了一个INP文件，再通过循环生成任意数量的INP文件。

从上面的叙述可知，对于文档的提取和重写是批量生成INP文件的重点，所以本节主要介绍python处理字符串的相关技巧。所用到的python库主要有 numpy、random和scipy库。

①读取INP文件
在建立一次边坡模型后，把INP文件放在Python的工作目录下，
首先建立一个INP文件，内容如下：

利用with open语句(上下文管理器)打开INP文件，并利用readlines语句将INP文件中的每一行作为一个元素放在列表中。

with open("test5-1.inp",'r',encoding='ANSI') as txtfile:line=txtfile.readlines()print(line)

结果：

['11111111111111\n', '22222222222222\n', '33333333333333']

可见通过以上操作，将INP文件中的每一行字符串都进行读取并放在了列表中，每行字符串占据列表的一个元素。

2.提取相关信息
将INP文件读取完毕后，需要提取出我们需要的内容，要用到字符串的一系列方法。本小节用到的方法主要有：enumerate find strip 和列表的append方法，enumerate是用来遍历列表内容的，find是用来查找指定字符串，strip是用来处理字符无用字符的，append用来生成新的列表。

以下是相关代码

with open("test5-1.inp",'r',encoding='ANSI') as txtfile:line=txtfile.readlines()list2=[]list4=[]for i,rows in enumerate(line):if rows.find("*End Assembly")!=-1: print(f"第{i+1}行出现了'*End Assembly'")line_1=i+1elif rows.find(" 7.02097,0.")!=-1:print(f"第{i+1}行出现了' 7.02097,0.'")line_2=i+1elif rows.find("*End Step")!=-1:print(f"第{i+1}行出现了'*End Step'")line_3=i+3else: #print(f"第{i+1}行没有出现")continuefor i,rows in enumerate(line):if i in range(2,line_1):list2.append(rows.strip('\n'))elif i in range(line_2,line_3):list4.append(rows.strip('\n'))
print(list2)
print(list4)

计算结果如下：

第5行出现了'*End Assembly'
第10行出现了' 7.02097,0.'
第15行出现了'*End Step'
['22222222222222', '33333333333333', '*End Assembly']
['11111111111111', '22222222222222', '33333333333333', '------------------', '*End Step']

从结果可知,最终将文档中 *End Assembly行之前的所有行存在列表list2中，将数据7.02097,0.之后的内容放在list4中，这两个列表中的内容将在后续使用。并且将特定字符串所在的行数提取出来。
文档内容如下：


11111111111111
22222222222222
33333333333333
*End Assembly
11111111111111
22222222222222
33333333333333
------------------7.02097,0.
11111111111111
22222222222222
33333333333333
------------------
*End Step

将不变的部分提取出来以后，下一步就是不断生成新的部分，将两部分信息合并，输出为INP文件。当然，新的部分在这里也就是指材料的参数部分。
代码如下：

str(np.round(np.arctan((np.tan(chi1_phi[i]*(np.pi/180))/1))*(180/np.pi),10))+','+'  '+str(0.1)+','+' '+','+' '+'1'+'\n'

通过复制粘贴修改此语句，将0.4-3的折减系数都变为字符串形式，通过改变chi1_phi的值进行材料内摩擦角的给定。粘聚力操作方法类似。

2.关于强度折减法土体随机参数的生成

Abaqus边坡稳定性计算方法为有限元强度折减法，需要不断变化的土体抗剪强度参数值。本部分内容也主要用到字符串相关语法。
参数的生成流程如下：
①.首先利用拉丁超立方抽样方法抽取独立标准正态分布样本点
②.将独立标准正态分布样本点转化为相关正态分布样本点
③.将上面的到的样本点转化为服从一定均值和标准差的实际样本点

3.服从对数正态分布参数的生成

首先生成服从正态分布的样本点

##生成土层一的土体参数，参数的个数和随机响应面的配点数一致
co_1_norm=0.5
c1=LHS(Pei_number)
phi1=LHS(Pei_number)
tuceng1=np.vstack((c1,phi1))           #将两组样本点组合在一起
co_mat1=[[1,co_1_norm],[co_1_norm,1]]  #生成相关系数矩阵
down_mat1=np.linalg.cholesky(co_mat1)  #对相关系数矩阵进行乔列斯基分解
#利用乔列斯基分解法将独立标准正态分布转化为相关标准正态分布
chi1_1=np.dot(down_mat1,tuceng1)
chi1_c=c1_sigma_norm*chi1_1[0]+c1_mean_norm
chi1_phi=phi1_sigma_norm*chi1_1[1]+phi1_mean_norm

由于抗剪强度参数一般为正值，故对数正态分布更适合作为其分布类型。如果想要生成服从对数正态分布的样本点，就必须知道正态分布和对数正态分布的关系。
一般来说，假如X服从正态分布，那么e^X就服从对数正态分布，也就可以得到对数正态分布与正态分布统计值之间的关系如下：

其中σ_xi和μ_xi为对数正态分布的标准差与均值，σ_lnxi和μ_lnxi为正态分布的标准差与均值。
先将对数正态分布的统计参数（均值和标准差）转化为正态分布的统计参数，再通过此均值与标准差生成一系列正态分布样本点X，再考虑参数之间的互相关性，之后用e^X将正态分布转化为对数正态分布。
上述过程中需要将对数正态分布的相关系数先转化为正态分布的相关系数，用到了Nataf变化的相关知识，在此不赘述，具体可以参考李典庆老师的《边坡可靠度非侵入式随机分析方法》，对数正态分布与正态分布相关系数转换关系如下：

ρ_0ij为正态分布互相关系数。
具体代码如下：

##生成土层一的土体参数，参数的个数和随机响应面的配点数一致
c1=LHS(Pei_number)
phi1=LHS(Pei_number)
tuceng1=np.vstack((c1,phi1))
co_1_norm=np.log(co_1*cov_c1_ln*cov_phi1_ln+1)/(np.sqrt(np.log(1+cov_c1_ln**2)*np.log(1+cov_phi1_ln**2)))
co_mat1=[[1,co_1_norm],[co_1_norm,1]]
down_mat1=np.linalg.cholesky(co_mat1)
#利用乔列斯基分解法将独立标准正态分布转化为相关标准正态分布
chi1_1=np.dot(down_mat1,tuceng1)
#土层一的粘聚力值(正态)
c1_mean_norm,c1_sigma_norm=lnnorm_norm(c1_mean_ln,c1_sigma_ln)
#土层一的内摩擦角值(正态)
phi1_mean_norm,phi1_sigma_norm=lnnorm_norm(phi1_mean_ln,phi1_sigma_ln)
#对(相关/不相关)正态分布样本点取对数就得到了对数正态分布的样本点，这里是直接用np，对向量中的值进行计算。
chi1_c=np.exp(c1_sigma_norm*chi1_1[0]+c1_mean_norm)
chi1_phi=np.exp(phi1_sigma_norm*chi1_1[1]+phi1_mean_norm)

5.利用for循环批量生成INP文档

通过上述过程，一方面提取出INP文件中的固定内容，另一方面生成了服从对数正态分布的样本点，将他们进行字符串拼接，就可以得到一个新的INP文件。利用for循环生成多个INP文件即可。代码大意如下：

numbers = input("请输入需要生成INP文件的个数：")
for i in range(numbers):**************************1.提取固定内容模块********************************************************2.生成对数正态分布模块**************************************************************3.生成材料参数字符串模块***********************************************************************4.拼接字符串并输出为inp文件模块***************************************

至此，就完成了INP文件的批量生成。为可靠度计算提供了计算基础。

六、批量计算INP文件

得到INP文件后，需要启动Abaqus进行批量计算，此过程需要用到批处理功能，本节将对此进行介绍。
首先在INP文件所在路径建立一个txt文件，将以下内容写入：

echo ------Job begins!
echo. #echo.表示回车
echo.echo ------Job 1 .............................................. Job 1
echo.call abaqus job=1 cpus=1  interactive #使用abaqus内核计算echo.
echo.
echo ------Job 2 .............................................. Job 2
echo.call abaqus job=2 cpus=1  interactive
echo.
echo.

此时我的INP文档里只有两个INP文件

再将txt文件后缀改为.bat，双击此文件即可运行批处理功能。
书写BAT文档可以用python字符串拼接完成，在这里不赘述，遇到具体问题可以留言。

七、批量提取ODB文件并进行可靠度计算

经过上述操作会得到计算结果文件，也就是ODB文件，ODB文件中包含计算以后得到的重要内容，包括应力应变，场变量等相关信息。我们在此部分讲解如何提取出安全系数这一结果。首先说明我提取结果文件的思路，这也许不是最好的方法，如果大家有更好的方法，还请在评论区说一下供大家学习。
思路
首先我利用一个py文件完成了一次ODB文件的提取，之后又创建了N个py文件，并用一个新的py文件逐个启动之前建立的N个py文件，让他们依次提取ODB中的安全系数。

1.利用py文件提取odb中的场变量

我这里只进行两个ODB文件的提取，将ODB文件剪切到另一个文件夹中，文件夹中有py文件和runall.py文件，如下图所示

1.py文件中的代码如下：

#!/user/bin/python
# -*- coding:UTF-8 -*-
from odbAccess import*
from odbSection import*
import os
myodb=openOdb(path="1.odb") #odb文件的名称
cpFile=open("1.txt","w") #建立新的txt文件
##提取场变量值，注意分析步名称！！！！！要和模型中的分析步名称一致，
##我们要提取的是强度折减分析步的场变量。
RS=myodb.steps['Step-2'].frames[-1].fieldOutputs['FV1'].values
#将场变量写入txt文档中
for i in range(1):cpFile.write('%s\n'%(RS[i]))
else:cpFile.close()

runall.py文件中的代码如下：

#!/user/bin/python
# -*- coding:UTF-8 -*-
import os
os.system("python ./1.py")
os.system("python ./2.py")

2.批量提取ODB结果文档

在cmd命令中，先进入odb文档所在路径，然后执行命令行即可：

提取后的文件夹中出现了1.txt和2.txt两个文件其中的内容是

其中1.txt文件中的内容如下：

({'baseElementType': 'CPE4P', 'conjugateData': None, 'conjugateDataDouble': 'unknown', 'data': 1.04779052734375, 'dataDouble': 'unknown', 'elementLabel': 1, 'face': None, 'instance': 'OdbInstance object', 'integrationPoint': 1, 'inv3': None, 'localCoordSystem': None, 'localCoordSystemDouble': 'unknown', 'magnitude': None, 'maxInPlanePrincipal': None, 'maxPrincipal': None, 'midPrincipal': None, 'minInPlanePrincipal': None, 'minPrincipal': None, 'mises': None, 'nodeLabel': None, 'outOfPlanePrincipal': None, 'position': INTEGRATION_POINT, 'precision': SINGLE_PRECISION, 'press': None, 'sectionPoint': None, 'tresca': None, 'type': SCALAR})

‘data’: 1.04779052734375即为安全系数。
接下来只要对txt文档进行合并，并且剔除多余内容，仅留下安全系数内容，就完成了ODB文档的提取过程。

八、可靠度结果后处理

可靠度结果的后处理，实际上就是对计算结果进行统计，当我们使用蒙特卡罗方法计算可靠度时，会生成一些列的安全系数，假设有1000组数据，得到了1000个安全系数，我们通过对这1000个安全系数进行统计分析，就会发现一些数据上的规律。本文用Matlab进行后处理。
常见的统计特征有前四阶矩，他们分别是均值、标准差、偏度和峰度。我会在下面展示这前四阶矩的相关代码(Matlab)：

Z=[1,2,3,4,5,6,7,8,9]jun_Zhi=mean(Z)biao_zhun_cha=std(Z)pian_du=skewness(Z)feng_du=kurtosis(Z)

所得结果

jun_Zhi=5
biao_zhun_cha=2.7386
pian_du=0
feng_du=1.7700

1.CDF、PDF曲线绘制

画累积分布曲线

Z1=randn(2000000,1)
figure(1);
h2=cdfplot(Z1);

画概率密度曲线

Z1=randn(2000000,1)
hist_number = 100
[N,Y]=hist(Z1,hist_number);
N1=N./sum(N);
plot(Y,N1);

对于概率密度曲线来说，数据越多，图像就越平滑。

2.失效概率计算

可靠度计算相关代码（Matlab）：

FS1=[1.1,0.9,1.2,0.8,1.3,0.7,1.5,0.6];
Function=[];
for i = 1:length(FS1)if FS1(i)-1 <=0Function1(i)=1;elseFunction1(i)=0;end
end
Pf4=vpa((sum(Function1)/length(FS1)),10)

结果：

Pf4 = 0.5 #失效概率为50%

后记

基于Abaqus 的边坡可靠度计算到这里就介绍的差不多了，今后会进行小的错误修改和更新，希望转载的同学表明出处支持原创。
此外，蒙特卡洛方法进行边坡可靠度计算虽然结果较为精确，但是其代价是耗费大量时间，所以有学者提出了代理模型，我所熟悉的代理模型是随机响应面方法，在今后的文章中会进行详细介绍。我还编写了自动生成随机多项式的程序辅助最小二乘拟合。

另外，基于空间变异性的随机场模型可以更加客观的讨论边坡可靠度的实际情况，通过K-L级数离散随机场，对每个单元进行参数赋值，实现模型的计算，如果有机会，我也会对此方法进行介绍。

最后感谢大家的耐心阅读。我也会进行下一阶段的学习，写出质量更多的教程。

基于Abaqus的边坡可靠度计算相关推荐

python自动化算法_基于Python语言和Abaqus平台的边坡可靠度计算自动化算法开发
2.2 求解过程求解过程包括7步,如图 2 图 2 自动化程序的求解过程 Fig. 2 The solving part of the automation program 1) Python形成 ...
matlab计算边坡可靠度,基于最优化方法的边坡稳定可靠度计算分析
0引言岩土性质在时间.空间上存在着可变性,传统的确定性分析方法只考虑了岩土材料各参数的均值,无法准确地描述各参数实际存在的不确定性和复杂性.而基于可靠度理论的极限状态分析方法采用了概率论和数理统计的方 ...
abaqus编写本构方程vumat_基于ABAQUS的木材本构关系数值模拟方法与流程
本发明属于木材本构研究技术领域,提供了木材在复杂应力状态下应力-应变关系的数值模拟方法. 背景技术: 木材是各向异性材料,其本构的复杂主要表现为在受压作用下发生塑性变形,而在拉.剪作用下发生脆性破坏, ...
基于abaqus的各向异性材料的抗拔力学性能分析
基于abaqus的各向异性材料的抗拔力学性能分析最近在做有关木材的抗拔性能的分析,遇到了一系列的问题,比如各向异性材料参数填写,材料方向分配,基于HILL屈服准则的材料参数等,经过一段时间的文献调研 ...
abaqus编写本构方程vumat_基于Abaqus子程序的高分子材料本构关系实现
摘要: 对于高分子材料的仿真,业界一般使用经典的弹塑性本构模型来描述其应力应变关系, 但其真实的应力应变关系与经典的弹塑性本构模型存在一定差异,从而导致仿真与实际测试之间的差异.Abaqus提供UMA ...
abaqus截面惯性矩_基于ABAQUS的路灯灯杆抗风能力校核
基于ABAQUS的路灯灯杆抗风能力校核摘要:出现台风等大风灾害时,因为路灯灯杆抗风能力不够导致路灯破坏会给人们的生命财产安全带来巨大的危害.因此根据路灯安装地区的风力情况对设计的路灯灯杆进行抗风能 ...
用随机场理论和卷积神经网络（CNN）分析边坡可靠度
近期,深度学习变得越来越火热,用深度学习跨领域研究一些其他领域的问题可以开拓研究思路.比如说用深度学习中的人工神经网络(ANNs)和卷积神经网络(CNN)去分析边坡可靠度的问题.边坡可靠度的研究是在假 ...
基于Abaqus的随机纤维增强复合材料拉伸试样建模插件2.0
复合材料研究是目前一个较为热门的方向,复合材料主要分为:①纤维增强复合材料②夹层复合材料③颗粒复合材料④混杂复合材料:对于纤维增强复合材料来说,又分为连续增强复合材料.短纤维增强复合材料.短纤维增强复 ...
基于ABAQUS二次开发的仿真分析平台
✨基于ABAQUS二次开发的仿真分析平台✨ 随着近年来计算机领域里程碑式的进步,计算机软件市场的迅速扩张,推出了许多功能强大的计算机仿真软件.ABAQUS有限元仿真分析软件则是其中的翘楚,作为应用广泛 ...

基于Abaqus的边坡可靠度计算

文章目录

基于Abaqus的边坡可靠度计算

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

1.在CAE中建立边坡模型

2.将CAD模型导入CAE

3.建模过程

4.后处理

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

1.关键字

2.数据

3.强度折减法关键字

三、抽样方法

1.随机数的生成

2.拉丁超立方抽样

四、参数相关性

1.相关系数矩阵

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

五、通过Python批量生成INP文件

1.Python字符串处理基础知识

2.关于强度折减法土体随机参数的生成

3.服从对数正态分布参数的生成

5.利用for循环批量生成INP文档

六、批量计算INP文件

七、批量提取ODB文件并进行可靠度计算

1.利用py文件提取odb中的场变量

2.批量提取ODB结果文档

八、可靠度结果后处理

1.CDF、PDF曲线绘制

2.失效概率计算

后记

基于Abaqus的边坡可靠度计算相关推荐

最新文章

热门文章