CN113705047A - 一种拉伸载荷下cfrtp干涉连接孔周应力分布计算方法 - Google Patents

Abstract

本发明提供了一种拉伸载荷下CFRTP干涉连接孔周应力分布计算方法，建立CFRTP的本构模型和损伤准则，建立CFRTP干涉链接模型实现数值分析，将所得结果拟合成函数形式的孔周位移约束条件；然后，将问题分解为干涉状态和拉伸状态两部分，将所得位移约束条件进行孔周应力求解，并建立在载荷传递模型计算拉伸载荷和位移。本申请提出有限元数值建模与解析计算相结合的方法，精确预测CFRTP的孔周应力分布。通过VUMAT用户子程序实现CFRTP干涉连接拉伸的数值分析，获得孔周位移边界条件的精确数值解，并将其拟合为三角函数的形式。将问题分解为干涉状态和拉伸状态，进行求解。建立载荷传递模型关联孔周应力和拉伸力，并进行CFRTP单搭接件的拉伸试验与计算结果对比验证其正确性。

Description

一种拉伸载荷下CFRTP干涉连接孔周应力分布计算方法

技术领域

本发明涉及碳纤维增强树脂基复合材料(CFRTP)结构设计和寿命预测技术领域，特别地涉及一种拉伸载荷下CFRTP干涉连接孔周应力分布计算方法。

背景技术

碳纤维增强树脂基复合材料(CFRTP)具有比强度高、比模量高、疲劳性能好、耐腐蚀等特点，被广泛应用于航空航天领域。当CFRTP作为飞机结构件时，不可避免要进行连接。而目前飞机装配中常用的螺栓连接会使螺栓杆与孔壁产生干涉，挤压孔壁导致应力分布不均匀造成结构损伤，给CFRTP结构设计和寿命预测带来了很多问题。因此，CFRTP干涉连接孔周应力分布是解决这些问题基础。

发明内容

针对上述现有技术中，CFRTP干涉连接孔周应力分布计算效率和计算准确度的问题，本申请提出了一种拉伸载荷下CFRTP干涉连接孔周应力分布计算方法，建立CFRTP的本构模型和损伤准则，建立CFRTP干涉链接模型实现数值分析，将所得结果拟合成函数形式的孔周位移约束条件；然后，将问题分解为干涉状态和拉伸状态两部分，将所得位移约束条件进行孔周应力求解，并建立在载荷传递模型计算拉伸载荷和位移。

本申请提出有限元数值建模与解析计算相结合的方法，精确预测CFRTP的孔周应力分布。通过VUMAT用户子程序实现CFRTP干涉连接拉伸的数值分析，获得孔周位移边界条件的精确数值解，并将其拟合为三角函数的形式。将问题分解为干涉状态和拉伸状态，进行求解。建立载荷传递模型关联孔周应力和拉伸力，并进行CFRTP单搭接件的拉伸试验与计算结果对比验证其正确性。

优选的，CFRTP的本构模型包括单向CFRTP单层板的应力应变关系；损伤准则包括定义损伤起始判据。

优选的，单向CFRTP单层板的应力应变关系为：

优选的，损伤起始判据包括：

纤维拉伸损伤(σ₁≥0)

纤维压缩损伤(σ₁<0)

面内基体压缩损伤(σ₂≥0)

面内基体拉伸损伤(σ₂<0)

面外基体压缩损伤(σ₃≥0)

面外基体拉伸损伤(σ₃<0)

中的至少一项；

式中，σ₁,σ₂,σ₃,τ₁₂,τ₁₃,τ₂₃代表各方向的正应力的剪应力，X_T,X_C表示纵向抗拉/抗压强度，Y_T,Y_C表示横向抗拉/抗压强度，Z_T，Z_C表示面外抗拉/抗压强度，S₁₂,S₂₃,S₁₃表示对应方向的剪切强度。

上述技术特征可以各种适合的方式组合或由等效的技术特征来替代，只要能够达到本发明的目的。

