CN109002635B - 一种考虑连接非线性的螺栓连接等效建模方法 - Google Patents

Abstract

本发明提供了一种考虑连接非线性的螺栓连接等效建模方法，根据待分析对象的几何特征建立三维实体网格，确定螺栓所在位置需要赋予螺栓连接非线性特征的单元，定义基于Iwan模型的弹塑性随动强化材料，在有限元软件中，使用自定义本构的相关模块定义该材料，得到可使用的材料属性，取出螺栓所在位置希望赋予非线性连接特性的网格，赋予所定义的材料属性，其他部分仍赋予原定的材料属性，根据实际分析需求进行有限元计算。本发明基于Iwan模型力学关系，通过商用有限元软件提供的二次开发接口定义基于Iwan模型的随动强化塑性材料本构，再将该材料赋予连接部位单元以实现连接非线性滞回特性简化建模。本发明方法能够适用于螺栓承受较大载荷、变形等出现界面微滑移状态，计算操作简便。

Description

一种考虑连接非线性的螺栓连接等效建模方法

技术领域

本发明涉及一种非线性连接的螺栓，具体涉及一种考虑连接非线性的螺栓连接等效建模方法。

背景技术

对大型、复杂结构进行静、动力学分析一般需要通过有限元方法，即利用商用有限元软件完成，如ABAQUS、ANSYS、Nastran等。进行有限元分析时，如果没有具体分析螺栓力学特性的必要，则往往将螺栓连接部分简化为线性单元(如单向弹簧、多向弹簧、薄层单元等)，以减少有限元分析的计算压力。

螺栓连接部位具有复杂的摩擦接触关系，因此具有很强的非线性特征。线性的简化方法参数少、关系简单且操作简便，可以通过参数修正等方法达到一定精度。但这些特点也使得该方法难以在需要考虑螺栓连接受较大荷载、力-位移关系产生显著的非线性的情况下应用。

现有文献中已提出了多种考虑螺栓连接非线性的简化建模方法，其中比较著名的为基于Iwan模型的有限元模拟方法。在应用这些建模方法时，这些文献往往选择通过数值计算、自行编写有限元计算程序等方法实现，目前缺乏一种简便易行、且适用于商业有限元软件分析的建模方法，也就无法将非线性建模方法从概念推广到复杂实际结构分析中。

发明内容

发明目的：本发明的目的在于针对现有技术的不足，提供一种考虑连接非线性的螺栓连接等效建模方法，可通过商用有限元软件(如ABAQUS、ANSYS、Nastran)实现对螺栓连接结构的非线性有限元分析，计算效率和计算精度高。

技术方案：本发明提供了一种考虑连接非线性的螺栓连接等效建模方法，包括以下步骤：

(1)根据待分析对象的几何特征建立三维实体网格；

(2)确定螺栓所在位置需要赋予螺栓连接非线性特征的单元；

(3)定义基于Iwan模型的弹塑性随动强化材料；

(4)在所用的有限元软件中，使用自定义本构的相关模块，定义步骤(3)的材料，得到可使用的材料属性；

(5)取出螺栓所在位置希望赋予非线性连接特性的网格，薄层或一般实体单元均可，赋予步骤(4)所定义的材料属性，其他部分仍赋予原定的材料属性；

(6)根据实际分析需求进行有限元计算。

进一步，步骤(3)包括以下步骤：

(31)将Iwan模型的力-位移关系式前段近似线性的部分改写为线性形式；

(32)将Iwan模型的力-位移关系式转化为应力-应变关系式；

(33)利用弹塑性随动强化材料本构关系，求解各参数并带入本构关系方程，得到的最终表达式即为基于Iwan模型的弹塑性随动强化材料。

进一步，步骤(31)力-位移关系式如下：

其中，x为位移，F(x)为受力，K为受力方向刚度，α为完全滑移后的残余刚度，f_q为完全进入滑移时对应拉力类比为屈服拉力，x_s＝2f_q/K为力-位移关系转捩点；

式中0≤x≤x_s部分的前段近似为直线，因此可确定一近似直线段，再将其取为直线：记x＝x_l、y＝F(x)、y＝0三条曲线包围的面积为S₀，x＝x_l、y＝F(x_l)·x/x_l、y＝0三条曲线包围的面积为S₁，求解S₁＝95％S₀时x_l的值作为对应x_l的解，将解得的x_l作为材料的弹性极限，此时

