复合材料壁板螺栓连接的非线性分析方法与流程

本发明主要属于非线性有限元分析，具体涉及复合材料壁板螺栓连接的非线性分析方法。

背景技术：

目前，航空制造业广泛采用复合材料用于民用飞机的设计与制造。复合材料是非金属材料，并具有各向异性的特点。因此，广泛应用于合金材料结构件连接的传统铆钉铆接技术不适用于复材结构连接。复合材料结构件选择各种螺栓，例如高锁螺栓或沉头螺栓，进行连接。这给民机结构件的连接设计带来了新的挑战。机翼和机身结构中常见的复合材料壁板连接设计，需要准确计算通过螺栓所传递的剪切载荷，进而准确计算复合材料壁板连接区应力分布，确定应力集中区域并进行强度校核。确保在极限载荷情况下，复合材料壁板连接结构安全不破坏。需要指出的是，螺栓连接的特殊性在于其通过连接件的相互挤压和变形，即复合材料壁板与螺栓之间的接触传递载荷，是一种强非线性问题。

现有技术方案中，广泛应用的技术方案一将螺栓连接简化为固定连接，忽略接触。采用nastan中的一维钉单元或者abaqus中的一维cbush单元同时模拟螺栓弯曲刚度以及螺栓、复合材料壁板之间接触刚度的综合效应。总计6个刚度分量均由用户自己定义。通常，使用实验手段获得的经验公式确定上述刚度，最常用的公式为huth公式。但是该公式只能确定剪切刚度k2、k3，无法确定扭转刚度。显然，该技术方案对原始问题进行了过度简化，适用范围有限。该技术方案将螺栓连接简化为固定连接，可以实现复合材料壁板连接的快速应力分析。由于忽略了复合材料壁板开孔，预测出的应力分布是不准确的，也无法确定孔边应力状态。即使快速完成了相关的应力校核工作，该结果也需要进行相当保守的修正后才能用于民机的结构设计。

现有技术方案二采用abaqus中的三维实体单元或二维壳单元建立复合材料壁板的精细模型，采用三维实体单元建立螺栓、螺帽的精细数值模型。严格定义螺栓与复合材料壁板之间的接触面、复合材料壁板之间的接触面、螺帽与复合材料壁板之间的接触面，将其提取后定义成接触对。每一个接触对由接触单元描述。abaqus中的接触单元将接触条件作为一类特殊的不连续约束来处理。当两个表面接触时应用接触条件，并传递接触载荷，对载荷大小未做任何限制。当接触载荷变为零或负值时，接触面分离，约束被撤销。这种处理接触的方式被称为“硬”接触，忽略了两个物体接触时相互之间的材料渗透(inner-penetration)。这种渗透是一种真实物理现象，将会缓冲接触时的冲击载荷。虽然这种方案能够完成螺栓连接的应力分析，由于接触单元的引入，使问题变得非常复杂，需要耗费大量的计算机资源。在目前的民机设计中，还无法实现大规模的接触计算，广泛应用的仍然是现有技术方案一。

技术实现要素：

针对上述问题，本发明提供了复合材料壁板螺栓连接的非线性分析方法。

本发明的分析方法建立了hertz接触单元，用接触单元模拟复合材料壁板螺栓连接中的接触现象，解决技术方案一中复合材料壁板螺栓连接的过度简化问题，解决技术方案二中复合材料壁板螺栓连接的过度复杂、计算效率低的问题。

本发明是通过以下技术方案实现的：

一种复合材料壁板螺栓连接的非线性分析方法，所述方法包括以下步骤：

(1)分别对复合材料壁板和螺栓建立几何精确壳单元和梁单元模型；

(2)建立hertz接触单元，所述hertz接触单元模拟复合材料壁板壳单元和螺栓梁单元之间的接触；根据复合材料壁板壳单元和螺栓梁单元接触时的材料相互渗透量计算接触载荷，以接触载荷分别作为复合材料壁板壳单元和螺栓梁单元的受力载荷分析复合材料壁板壳单元的变形,实现对复合材料壁板螺栓连接结构的非线性分析,可进一步得到该螺栓连接情况下复合材料壁板的各处应力分布及螺栓的传承比例。进一步地，所述步骤(2)具体为：

