一种基于ABAQUS的工程结构响应面随机有限元分析处理方法与流程

本发明涉及工程数据处理应用技术领域，更为具体地，涉及一种基于ABAQUS的工程结构响应面随机有限元分析处理方法。

背景技术：

计算机辅助工程(Computer Aided Engineering，简称为CAE)分析是以有限元法为基础发展起来的一种十分有效的计算机数值仿真与优化设计技术。有限元法是一种高效能、常用的数值计算方法，其基于确定性的参数进行计算，其中确定性的参数是指结构的材料性能、载荷及边界条件等都是确定的。但在工程实际中，结构的材料性能、载荷及边界条件等都存在着很大的统计随机性，而这些结构参数的统计随机性，会对结构的临界性能和可靠性有较大的影响，给分析结果带来误差。

技术实现要素：

本发明的目的在于克服现有技术的不足，提供一种能够减小误差，提高工程结构可靠度和灵敏度的基于ABAQUS的工程结构响应面随机有限元分析处理方法。

本发明提供了一种基于ABAQUS的工程结构响应面随机有限元分析处理方法，具体包括如下步骤：

步骤S1：建立工程结构参数化有限元模型；

步骤S2：根据工程结构参数化有限元模型，通过数据采集单元采集工程结构中存在随机分布特性的参数，形成参数集合X＝(x₁，x₂，...，x_n)，并根据工程标准、试验数据和实际经验，确定所述参数集合中每个参数的随机分布特性；

步骤S3：采用序列响应面法拟合极限状态方程g(X)：

其中，a、b_i、c_i分别为极限状态方程的二次项系数，x_i为原始空间的随机分布参数；

步骤S4：取参数集合在原始空间中的均值点为验算点的初值并通过Rackwits-Fiessler算法将随机分布参数x_i当量正态化；其中，在验算点处，当量正态随机变量的累积分布函数与原随机变量的累积分布函数相等，以及当量正态随机变量的概率密度函数与原随机变量的概率密度函数相等；

式(1)中，为标准正态空间的随机分布参数；

式(2)中，为标准正态概率密度函数；

式(3)中，和分别为随机分布参数x_i对应的近似正态分布函数的均值和方差；

根据式(1)～(3)，求出近似正态分布函数的均值和方差

步骤S5：将极限状态响应面方程函数转化到标准正态空间变为：

其中，g'(U)为转化到标准正态空间的极限状态方程函数，α₀、r_i、λ_i分别为极限状态方程函数的二次项系数，u_i为标准正态空间的随机分布参数；

步骤S6：取标准正态空间中拟合点的参数δ＝0.1～0.5，形成标准正态空间中的2n+1个拟合点，2n+1个拟合点分别为：

步骤S7：将标准正态空间中的2n+1个拟合点换算成原始空间中的拟合点，原始空间的拟合点为：

步骤S8：将原始空间的拟合点代入参数化有限元模型获得极限状态响应值g₁、g₂、…、g_2n、g_2n+1；

步骤S9：将响应值代入标准正态空间响应面拟合方程获得联立方程式，求出联立方程式的系数a₀、r₁、λ₁、…、r_n、λ_n；其中，

步骤S10：将联立方程式的系数a₀、r₁、λ₁、…、r_n、λ_n分别带入标准正态空间响应面拟合方程，获得响应面方程式：

步骤S11：对响应面方程式采用JC算法求得工程结构的可靠度指标β、结构可靠度对随机分布参数的灵敏度系数α_i和新验算点的初值其中，

JC算法的计算过程是先将极限状态方程展开为泰勒级数，在分别计算，由于JC算法是现有技术，故计算过程的细节在此不过多说明；

步骤S12：根据可靠度指标β和新验算点的初值重复步骤S7～S11直至可靠度指标β收敛，其中可靠度指标β收敛是指上次求出的β至上次求出的β≤允许误差ε，而允许误差ε根据工程结构的要求确定。