据此将式(1)改写为

式中0≤x≤xl段为弹性段、其余为塑性段。

进一步，步骤(32)虚设或根据步骤(2)中单元域的尺寸设定一个满足材料力学假设的梁，令A和l分别为该虚设梁的长度及截面积，有

式中，σ,ε,ε^p分别为应力、全应变、塑性应变；

将式(4)带入式(3)，并将力变量F替换为应力σ、位移x替换为塑性应变ε^p，即通过力-位移关系式F＝F(x)得到应力-应变关系式σ＝σ(ε^p)及其一、二阶导数，这三个关系式即为Iwan模型的特征关系式。

进一步，步骤(33)在有限元软件计算中使用增量式的本构关系，得到增量率型本构关系方程：

式中使用张量标记，双下标为爱因斯坦求和约定记号，为应力率，为应变率，应力率和应变率为因、自变量，λ^*、μ^*为过渡参数；式中已知参数有：材料泊松比μ、体积弹性模量克罗内克符号δ_ij、交错单位张量η_ij及η_kl，其中E为材料弹性模量；未求解的参数有：屈服应力σ_Y、塑性模量E_p、强化方向η_ij、有限元单元等效应力有限元单元等效应变及塑性模量E_p关于等效应变的导数H；

其中，屈服应力σ_Y为式(3)的力F(x_l)通过式(4)计算得到的应力；为通过牛顿迭代法确定的相邻两步等效应变与之差，为迭代法上一步算得的值；

参数E_p、η_ij、通过以下方法获得：

令该计算增量步的有限元单元等效应力为塑性模量E_p为关系式的导数，这里的等价于步骤(32)中的σ＝σ(ε^p)，即等同于σ＝σ(ε^p)：

当单元等效应力超过屈服点σ_Y时，求解屈服条件：

由于E_p、是等效塑性应变的函数，通过牛顿迭代法求解式(7)，解得E_p、和代入下式解得该单元应力σ_ij和强化方向η_ij：

式中，为偏应力张量，α_ij为本迭代步算得的背应力，Δα_ij为本迭代步的背应力增量，为上一迭代步算得的α_ij，为单元应变增量，σ_ij为单元应力，为通过爱因斯坦求和约定表示的单元塑性应力；

H为E_p关于的导数，即：

有益效果：本发明基于Iwan模型力学关系，通过商用有限元软件提供的二次开发接口定义基于Iwan模型的随动强化塑性材料本构，再将该材料赋予连接部位单元以实现连接非线性滞回特性简化建模。同传统线性简化等效建模方法相比，该方法能够适用于螺栓承受较大载荷、变形等出现界面微滑移状态，该方法兼容商用有限元软件、计算操作简便。

附图说明

图1为赋予螺栓连接非线性区域举例(连接区域)；

图2为有限元螺栓连接非线性模拟方法流程图；

图3为ABAQUS UMAT实现流程图。

具体实施方式

下面对本发明技术方案进行详细说明，但是本发明的保护范围不局限于所述实施例。

一种考虑连接非线性的螺栓连接等效建模方法，如图2所示，结合ABAQUS及其用户定义子程序(User Subroutine in ABAQUS)进行进一步说明：

1、根据需分析的对象的几何特征建立3D实体网格。

2、确定螺栓所在位置需要赋予螺栓连接非线性特征的单元。

3、定义基于Iwan模型的弹塑性随动强化材料本构关系方程。

31)改写Iwan模型力-位移关系式：

Iwan模型描述了一种模型材料受拉力作用后施力点位移x与拉力大小F的关系，该关系也可以推广到一般的受力-受力点位移关系。力-位移关系式如下：

其中，K为受力方向刚度，α为完全滑移后的残余刚度，f_q为完全进入滑移时对应拉力(类比为屈服拉力)，x_s＝2f_q/K为力-位移关系转捩点、可类比为材料屈服点。

该式不存在一般材料定义中具备的线弹性段，但函数0≤x≤x_s部分的前段近似为直线，因此可确定一近似直线段，再将其取为直线。记x＝x_l、y＝F(x)、y＝0三条曲线包围的面积为S₀，x＝x_l、y＝F(x_l)·x/x_l、y＝0三条曲线包围的面积为S₁，求解S₁＝95％S₀时x_l的值作为对应x_l的解。将解得的x_l作为材料的弹性极限，此时

据此将式(1)改写为