21)所述hertz接触单元将复合材料壁板和螺栓的接触简化为壳单元与梁单元之间局部的二维平面接触；

在所述二维平面，螺栓的截面为梁单元包括梁节点的刚性截面圆，所述复合材料壁板的可能接触面为多个壳单元边节点构成的弹性圆；根据hertz接触理论，将复合材料壁板和螺栓的接触认为是弹性接触，并随时间变化；接触初始时刻，假设螺栓与复合材料壁板孔边接触但无渗透，所以无接触力；此时，描述复合材料壁板孔边形状的所有壳单元边节点均为可能接触点；

22)根据复合材料壁板螺栓连接件变形过程中各可能接触点相对所述螺栓刚性圆截面中心的位移来判断接触是否发生；如果接触发生则计算复合材料壁板和螺栓之间的材料相互渗透距离；

23)基于hertz公式，利用材料相互渗透距离计算得到接触面之间的接触力、接触法向量、接触刚度和分别作用到壳单元和梁单元相关节点上的接触载荷；

24)以接触载荷分别作为壳单元和梁单元的受力载荷，根据有限元原理将壳单元或梁单元载荷、单元刚度以及整个复合材料壁板螺栓连接件所受外载荷分别集成，得到包含接触载荷和接触刚度在内的总体刚度矩阵和总体载荷矢量，形成复合材料壁板螺栓连接的变形控制方程；

25)通过求解所述变形控制方程，计算获得复合材料壁板螺栓连接的变形，进一步可获得由螺栓传递的剪切载荷以及当前计算步骤所有可能接触点相对螺栓刚性截面中心的位移矢量；

26)完成当前步骤后，计算得到的复合材料壁板螺栓连接的变形构型将作为下一时间步骤计算的初始值，在时间域内重复步骤23)-25)直至时间域内的计算全部完成为止。

进一步地，所述步骤(1)中

所述几何精确壳单元和梁单元，是几何非线性结构有限元单元；在推导所述几何精确壳单元和梁单元变形控制方程时，仅引入线性应力-应变关系一个假设，严格定义了应变-位移关系，能够描述薄壁结构的大变形，具有强非线性特征；

所述几何精确壳单元和梁单元采用增量表达式理论处理了任意大变形和有限转动问题；在理论推导中，定义了三种不同构型描述其在时间域内的非线性变形，分别为参考构型、初始构型和变形构型；其中，参考构型、初始构型为已知构型、变形构型为未知构型；数值计算时仅仅预测从初始构型到变形构型的变形增量和旋转增量；依据建模步骤完成当前时间步的计算后，当前的变形构型将作为下一时间步骤计算的初始值，在时间域内重复建模直至时间域内的计算全部完成为止。

进一步地，

所述几何精确壳单元和梁单元的建模过程遵循经典有限元理论的分析流程，包括结构有限元单元网格划分、赋予结构有限元单元材料属性、单元载荷向量和单元刚度矩阵生成、单元载荷向量和刚度矩阵总体集成、边界条件处理、外部载荷加载、基于newton-raphson迭代算法的变形计算、计算结果的后处理获得结构变形、应变和应力分布；所述结构有限元单元网格划分时，需要标明那些几何精确壳单元和梁单元将组成“软”接触模型；在所述单元载荷向量和单元刚度矩阵生成时，对将组成“软”接触模型的几何精确壳单元和梁单元做处理：根据“软”接触模型判断接触是否发生，如发生接触，额外计算接触载荷和接触刚度，并将接触载荷和接触刚度添加到单元载荷向量和单元刚度矩阵中。

进一步地，所述构建复合材料壁板壳单元具体为几何精确壳单元，能够处理任意大变形和任意转动；

所述几何精确壳单元具有静力学、模态分解和动力学分析功能；