本发明的基于ABAQUS的工程结构响应面随机有限元分析处理方法，将可靠性计算方法与有限元法相结合建立的，能够求解工程复杂结构可靠度的问题，该方法不受随机参数变异性大小的限制，计算较为准确，具有高效、方便的特点。

附图说明

图1为基于ABAQUS的工程结构响应面随机有限元分析处理方法的流程示意图；

图2为悬臂梁的力学模型示意图；

图3为轴承座有限元网格模型；

具体实施方式

下面详细说明本发明的具体实施，有必要在此指出的是，以下实施只是用于本发明的进一步说明，不能理解为对本发明保护范围的限制，这些方面指示的仅仅是可使用本发明的原理的各种方式中的一些方式，本发明旨在包括所有这些方面以及它们的等同物，并且该领域技术熟练人员根据上述本发明内容对本发明做出的一些非本质的改进和调整，仍然属于本发明的保护范围。

本发明提供了一种基于ABAQUS的工程结构响应面随机有限元分析处理方法，首先，ABAQUS(有限元分析)是工程模拟的有限元软件，其解决问题的范围从相对简单的线性分析到复杂的非线性问题。本发明通过结合ABAQUS的二次开发技术，基于Python语言编程，实现工程结构的参数化建模，通过调用有限元分析结果建立工程结构的极限状态响应面方程，结构可靠度计算方法(优选为JC算法)来计算结构的可靠度和分析参数的敏感度，从而建立一种能有效地解决实际工程复杂结构可靠性分析问题的响应面随机有限元分析处理方法。

图1示出了根据本发明的基于ABAQUS的工程结构响应面随机有限元分析处理方法的流程示意图，具体包括如下步骤：

步骤S1：建立工程结构参数化有限元模型；

由于ABAQUS软件程序的主体框架均由Python语言构成，利用Python才可以对ABAQUS软件进行二次开发，所以需要基于Python语言建立参数化有限元模式；

步骤S2：根据工程结构参数化有限元模型，通过数据采集单元采集工程结构中存在随机分布特性的参数，形成参数集合X＝(x₁，x₂，...，x_n)，并根据工程标准、试验数据和实际经验，确定所述参数集合中每个参数的随机分布特性；

采集工程结构中存在随机分布特性的参数，为采集工程结构的尺寸、材料参数、载荷及边界条件等存在着随机分布特性的参数，参数的随机分布特性是根据工程实际经验或试验测试统计获得的。确定随机分布特性是为了掌握参数的随机分布情况，在此基础上才能进行随机有限元分析，计算可靠性。

步骤S3：采用序列响应面法拟合极限状态方程g(X)：

其中，a、b_i、c_i分别为极限状态方程的二次项系数，x_i为原始空间的随机分布参数；

步骤S4：取参数集合在原始空间中的均值点为验算点的初值并通过Rackwits-Fiessler算法将随机分布参数x_i当量正态化；其中，在验算点处，当量正态随机变量的累积分布函数与原随机变量的累积分布函数相等，以及当量正态随机变量的概率密度函数与原随机变量的概率密度函数相等；

式(1)中，为标准正态空间的随机分布参数；

式(2)中，为标准正态概率密度函数；

式(3)中，和分别为随机分布参数x_i对应的近似正态分布函数的均值和方差；

根据式(1)～(3)，求出近似正态分布函数的均值和方差

步骤S5：将极限状态响应面方程函数转化到标准正态空间变为：

其中，g'(U)为转化到标准正态空间的极限状态方程函数，α₀、r_i、λ_i分别为极限状态方程函数的二次项系数，u_i为标准正态空间的随机分布参数；

步骤S6：取标准正态空间中拟合点的参数δ＝0.1～0.5，形成标准正态空间中的2n+1个拟合点：

步骤S7：将标准正态空间中的2n+1个拟合点换算成原始空间中的拟合点，原始空间的拟合点为：

步骤S8：将原始空间的拟合点代入参数化有限元模型获得极限状态响应值g₁、g₂、…、g_2n、g_2n+1；