本发明提供的一种拉伸载荷下CFRTP干涉连接孔周应力分布计算方法，与现有技术相比，至少具备有以下有益效果：在求解CFRTP孔周应力分布方面，数值与解析结合法具有较高的计算效率和计算准确度。

附图说明

在下文中将基于实施例并参考附图来对本发明进行更详细的描述。其中：

图1显示了本发明干涉连接数值分析模型图；

图2显示了本发明1500N拉伸载荷孔周位移约束曲线图；

图3显示了本发明拉伸载荷下干涉连接受力示意图；

图4显示了本发明载荷传递模型图；

图5显示了本发明单搭接拉伸试验示意图；

图6显示了本发明载荷位移曲线与DIC实验结果图；

图7显示了本发明孔周应力分布图。

具体实施方式

下面将结合附图对本发明作进一步说明。

本发明的一种拉伸载荷下CFRTP干涉连接孔周应力分布计算方法，提出有限元数值建模与解析计算相结合的方法，精确预测CFRTP的孔周应力分布。通过VUMAT用户子程序实现CFRTP干涉连接拉伸的数值分析，获得孔周位移边界条件的精确数值解，并将其拟合为三角函数的形式。将问题分解为干涉状态和拉伸状态，进行求解。建立载荷传递模型关联孔周应力和拉伸力，并进行CFRTP单搭接件的拉伸试验与计算结果对比验证其正确性。

1位移边界条件的数值分析；

1.1CFRTP本构模型和渐进损伤准则；

单向CFRTP单层板的应力应变关系可以写为：

式中Cij可用CFRTP各方向的弹性模量、剪切模量、泊松比计算得到。由于CFRTP干涉连接的干涉量不易不超过1.5％，层间损伤较少。定义损伤起始判据为：

纤维拉伸损伤(σ₁≥0)

纤维压缩损伤(σ₁<0)

面内基体压缩损伤(σ₂≥0)

面内基体拉伸损伤(σ₂<0)

面外基体压缩损伤(σ₃≥0)

面外基体拉伸损伤(σ₃<0)

式中，σ₁,σ₂,σ₃,τ₁₂,τ₁₃,τ₂₃代表各方向的正应力的剪应力，X_T,X_C表示纵向抗拉/抗压强度，Y_T,Y_C表示横向抗拉/抗压强度，Z_T，Z_C表示面外抗拉/抗压强度，S₁₂,S₂₃,S₁₃表示对应方向的剪切强度。

刚度退化矩阵定义为：

式中，G_fc,G_ft,G_mc,G_mt,G_zc,G_zt分别为纤维压缩/拉伸、面内基体压缩/拉伸、面外基体压缩/拉伸的临界应变能释放率。L^C是单元的特征长度，通过VUMAT中“CharLength”控制的初始特征单元长度L_initial，以及单元在全局坐标下的尺寸l_z控制。

1.2CFRTP螺栓连接数值分析：

如图1所示，为了获得拉伸载荷下准确的孔周位移约束条件，在ABAQUS中建立铺层为[0]_18s的CFRTP干涉连接的模型，并将上述本构模型嵌入用户定义材料属性子程序VUMAT中实现数值分析。本申请的关注重点是螺栓杆和CFRT P孔壁的干涉区域，因此将螺栓外形简化，并在接触区域划分0.1mm×0.1mm×0.1mm的C3D8R网格。螺栓中心处易出现单元扭曲变形导致计算不收敛，因此在中心处划分更精细的0.05mm×0.05mm×0.05mm的C3D8R网格。螺栓与CFRTP孔壁采用小滑移接触，定义螺栓靠近CFRTP板下表面的横截面为预加载面，采用ABAQUS的“bolt load”功能实现螺栓的预紧力加载。设置拉伸载荷范围300N～3700N，分为十组，分别加载在干涉连接模型上。CFRTP板的材料参数见表1。

表1 CFRTP材料参数

采用MATLAB中的“fittype”函数，将十组仿真数据拟合成数值边界条件表达式：