所述几何精确壳单元依据弹性力学理论中green-lagrange应变张量的定义，通过高阶项的近似处理，获得适合于任意大位移、有限转动和大变形情况的应变-位移关系；同时，所述几何精确壳单元假设变形过程中，与壳中面垂直的材料线变形后仍然保持直线，不发生拉伸变形；

所述几何精确壳单元由壳的运动控制方程结合有限元理论得出；

其中壳的运动控制方程为

式中f和m代表壳中面单位面积上作用的外部力和力矩，h和g代表动量矩，ni，i,n3和mi，i，i＝1，2，为弹性力和力矩，b为转换矩阵，为单位矢量的旋转对称矩阵，rr0为旋转张量，符号(^·)代表相对于时间的一阶导数。

进一步地，所述构建螺栓的梁单元具体为考虑翘曲变形的几何精确梁单元，能够处理大变形和任意转动；该单元能够根据预测出的一维应变结果，反向运算可获得螺栓截面上的二维应变和应力分布；半离散形式的运动控制方程为：

式中f和t代表梁单元单位展长上作用的外部力和力矩，h和g代表动量矩，n1和m1为弹性力和力矩，ω为试函数，u0和u为梁截面上任一点的初始位置和当前位移，符号()′代表相对于位置坐标的一阶空间导数。进一步地，所述精确模拟复合材料壁板和螺栓之间接触现象的单元具体为hertz接触单元，是一种“软”接触模型，能够精确计算接触发生时的相互渗透量，准确预测接触刚度和接触载荷，进而预测复材壁板孔边应力，计算精度好，计算效率高。

根据hertz接触理论，假设在初始构型，螺栓与复合材料孔边接触但没有接触力,在变形构型中螺栓与复合材料孔边发生接触，进而产生相互渗透；接触产生的渗透可由下式计算q＝|ur|-|u0|，a＝-q；其中，u0为从螺栓随体坐标系原点指向该接触点的矢量；ur为局部坐标系中定义的相对位移矢量

公式中us为由于初始接触位置移动而产生的位移矢量；

根据非线性hertz接触理论，当q＜0时接触发生，法向接触力以及法向接触力虚功可写为

式中，k为接触刚度，可根据发生接触的两个物体的真实材料属性严格确定；为矢量us在螺栓随体坐标系中的映射；

根据虚功原理，可以获得如下接触载荷矢量

当q≥0时，接触载荷矢量fc＝0，此时认为螺栓与复合材料壁板孔边无接触。

进一步地，

利用步骤25)计算结束时得到的复合材料壁板螺栓连接件变形后的最终构型，获得所有离散单元节点上的位移和变形，以及接触单元上的材料渗透和接触载荷；后处理过程中，可以进一步获得复合材料壁板孔边的应力分布以及螺栓的传载比例。

本发明的有益技术效果：

(1)本申请的技术方案建立接触单元，精确模拟复合材料壁板螺栓连接中的接触现象，解决技术方案一中复合材料壁板螺栓连接的过度简化问题。

(2)本申请的技术方案基于penalty理论中的非线性hertz模型，准确确定连接刚度，解决技术方案一中无法准确确定连接刚度的问题。

(3)本申请的技术方案中开发的接触单元计入两个物体接触时相互之间的材料渗透，并准确计算接触载荷，同时考虑了接触引起的能量耗散，解决技术方案二中未考虑两个物体接触时相互之间的材料渗透，也未对接触载荷大小做任何限制的问题。

(4)本申请的技术方案开发几何精确壳单元建立复合材料壁板的有限元模型，计入连接区开孔，准确预测孔边应力分布，解决技术方案一中无法预测孔边应力分布的问题。

(5)本申请的技术方案采用计入翘曲变形的几何精确梁单元建立螺栓的有限元模型，与上述其他改进措施一起，解决现有技术方案一计算精度低，现有技术方案二计算效率低的问题。

附图说明

图1、本发明流程框架图；

图2、四边固支平板在剪切载荷作用下的变形示意图；