步骤S9：将响应值代入标准正态空间响应面拟合方程获得联立方程式，求出联立方程式的系数a₀、r₁、λ₁、…、r_n、λ_n；其中，

步骤S10：将联立方程式的系数a₀、r₁、λ₁、…、r_n、λ_n分别带入标准正态空间响应面拟合方程，获得响应面方程式：

步骤S11：对响应面方程式采用JC算法求得工程结构的可靠度指标β、结构可靠度对随机分布参数的灵敏度系数α_i和新验算点的初值其中，

步骤S12：根据可靠度指标β和新验算点的初值重复步骤S7～S11直至可靠度指标β收敛，其中可靠度指标β收敛是指上次求出的β至上次求出的β≤允许误差ε，而允许误差ε根据工程结构的要求确定。

上述内容详细说明了本发明提供的基于ABAQUS的工程结构响应面随机有限元分析处理方法的流程，下面将以两个实例对上述方法可取得的有益效果进行说明。

实例一

以某悬臂梁为例，如图2所示，悬臂梁长为L，梁的截面形状为矩形，尺寸为a×b。悬臂梁上表面受均布载荷，载荷大小为F，材料屈服应力为σ_s。分析中考虑结构尺寸a、b、L及载荷F、屈服强度σ_s五个变量的随机性，具体分布如表1所示：

表1随机变量的分布

根据材料力学知识，该悬臂梁的极限状态方程为：

采用序列响应面拟合标准正态空间的极限状态方程为：

g＝a₀+r₁A+λ₁A²+r₂B+λ₂B²+r₃l+λ₃l²+r₄f+λ₄f²+r₅s+r₅s²；

其中，A、B、l、f、s均为标准正态随机向量。

将最后求得的β值与Monte-Carlo抽样法或其它解析方法计算的β值进行比较，以Monte-Carlo抽样法的β值为精确解，计算允许误差，比较结果的分析如表2所示：

表2可靠度指标β的结果分析

通过表2可以看出，由表中的数据可知，通过本发明提供的方法计算得到的β与精确解相比，允许误差只有1.997％，虽然比二次二阶矩法的精度低，但比一次二阶矩法的均值法要精确，而1.997％的相对误差在大多数工程实际中是可以接受的。同时，该方法只需根据变异性设定随机变量参数即可，不受随机变量变异性大小的限制，适用性广。

实例二

某轴承座通过4个安装孔进行固定，轴承孔的下半部分承受由轴传来的径向压力载荷P，轴承孔圆周上承受推力载荷P/5。轴承座材料为钢，弹性模量为E，泊松比μ＝0.3，屈服强度为σ_s。分析中考虑径向压力载荷P₁、推力载荷P₂、弹性模量E、屈服强度σ_s四个变量的随机性。

根据轴承座的几何模型建立轴承座的有限元模型，建立的轴承座的有限元模型如图3所示。

将最后求得的β值与Monte-Carlo抽样法的β值进行比较，以MCSFEM法的β值为精确解，计算允许误差，比较结果的分析如表3所示：

表3可靠度指标β的结果分析

由表中的数据可知，本发明提供的方法与MCSFEM法(10000次抽样实验)计算得到的精确解相比误差只有1.93％，而计算时间只有MCSFEM法的1.5％，证明本发明通过的方法计算较为准确，计算效率很高，可用于工程实际。

尽管为了说明的目的，已描述了本发明的示例性实施方式，但是本领域的技术人员将理解，不脱离所附权利要求中公开的发明的范围和精神的情况下，可以在形式和细节上进行各种修改、添加和替换等的改变，而所有这些改变都应属于本发明所附权利要求的保护范围，并且本发明要求保护的产品各个部门和方法中的各个步骤，可以以任意组合的形式组合在一起。因此，对本发明中所公开的实施方式的描述并非为了限制本发明的范围，而是用于描述本发明。相应地，本发明的范围不受以上实施方式的限制，而是由权利要求或其等同物进行限定。