式中0≤x≤x_l段为弹性段、其余为塑性段。

32)将Iwan模型力-位移关系式，转化为应力-应变关系式。虚设(或根据步骤2中选定的单元域的尺寸设定)一满足材料力学假设、两端受力F的梁，令A和l分别为该虚设梁的长度及截面积，有

式中，σ,ε,ε^p分别为应力、全应变、塑性应变，F为式(3)所述Iwan模型所受力。

将式(4)代入式(3)，将式(3)的因变量F、自变量x分别置换为σ、ε^p，即可得到应力-应变关系式

其导数为

二阶导数为

这三个关系式即为Iwan模型的特征关系式。此处的σ＝σ(ε^p)即等价于步骤33)的

33)利用通用的弹塑性随动强化材料本构关系，构造基于Iwan模型的弹塑性随动强化材料。商用有限元软件一般不提供直接输入复杂随动强化材料本构关系的选项，但提供编程接口定义复杂本构关系。有限元软件计算中使用增量式的本构关系，常使用的增量率型本构关系方程表达为：

式中使用张量标记，式中双下标为爱因斯坦求和约定记号、不作为变量区分依据。为应力率，为使用爱因斯坦求和约定记号的不同表述形式的应变率，应力率、应变率两者为因、自变量。

式中已知参数有：材料泊松比μ、体积弹性模量(E为材料弹性模量)、克罗内克符号δ_ij、交错单位张量η_ij及η_kl。

未求解的参数有：屈服应力σ_Y、塑性模量E_p、强化方向η_ij、有限元单元等效应力有限元单元等效应变及塑性模量E_p关于等效应变的导数H。

其中，H为塑性模量E_p关于的导数，即

σ_Y为材料屈服点，将式(3)的力F(x_l)通过式(4)变换得到的应力σ即为σ_Y。

为通过牛顿迭代法确定的相邻两步等效应变与之差，为迭代法上一步算得的值。

参数E_p、η_ij、可通过以下方法获得：令该计算增量步的有限元单元等效应力(常用Mises应力表达)为塑性模量E_p为关系式(等价于步骤(32)描述的σ＝σ(ε^p)，即等同于σ＝σ(ε^p))的导数：

当单元等效应力超过屈服点σ_Y时，求解屈服条件

由于E_p、是等效塑性应变的函数，通过牛顿迭代法求解式(11)，解得E_p、和代入下式

解得该单元应力σ_ij和强化方向η_ij；式中为偏应力张量，α_ij为本迭代步算得的背应力、Δα_ij为本迭代步的背应力增量、为上一迭代步算得的α_ij，为单元应变增量，σ_ij为单元应力，为通过爱因斯坦求和约定表示的单元塑性应力。

34)步骤32)式(5)描述的σ(ε^p)即为式(10)的那么E_p、H的表达式可分别由式(10)、式(9)得到。将η_ij、E_p及H的表达式代入步骤33)式(8)，即得到最终的基于Iwan模型的弹塑性随动强化材料。

4、在所用的有限元软件中，使用自定义材料本构关系的相关模块，定义步骤3所述的材料、得到可使用的材料属性。在ABAQUS中，操作如下：

41)预先配置好可以运行用户定义子程序的环境。

42)新建一个Fortran文件，并使用Fixed Form格式进行编程。

43)基于步骤31)、32)所述理论构造函数(为称呼方便，此处命名为IWANHARD)。该函数接口参数包括：屈服应力、E_p、H、等效塑性应变Iwan模型参数(α、f_q、A、l、材料弹性模量E)。前三项为输出项，其余为输入项。

44)基于步骤33)所述理论构造用户定义材料子程序(UMAT)函数。该函数接口、结构等应遵循ABAQUS的相关规范。应在UMAT函数中声明需要输入的Iwan模型参数。在需要Iwan模型特征参数(即步骤43)所述屈服应力、E_p、H、等效塑性应变)时，通过调用IWANHARD函数、输入参数后获取，算法流程见图3。

45)将上述两个函数前后保存在步骤42)建立的文件中，确保格式无误。

5、取出螺栓所在位置希望赋予非线性连接特性的网格，赋予步骤4所定义的材料属性。赋予材料属性的网格应处于连接区域、类型不拘，可直接施加到连接处的一般实体单元上(如图1灰色区域)，其他部分仍赋予原定的材料属性。在ABAQUS的前处理模块中，需按步骤4定义的参数顺序输入参数，并给出每一增量步间需传递的参数个数。