式中，P是拉伸载荷，E₁，E₂为对应方向的弹性模量，h为约束平面距加载面的距离，θ是孔周沿X方向的角度。如图2所示，取拉伸载荷为1500N，绘制孔周位移约束曲线。

2基于复函数理论孔周应力建模

以仿真结果作为位移边界条件约束，采用解析法分析孔周应力。如图3所示，参考载荷分析方法，将变形区划分为干涉配合作用和拉伸载荷作用两部分的叠加。对于受载形式未知，仅知合力大小的正交各向异性材料孔，其孔周平面应力分量的一般形式为：

式中，

是特定载荷下的应力函数，复变量z_k＝x+μ_ky(k＝1,2)，μ_k是式(15)的虚根。式(15)通过平衡方程、应变协调方程、材料本构关系联立求得。

因此，求得孔周应力的关键在于求解应力函数。

在式(13)形式下的位移边界条件基础上，拉伸载荷下孔周应力函数可表示为式(16)的形式，

式中，W₁,W₂,V₁,V₂由式(13)求得。A,B,D,p₁,p₂,q₁,q₂与材料参数相关，由式(17)，式(18)求得。ζ_k(k＝1,2)是保角映射得到的表达式，由式(19)求得。

上述三式中a_ij为复合材料的柔度系数。r₀是复合材料孔的半径。

求解正交各向异性材料干涉配合下的应力函数，

式中Δ是平均干涉量。因此，拉伸载荷下的干涉连接应力函数可写为：

对式(21)求导，

将式(22)代入式(14)，即可求得拉伸载荷下的孔周应力分布。

如图4所示，为了获得拉伸载荷和孔周应力分布的关系，将板划分为三个区域，其中区域2上下部分为载荷传递区域B，孔周两侧区域分别为A₁,A₂。将各部分视为等效弹簧，不考虑螺栓变形，将其视为固支连接的不可压缩杆。那么，可以建立以下载荷传递模型。

式中，F_A2，F_A1，F_B为各区域传递的载荷，K_i＝E_xxw_it_i/l_i为等效弹簧的刚度系数，E_xx表示沿长度方向的弹性模量，w_i，t_i，l_i分别是对应区域的等效宽度，厚度和等效长度。通过式(24)，可以得出以下结论：

其中，Σ表示A₂挤压区域。

3CFRTP单搭接拉伸实验；

为了验证所得孔周应力分布的正确性，按照ASTM D 5961/5961M-05标准对铺层为[0]_18s的碳纤维增强环氧树脂基复合材料进行单搭接拉伸试验，并对孔周变形区进行DIC测量。所用材料为光威公司提供的单向增强碳纤维增强环氧树脂基复合材料，其碳纤维含量为60％，铺层为[0]_18s，具体力学性能见表1。所用螺栓为直径6mm的钛合金螺栓。通过改变拧紧力矩控制干涉量为0.2％。具体实验布置见图5，采用CTM2050万能试验机进行拉伸试验，控制拉伸速度1mm/min准静态加载，达到预定载荷后停止，连接计算机记录拉力和位移。在拉伸件的孔周随机标定散斑，采用Phantom-v212高速相机拍摄拉伸全过程，经软件处理后得到孔周应变分布。

结果与讨论：

由式(25)得出理论载荷位移曲线，与实验结果对比如图6所示。理论结果的数值与实验结果基本吻合，最大误差14％。理论计算的载荷位移曲线为线性增长，而实验结果表明在拉伸初始阶段载荷增长迅速，在900N附近由于材料的损伤导致载荷下降，随后继续线性增长。DIC测试结果表明1方向与2方向的正应变在挤压区域较大，而剪应变在90°方向较大。将DIC测试得到的应变分布结果通过CFRTP柔度矩阵转化为应力分布，再通过坐标变换公式转为极坐标下的应力分布。提取仿真结果中试件表面的应力分布，也通过坐标变换公式转为极坐标下的应力分布。