图3a、复合材料平板1阶模态；

图3b、复合材料平板2阶模态；

图3c、复合材料平板3阶模态；

图3d、复合材料平板4阶模态；

图3e、复合材料平板5阶模态；

图3f、复合材料平板6阶模态；

图4、接触模型的运动学描述示意图；

图5a、复合材料板单剪试验件俯视示意图；

图5b、复合材料板单剪试验件正视示意图；

图6a、hertz接触模型复合材料壁板单剪试验件变形分布；

图6b、mpc模型复合材料壁板单剪试验件变形分布。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细描述。应当理解，此处所描述的具体实施例仅仅用于解释本发明，并不用于限定本发明。

相反，本发明涵盖任何由权利要求定义的在本发明的精髓和范围上做的替代、修改、等效方法以及方案。进一步，为了使公众对本发明有更好的了解，在下文对本发明的细节描述中，详尽描述了一些特定的细节部分。对本领域技术人员来说没有这些细节部分的描述也可以完全理解本发明。

实施例1

(一)对复合材料壁板构建壳单元：

根据复合材料壁板外形尺寸、材料属性等构建复合材料壁板有限元离散模型，与螺栓相比，复合材料壁板柔性更大、变形明显。在外载荷作用下，复合材料壁板的变形具有非线性特征，具体的理论为：

几何精确壳单元理论

本申请将定义三种不同构型描述壳单元在时间域内的非线性变形。当t＝0时为参考构型，此时壳单元无变形；当t＝ti时为初始构型，作为当前步骤变形计算的初始输入；当t＝tf时为变形构型，为未知构型，需要在当前步骤的计算中确定。

在参考坐标系b中给定壳的厚度h和中面面积ω，以及壳中面上任意材料点的位置坐标x0和沿中面法向方向的材料坐标ξ。壳中任意材料点的位置坐标x可写为

x(α1，α2，ξ)＝x0(α1，α2)+s(α1，α2，ξ)(1)

其中，α1和α2是壳中面的材料点坐标，s(α1，α2，ξ)是垂直于中面的向量。坐标α1,α2和ξ称之为曲线坐标。

相应的，变形构型中材料点的位置矢量可写为

x(α1，α2，ξ)＝x0(α1，α2)+s(α1，α2，ξ)＝x0(α1，α2)+u(α1，α2)+s(α1，α2，ξ)(2)

式中，x0(α1，α2)壳中面材料点位置坐标，u(α1，α2)壳中面位移矢量。几何精确壳单元理论推导中，需要用到如下reissner-midlin假设：与壳中面垂直的材料线变形后仍然保持直线，不发生拉伸变形。在该假设中未约束壳变形后材料线仍然垂直中面，因此壳会发生剪切变形。

依据弹性力学理论中green-lagrange应变张量的定义，通过高阶项的近似处理，可得适合于任意大位移、有限转动和大变形情况的新型应变-位移关系

上述表达式中，和为曲率应变的二次高阶项，上标i代表当前计算步骤的初始值。即上述应变-位移关系式采用高阶曲率应变的初始值近似当前计算步骤的应变值。该表达式为本次申请首次提出，在简化了壳的应变-位移表达式同时又计入了高阶应变项，具有处理大变形的能力,同时又提高了计算效率。壳的运动控制方程由hamilton原理推导。在时间域内通过对任意材料点的位置矢量，方程(2)，在时间域内进行求导，可得速度矢量

由此，获得壳的动能为

动能变分可写为

式中h和g分别代表线性动量矩矢量和角动量矩矢量。

壳的应变能变分是

其中τ^*ij，i＝1，2，j＝1，2，3，为随体坐标系中描述的cauchy应力张量。将应变分量带入应变能变分表达式并沿厚度积分可得应变能变分的等价形式

外力虚功可表示为

其中f和m代表壳中面单位面积上作用的力和力矩。将动能变分、应变能变分和外力虚功带入hamilton原理，最终得到壳的运动控制方程