6、根据实际分析需求进行有限元计算。在ABAQUS任务提交模块，选择步骤4定义好的子程序Fortran文件。

Claims (3)

1.一种考虑连接非线性的螺栓连接等效建模方法，其特征在于：包括以下步骤：

(1)根据待分析对象的几何特征建立三维实体网格；

(2)确定螺栓所在位置需要赋予螺栓连接非线性特征的单元；

(3)定义基于Iwan模型的弹塑性随动强化材料；

(4)在所用的有限元软件中，使用自定义本构的相关模块，定义步骤(3)的材料，得到可使用的材料属性；

(5)取出螺栓所在位置希望赋予非线性连接特性的网格，薄层或一般实体单元均可，赋予步骤(4)所定义的材料属性，其他部分仍赋予原定的材料属性；

(6)根据实际分析需求进行有限元计算；

具体的，步骤(3)包括以下步骤：

(31)将Iwan模型的力-位移关系式前段近似线性的部分改写为线性形式；

(32)将Iwan模型的力-位移关系式转化为应力-应变关系式；

(33)利用弹塑性随动强化材料本构关系，求解各参数并带入本构关系方程，得到的最终表达式即为基于Iwan模型的弹塑性随动强化材料；

步骤(31)力-位移关系式如下：

其中，x为位移，F(x)为受力，K为受力方向刚度，α为完全滑移后的残余刚度，f_q为完全进入滑移时对应拉力类比为屈服拉力，x_s＝2f_q/K为力-位移关系转捩点；

式中0≤x≤x_s部分的前段近似为直线，因此可确定一近似直线段，再将其取为直线：记x＝x_l、y＝F(x)、y＝0三条曲线包围的面积为S₀，x＝x_l、y＝F(x_l)·x/x_l、y＝0三条曲线包围的面积为S₁，求解S₁＝95％S₀时x_l的值作为对应x_l的解，将解得的x_l作为材料的弹性极限，此时

据此将式(1)改写为

式中段为弹性段、其余为塑性段。

2.根据权利要求1所述的考虑连接非线性的螺栓连接等效建模方法，其特征在于：步骤(32)虚设或根据步骤(2)中单元域的尺寸设定一个满足材料力学假设的梁，令A和l分别为假设的梁的长度及截面积，有

式中，σ,ε,ε^p分别为应力、全应变、塑性应变；

将式(4)带入式(3)，并将力变量F替换为应力σ、位移x替换为塑性应变ε^p，即通过力-位移关系式F＝F(x)得到应力-应变关系式σ＝σ(ε^p)及其一、二阶导数，这三个关系式即为Iwan模型的特征关系式。

3.根据权利要求2所述的考虑连接非线性的螺栓连接等效建模方法，其特征在于：步骤(33)在有限元软件计算中使用增量式的本构关系，得到增量率型本构关系方程：

式中使用张量标记，双下标为爱因斯坦求和约定记号，为应力率，为应变率，应力率和应变率为因、自变量，λ^*、μ^*为过渡参数；式中已知参数有：材料泊松比μ、体积弹性模量克罗内克符号δ_ij、交错单位张量η_ij及η_kl，其中E为材料弹性模量；未求解的参数有：屈服应力σ_Y、塑性模量E_p、有限元单元等效应力有限元单元等效塑性应变及塑性模量E_p关于等效塑性应变的导数H；

其中，屈服应力σ_Y为式(3)的力F(x_l)通过式(4)计算得到的应力；为通过牛顿迭代法确定的相邻两步等效塑性应变与之差，为迭代法上一步算得的值；

参数E_p、通过以下方法获得：

令计算增量步的有限元单元等效应力为塑性模量E_p为关系式的导数，这里的等价于步骤(32)中的σ＝σ(ε^p)，即等同于σ＝σ(ε^p)：

当有限元单元等效应力超过屈服点σ_Y时，求解屈服条件：

由于E_p、是等效塑性应变的函数，通过牛顿迭代法求解式(7)，解得E_p、和代入下式解得单元应力σ_ij和强化方向n_ij：

式中，为偏应力张量，α_ij为本迭代步算得的背应力，Δα_ij为本迭代步的背应力增量，为上一迭代步算得的α_ij，为单元应变增量，σ_ij为单元应力，为通过爱因斯坦求和约定表示的单元塑性应力；

H为E_p关于的导数，即：