如图7所示，径向应力在0°到90°区域均为负值，与钉孔接触挤压区域状态相一致。在0°取最大值约为110MPa，在90°时取最小值约为20MPa。与有限元结果相比，数值与解析结合法得到的结果更接近于实验值，但仍存在11％左右的误差。这是因为螺栓遮盖了0.5mm的孔周区域，导致DIC测量有误差，另一方面是因为计算过程忽略了钛合金螺栓的变形，将其做刚体处理导致径向应力计算结果偏大。周向应力在0°到60°区域基本保持不变，之后突然增长，在90°达到最大值300MPa左右。与有限元结果相比，数值与解析结合法得到的结果更接近于实验值。由此可见，数值与解析结合法在复合材料结构力学分析方面有着准确、快速的优点，可为结构设计提供更可靠的力学依据。

结论：

1)通过对螺栓连接过程建模与数值分析，发现CFRTP干涉连接孔周位移约束的可采用三角函数形式的关系描述，并与拉伸载荷、弹性模量、应力平面所在位置有关。

2)数值与解析结合法得到的载荷位移曲线和应力分布与DIC实验结果吻合较好，存在14％和11％的误差。

3)在求解CFRTP孔周应力分布方面，数值与解析结合法较数值分析法提高了计算效率和计算准确度。

4)3700N拉伸载荷下，干涉量0.2％铺层[0]_18sCFRTP螺栓连接件的径向应力从0°到90°均匀增长，最大值110MPa出现在0°，周向应力在0°到60°区域基本保持不变，之后突然增长，最大值300MPa出现在90°。

虽然在本申请中参照了特定的实施方式来描述本发明，但是应该理解的是，这些实施例仅仅是本发明的原理和应用的示例。因此应该理解的是，可以对示例性的实施例进行许多修改，并且可以设计出其他的布置，只要不偏离所附权利要求所限定的本发明的精神和范围。应该理解的是，可以通过不同于原始权利要求所描述的方式来结合不同的从属权利要求和本申请中所述的特征。还可以理解的是，结合单独实施例所描述的特征可以使用在其他所述实施例中。

Claims (4)

1.一种拉伸载荷下CFRTP干涉连接孔周应力分布计算方法，其特征在于，建立CFRTP的本构模型和损伤准则，建立CFRTP干涉连接模型实现数值分析，将所得结果拟合成函数形式的孔周位移约束条件；然后，将问题分解为干涉状态和拉伸状态两部分，将所得位移约束条件进行孔周应力求解；最后，在载荷传递函数的基础上，将孔周应力与拉力相关联。

2.根据权利要求1所述的一种拉伸载荷下CFRTP干涉连接孔周应力分布计算方法，其特征在于，所述CFRTP的本构模型包括单向CFRTP单层板的应力应变关系；损伤准则包括定义损伤起始判据。

3.根据权利要求1-2任一项所述的一种拉伸载荷下CFRTP干涉连接孔周应力分布计算方法，其特征在于，单向CFRTP单层板的应力应变关系为：

式中C_ij可用CFRTP各方向的弹性模量、剪切模量、泊松比计算得到。

4.根据权利要求1-3任一项所述的一种拉伸载荷下CFRTP干涉连接孔周应力分布计算方法，其特征在于，损伤起始判据包括：

纤维拉伸损伤(σ₁≥0)

纤维压缩损伤(σ₁<0)

面内基体压缩损伤(σ₂≥0)

面内基体拉伸损伤(σ₂<0)

面外基体压缩损伤(σ₃≥0)

面外基体拉伸损伤(σ₃<0)

中的至少一项；

式中，σ₁,σ₂,σ₃,τ₁₂,τ₁₃,τ₂₃代表各方向的正应力的剪应力，X_T,X_C表示纵向抗拉/抗压强度，Y_T,Y_C表示横向抗拉/抗压强度，Z_T，Z_C表示面外抗拉/抗压强度，S₁₂,S₂₃,S₁₃表示对应方向的剪切强度。

Images