式中f和m代表壳中面单位面积上作用的外部力和力矩，h和g代表动量矩，ni，j,n3和mi，j，i＝1，2，为弹性力和力矩，b为转换矩阵，为单位矢量的旋转对称矩阵，rr0为旋转张量。忽略上述方程中的惯性力，即可获得静力学变形控制方程。依据上述推导出的公式并结合有限元理论，本申请的技术方案完成了9节点二阶精度壳单元的开发，该单元具有静力学、模态分解和动力学分析功能。

为了确保所开发壳单元的精度，针对该单元进行了静力学和模态分析测试。给定薄壁矩形板(h/l＝1/1000)宽度为a＝1.4英寸，厚度为h＝2.4×10^-3英寸。考虑固支边界条件和均布剪切压力，分别采用解析法、当前开发的几何精确壳单元和商用软件abaqus三种不同的方式预测了薄壁矩形板的变形。该矩形板的材料属性设置为：杨氏模量e＝2×10⁶帕，泊松比ν＝0.3。根据kirchhoff平板理论可以获得该矩形板质心变形(也即最大变形)的解析解如下

上述公式中，d为抗弯刚度。采用几何精确壳单元求解时，使用了18×18个本申请二阶精度壳单元。在质心集中剪切载荷p0＝1×10^-4磅作用下四边固支平板的面外变形情况参见图2。

理论分析和数值计算表明，所开发的几何精确壳单元能够准确的预测薄壁矩形板的面外变形。该计算结果与理论解、abaqus数值解均进行了比较，结果见表1。

表1矩形板的面外变形

模态分析测试中，建立了复合材料平板的几何精确有限元模型。复材板的尺寸为：长272毫米，宽48毫米，厚度5.2毫米。铺层数目为40层，铺层厚度为0.13毫米，铺层方向依次为[45,0,-45,90]。复合材料板的材料选用为碳纤维环氧树脂hta/6376，属性设置为：杨氏模量exx＝54.25兆帕，eyy＝54.25兆帕，ezz＝12.59兆帕，剪切模量gxy＝20.72兆帕，gxz＝4.55兆帕，gyz＝4.55兆帕，泊松比νxy＝0.309，νxz＝0.332，νyz＝0.332。采用本申请的壳单元技术、该复合材料板离散成34×6个本申请二阶精度壳单元。完成模态分析后，计算出前6阶结构模态，如图3a-图3f所示。

该计算结果和商用软件abaqus的计算结果进行了对比分析，参见表2。结果表明，所开发壳单元的模态分析具有足够的计算精度。

表2复合材料板前6阶模态

。对螺栓构建梁单元数值模型：

(二)根据螺栓外形尺寸、材料属性等构建螺栓有限元数值模型，相对于复合材料壁板，由航空级钛合金制作的螺栓刚性更大。在复合材料壁板螺栓连接件承受拉伸/压缩载荷时，螺栓会产生“跷跷板”式的旋转，螺栓截面变形可以忽略。因此，有限元建模中，可以采用一维梁单元模拟螺栓“跷跷板”式的旋转，具体的理论为：

给定梁参考线上任意一个点的位置向量，u0(α1)，梁截面上任意点的位置矢量可写为

x＝u0(α1)+s(α1，α2，α3)(12)

上式中α1,α2,α3为曲线坐标。相应的，变形后该材料点的位置矢量变为

x＝u0(α1)+u(α1)+s(α1，α2，α3)(13)

其中，u(α1)为梁参考线上任意点的位移。梁单元的应变-位移关系仍然根据green-lagrange应变定义推导获得，并能够处理大变形和任意转动。为了便于处理转动,将梁单元的应变-位移关系投影到随体坐标系中。相应的一维应变矢量可写为

其中为曲率应变矢量。另外，应变矢量包含轴向应变和切向应变

该表达式中rf为旋转张量。梁单元建模仍然采用hamilton原理，最终得到梁单元半离散形式的运动控制方程：

式中f和t代表梁单元单位展长上作用的外部力和力矩，h和g代表动量矩，n1和m1为弹性力和力矩，ω为试函数，u0和u为梁截面上任一点的初始位置和当前位移。

(三)接触单元：模拟螺栓和复合材料壁板孔边的真实接触，并预测接触时复材壁板和螺栓分别承受的接触载荷，具体理论为：

基于hertz接触模型的连接建模理论

本申请的技术方案从复合材料壁板螺栓连接设计的实际需求出发，利用所开发的几何精确壳单元，兼顾精度和效率，进一步发展了适用于螺栓连接快速分析的接触模型。该连接建模引入了hertz模型，能够准确预测接触区的微变形和接触反力，进一步确定带孔复合材料壁板的应力分布以及钉载分配比例。

复合材料壁板螺栓连接与常规接触问题不同，螺栓的刚度通常要大于复合材料壁板的刚度。螺栓预紧力的存在使得螺帽与复合材料壁板之间的接触区必须满足应变协调条件。本申请的技术方案为了提高计算效率，引入了如下假设：

(1)螺帽与复合材料壁板之间的接触区不发生面外位移；

(2)复合材料壁板采用几何精确壳单元建模，忽略厚度方向变形；

(3)变形前后，螺帽和复合材料壁板之间的接触区与复合材料壁板中面保持平行。

基于上述假设，接触模型可简化为局部的二维平面接触，参见图4。如图所4示，基于hertz接触模型的螺栓连接建模，在螺栓随体坐标系bb中引入了局部位移矢量us。假设在参考构型，螺栓与复合材料壁板孔边接触但没有接触力。定义从螺栓随体坐标系原点到该接触点的矢量为u0。随后，在变形构型中发生接触，进而产生相互渗透。此时，由于初始接触位置移动而产生的位移矢量为us。为了便于描述由于接触产生的渗透，引入新的位置矢量

接触产生的渗透可由下式计算

q＝|ur|-|u0|，a＝-q(18)

因为接触而产生的接触力，沿着螺栓接触面内法线方向作用在螺栓上，沿内法线反方向作用在复合材料壁板上，阻止二者之间的进一步接触。内法线定义为

根据非线性hertz接触理论，当q＜0时法向接触力以及法向接触力虚功可写为

式中，k为接触刚度，可根据发生接触的两个物体的真实材料属性严格确定。为矢量us在螺栓随体坐标系中的映射。根据虚功原理，可以得到接触载荷矢量为

当q≥0时，接触载荷矢量fc＝0。在数值计算中，时间域内的每一积分步都需要判断接触是否发生。

利用本申请提供的方案分析由三排钉连接的复合材料壁板单剪试验件的传载情况。图5a、5b给出了该单剪试验件的几何尺寸。

该单剪试验件包含上下复合材料平板以及三个螺栓。复合材料板的铺层以及材料属性可参考用于模态分析的复材平板参数设置。该复合材料板离散成几何精确壳单元。螺栓直径为8毫米，材料为钛合金，材料属性：杨氏模量e＝110兆帕，泊松比ν＝0.29。三个螺栓均由几何精确梁单元模拟。复合材料平板和螺栓之间的连接采用非线性hertz接触模型。

为了验证hertz接触模型的有效性，该单剪试验件同时采用多点约束近似模拟了连接作为对比。计算结果与本申请方法考虑接触的计算结果进行了比较。图6a、图6b给出了两种不同模型在相同拉伸载荷作用下的变形。由于单剪连接造成该连接件传载不共线，在拉伸载荷作用下，试验件会产生“二次弯曲”现象。图6a、图6b中两种模型的变形均捕捉到了这一现象，说明变形预测结果是正确的。获得复合材料平板变形分布后，根据线弹性材料本构关系，可获得带孔复合材料平板的应力分布用于强度校核。

另外，本申请的技术方案比较了上述两种模型的螺栓传载比例，见表3。结果表明，数值模型计入接触模型后，接近真实物理模型，传载比例更接近于均匀分布。

表3单剪试验件螺栓传载比例