2. 郑州航空工业管理学院 无人机研究院, 河南 郑州 450046;
3. 西北工业大学 航空学院, 陕西 西安 710072
金属材料广泛应用于工程领域中，如航空、航海和汽车等领域。疲劳失效是金属构件最常见的一种失效形式，美国土木工程师协会(ASCE)研究表明金属构件80%~90%的失效都与疲劳有关[1]。随着现代工程结构对总体性能要求的不断提高，结构构件不断向精细化设计和整体化制造发展，这使结构疲劳问题越来越突出。另外，疲劳寿命存在很大的分散性，尤其是对于高周疲劳，使得在结构可靠性安全寿命设计过程中存在很大的难题。工程上对于疲劳寿命预测的研究主要是通过标准件疲劳实验建立应力、应变及能量和确定性寿命(中值寿命或某一可靠度下的寿命)的关系，如S-N曲线、ε-N曲线和Darvaeux模型等。确定性寿命不能充分反映疲劳寿命的概率统计特征，如概率分布类型及标准差、峰度和偏度等，概率疲劳寿命应能表征其概率分布和相应的分布参数。目前，多数研究是基于这些模型的改进，如白恩军等[2]将不同级应力水平的疲劳寿命等效到最高级应力水平，实现了小样本下p-S-N曲线的拟合。刘红彬等[3]通过改进的平均应力公式预测了涡轮盘螺栓孔处的疲劳寿命。罗媛等[4]基于S-N曲线与Miner准则，构建了桥梁疲劳损伤的功能函数，用于简支梁桥的疲劳可靠度评估。马玉娥等[5]通过玻璃纤维铝合金层板(FMLs)疲劳实验建立了S-N曲线，为其疲劳寿命预测提供了材料基础性能和信息。Doudard等[6]在最弱连接理论框架下提出一种概率双尺度模型，用于预测S-N曲线的分散性。
上述方法相对简单易用，但对于变幅加载下的复杂工程构件，需结合损伤累积准则预测疲劳寿命[7]，且预测结果通常偏差较大，在工程应用上有一定局限性[8]。近年来，疲劳损伤形成和演化机理成为国内外学者的研究重点之一，连续介质损伤力学(continuum damage mechanics, CDM)[9]在此背景下逐渐发展成熟。Lemaitre等[10]基于损伤力学发展出一套疲劳损伤演化模型用于确定性疲劳寿命预测，但目前仍鲜见有基于损伤力学的概率疲劳寿命预测方法。本文旨在采用金属材料双尺度局域化准则和严谨的连续介质损伤力学理论，将Lemaitre疲劳损伤演化模型参数概率化，针对金属构件，建立预测高周疲劳寿命的概率方法。
1 双尺度疲劳损伤演化概率模型由于高周疲劳宏观无塑性，本节将首先建立宏细观双尺度桥联方程，用于计算细观应力应变场，作为Lemaitre疲劳损伤演化模型的输入，并将模型概率化，建立预测金属构件高周概率疲劳寿命的算法流程。
1.1 宏细观双尺度桥联方程基于Eshelby-Kroner局域化准则和损伤演化模型，Lemaitre等[10]于1999年提出了一种双尺度疲劳损伤演化模型，可以有效预测金属材料高周疲劳寿命。本节将对其进行简要介绍。
由于材料中存在各种缺陷, 疲劳损伤演化的微观机理极其复杂, 而微观尺度下大量微观缺陷的力学特性无法直接表征。微观缺陷总是劣化或减弱材料的力学性能, 从均匀化的力学角度出发, Eshelby-Kroner局域化准则将这些缺陷统一看作细观尺度下RVE内的一个弱夹杂, 再通过弱夹杂的力学行为分析损伤演化, 如图 1所示, 其中上标μ表示细观尺度下的变量, 后文均代表此意义。
|
| 图 1 Eshelby-Kroner局域化准则 |
其中, 联系宏细观应力应变的方程如(1)式所示
|
(1) |
(1) 式也称为Eshelby-Kroner局域化方程[11], 式中:
|
(2) |
结合Eshelby-Kroner局域化方程和经典弹塑性本构方程, 可得求解细观应力应变的方程组, 其矢量表达如(3)式所示
|
(3) |
式中: σμ, εμe和εμp是未知的细观应力、弹性应变和塑性应变; fμ是屈服函数; Xμ是随动强化背应力; Rμ是由于等向强化引起的屈服面半径的增量; Dμe是弹性模量; σ, ε和εp是已知的宏观应力、总应变和塑性应变。
对于实际工程构件, 需首先通过有限元分析其在载荷谱作用下的危险点宏观应力应变场序列, 作为(3)式的初始输入。为求解细观应力应变场序列, 由于其细观进入塑性, 因此需采用隐式迭代求解算法, 包括弹性预测和塑性修正, 算法流程如图 2所示。
|
| 图 2 宏细观双尺度桥联模型隐格式求解算法流程框图 |
基于细观应力应变场序列, 可代入Lemaitre疲劳损伤演化模型计算高周疲劳损伤累积, 如下
|
(4) |
式中: Ḋ为损伤演化率; pμ为细观等效塑性应变; pD为损伤累积门槛值; Dc为疲劳损伤临界值; S和s为损伤演化参数; Yμ为损伤能量释放率, 表达式如(5)式
|
(5) |
式中: σeqμ是von Mises等效应力; σmμ是平均球应力; E是弹性模量; ν是泊松比; Rvμ应力三轴度函数。
本文将按照Lemaitre文献对模型进行一些简化[10]: 忽略损伤和弹塑性本构方程的耦合; 疲劳损伤门槛值pD=0, 意味着进入塑性即发生损伤; 疲劳损伤临界值Dc=1;仅考虑材料随动强化, 且细观随动强化参数与宏观一致。模型参数E, ν, σ和C可通过材料静力拉伸实验或查手册获得, 对于金属材料, 局域化方程参数α通常取为0.4, 此时仅剩S和s未知, 因此疲劳寿命N是关于S和s的函数
|
(6) |
式中: 上标“+”代表上半拉伸循环；上标“-”代表下半压缩循环。
通常情况下, 模型参数可通过模型反演优化获得, 通过使模型预测寿命和实验寿命平方和最小, 建立如下优化方程
|
(7) |
式中: Nt为疲劳实验寿命; S和s的取值范围参考文献[12], 每个实验寿命可得到1组S和s的样本值。
由于相同应力水平加载下, 疲劳寿命具有很大分散性, 因此本文认为S和s是服从某种分布的随机变量, 通过(7)式优化获得S和s样本, 首先利用EDF拟合优度检验其分布类型, 并统计样本获得其概率特性。对于金属构件承受复杂加载情况, 如多级块谱或随机谱加载, 本文方法预测其高周疲劳寿命分为以下3个步骤:
1) 通过有限元计算金属构件危险点的宏观应力应变。对于高周疲劳, 宏观应力应变处于弹性状态, 因此在弹性状态下进行一次有限元分析, 确定构件危险细节并提取其应力应变, 乘以相对于载荷谱峰谷值比例系数即可获得应力应变峰谷值序列;
2) 按双尺度桥联模型隐格式算法进行细观应力应变滞回曲线的非线性计算, 得到相应的细观应力应变峰谷值序列;
3) 对S和s进行随机抽样, 每组样本代表 1个“试件”, 通过(4)式计算构件危险点在每个载荷循环的疲劳损伤量, 并累积损伤, 记录循环数; 当D≥Dc时, 结束计算, 得到总循环数, 即疲劳寿命N的1个样本。
上述双尺度疲劳损伤演化算法通过Fortran编程实现, 其中包括宏细观双尺度应力应变计算以及疲劳损伤累积计算, 作为非嵌入后处理器来计算每个循环内累积的损伤并预测疲劳寿命。
2 实验设计及模型验证本节设计完成了铝合金光滑圆棒在6组应力水平下的高周疲劳实验, 选取实验材料为工程中常用的铝合金2024-T3, 主要目的是通过铝合金2024-T3疲劳实验数据分别拟合得到细观尺度下损伤演化模型中的参数S和s。之后将利用双尺度疲劳损伤演化模型预测飞机翼身连接壁板疲劳寿命, 并与实验结果对比。
2.1 模型参数获取铝合金光滑圆棒实验件如图 3所示, 实验在6组不同应力水平下进行, 应力比均为-1, 每组8个实验件, 实验机为Instron8802电液伺服万能实验机, 按照《HB 5287-96金属材料轴向加载疲劳实验方法》进行。实验结果如表 1所示。
|
| 图 3 铝合金光滑圆棒疲劳实验件 |
| Sa/MPa | 寿命/周次 |
| 360 | 5 503, 7 355, 7 983, 6 016, 7 286, 5 298, 7 706, 6 161 |
| 315 | 14 542, 19 515, 12 181, 17 905, 14 202, 17 721, 14 857, 17 685 |
| 270 | 36 770, 39 729, 32 528, 53 710, 26 296, 35 577, 43 863, 26 752 |
| 225 | 103 802, 159 956, 149 417, 113 396, 188 604, 167 124, 129 531, 150 157 |
| 200 | 775 796, 421 346, 190 824, 253 911, 389 915, 204 056, 430 238, 241 385 |
| 180 | 2 004 958, 724 338, 905 673, 3 512 550, 1 551 041, 781 714, 1 921 325, 1 244 010 |
由于实验为光滑圆棒单轴拉伸, 且6级应力水平下试件未进入塑性, 因此可通过理论计算获得各应力水平下宏观应力应变数据, 作为双尺度桥联模型的输入, 计算细观应力应变场, 进而通过(7)式反演优化获得S和s样本, 每级应力水平下将得到8组S和s样本, 对其进行关于正态分布的EDF检验[13](95%置信度下), 并统计其均值和方差, 结果如表 2所示, 表中Dn, W2和A2为EDF统计量, 具体定义可参考文献[13]。由表可知, S和s均通过关于正态分布的EDF检验, 且不同应力水平下, S和s的均值和标准差变化很小, 因此本文认为其概率特性不随外加载荷变化, 统计6组应力水平下的S和s样本, 可得S均值和标准差分别为3.050 6和0.307 0, s的均值和标准差分别为2.353 1和0.159 9。
| Sa/MPa | 参数 | μ | σ | Dn | W2 | A2 |
| 360 | S | 3.063 4 | 0.287 4 | 0.396 7√ | 0.028 0√ | 0.201 1√ |
| s | 2.305 8 | 0.154 3 | 0.459 9√ | 0.044 8√ | 0.275 1√ | |
| 315 | S | 3.085 7 | 0.294 3 | 0.587 0√ | 0.046 5√ | 0.262 6√ |
| s | 2.387 4 | 0.157 7 | 0.548 8√ | 0.059 0√ | 0.392 3√ | |
| 270 | S | 3.047 1 | 0.301 8 | 0.556 0√ | 0.102 7√ | 0.577 8√ |
| s | 2.325 7 | 0.159 4 | 0.798 4√ | 0.107 9√ | 0.627 8√ | |
| 225 | S | 3.059 4 | 0.314 6 | 0.785 0√ | 0.111 8√ | 0.672 0√ |
| s | 2.358 9 | 0.162 9 | 0.679 7√ | 0.067 7√ | 0.385 0√ | |
| 200 | S | 3.036 9 | 0.320 8 | 0.429 6√ | 0.034 8√ | 0.257 4√ |
| s | 2.371 6 | 0.163 6 | 0.486 5√ | 0.032 3√ | 0.203 6√ | |
| 180 | S | 3.011 1 | 0.323 4 | 0.533 1√ | 0.055 6√ | 0.331 0√ |
| s | 2.369 1 | 0.161 7 | 0.511 7√ | 0.044 9√ | 0.265 5√ | |
| 注：“√”表示在95%置信度下, 通过关于正态分布EDF检验 | ||||||
飞机连接壁板取自某型客机, 为蒙皮-带板搭接壁板, 通过螺栓连接, 几何形状及螺孔位置分布如图 4所示, 图中未给出螺栓具体构型。模型关于中面对称。沿对称面从内到外的三排连接螺孔的直径分别为12.7, 12.7和11.11 mm。
|
| 图 4 蒙皮-带板搭接壁板 |
蒙皮和带板材料为铝合金2024-T3, 螺栓材料为30CrMnSiA。蒙皮-带板壁板常幅谱最大应力(典型截面应力)为80 MPa, 应力比R均为0.06。
2.2.1 搭接壁板有限元分析按照图 2所示的算法流程图, 本文首先对壁板构件进行有限元分析, 获取构件危险点的应力应变场数据, 作为本文算法的初始输入。
由于壁板含多个螺栓连接, 对所有螺栓进行实体建模并设置接触计算过于繁冗, 本文在保证精度的前提下通过整体-细节2级模型计算危险点应力应变: 先通过梁单元模拟螺栓, 并通过MPC连接梁单元和螺孔实际接触受力面模拟螺栓连接, 统计分析螺栓钉传载荷及旁路载荷选取危险细节; 再切取危险细节进行精细建模, 对螺栓进行实体建模并与螺孔设置接触, 从第一步提取危险细节边界节点力施加, 保证细节建模受力状态与整体一致, 进而计算孔边危险点应力应变。
利用Abaqus有限元分析软件对搭接壁板进行建模分析, 边界条件施加与实际实验保持一致。蒙皮-带板壁板的整体有限元计算结果如图 5所示。统计钉传载荷大小, 如表 3所示, 由于钉传载荷关于对称面大小相等, 因此只给出对称面一侧螺栓钉传载荷的结果。
|
| 图 5 蒙皮-带板整体有限元计算结果 |
| 排数 | 第一列 | 第二列 | 第三列 |
| 第一排 | 14 149.1 | 13 960.3 | 14 131.5 |
| 第二排 | 14 093.8 | 14 097.7 | |
| 第三排 | 12 185.6 | 12 181.9 |
由表 3可知, 对于带板来说, 第一排钉传、旁路都最大, 因此最危险; 对于蒙皮来说, 第三排钉处旁路载荷最大, 虽然钉传载荷没有第一、第二排大, 但其螺栓直径较小, 因此第三排钉最危险, 又由于此处带板的边不受约束, 因此螺栓对螺孔的挤压更严重, 产生应力集中更严重。因此对于蒙皮-带板壁板, 第三排钉处的蒙皮最危险。
提取第三排钉处的局部模型, 建立相应的螺栓, 螺母与蒙皮、带板建立接触, 螺栓与螺孔建立接触。细节边界面处网格与整体模型保持一致, 以便施加相应节点力, 且孔边网格要足够细。提取整体模型危险细节边界面上的节点力, 施加到细节模型上, 并约束任一单元以消除刚体位移。有限元计算结果如图 6所示, 可以看出带板和蒙皮危险点均在孔边, 且都未进入塑性, 而蒙皮孔边最大应力大于带板孔边最大应力, 说明蒙皮孔边最危险, 与整体有限元分析结果一致。
|
| 图 6 蒙皮-带板细节模型计算结果 |
由有限元分析结果分别提取蒙皮危险点处的应力应变, 乘以相应系数得到实验常幅谱应力应变峰谷值序列, 作为双尺度疲劳损伤演化模型非嵌入后处理程序的输入, 计算得到细观应力应变曲线, 如图 7所示。本文对模型参数S和s抽取5 000组样本, 进而利用宏细观双尺度损伤演化模型计算得到5 000个疲劳寿命样本, 在此基础上计算疲劳寿命概率分布参数, 结果如表 4所示。
|
| 图 7 搭接壁板危险点细观应力应变滞回环 |
本文完成了蒙皮-带板搭接壁板试件的常幅谱疲劳实验, 试件数量为8。加载水平与有限元分析一致, 应力比为0.06, 实验在Instron8802电液伺服疲劳实验机上进行, 搭接壁板疲劳破坏位置如图 8所示。
|
| 图 8 搭接壁板破坏位置 |
由图 8可知, 蒙皮-带板试件疲劳破坏位置为蒙皮第三排螺孔处, 与本文前述有限元分析预测的破坏位置相同。实验疲劳寿命均值和标准差如表 5所示, 对比表 4, 对于蒙皮-带板试件, 本文方法预测寿命均值相对误差为-8.70%, 标准差相对误差为-13.47%, 均处于工程可接受范围内。
| 实验寿命 | 均值 | 标准差 |
| 61 057, 134 320, 102 902, 79 167, 47 245, 55 988, 113 950, 154 064 | 93 587 | 38 997 |
假设疲劳寿命服从对数正态分布, 实验结果与本文预测的概率疲劳寿命分布如图 9所示, 实验寿命均位于预测概率寿命分布的95%置信度双侧分位数区间内, 从疲劳寿命可靠性角度验证了本文方法的有效性。
|
| 图 9 搭接壁板预测疲劳寿命概率分布 |
1) 本文建立了金属材料双尺度桥联方程, 将Lemaitre疲劳损伤演化模型概率化, 通过反演优化获取模型参数样本, 进而统计其概率特性。结合双尺度桥联方程和疲劳损伤演化概率模型, 对模型参数进行随机抽样, 通过损伤累积计算即可获得寿命样本, 统计获得寿命概率特性。
2) 本文对疲劳损伤演化模型进行了简化处理, 如忽略损伤和弹塑性本构方程的耦合、材料变形进入塑性阶段即发生损伤和材料随动强化等。对于晶体结构不同的金属材料, 其弹塑性本构方程和强化模型构建地越精确, 对疲劳损伤演化机理的表征和寿命预测也会越精确, 但相应的计算量也会大大增加, 因此本文对其进行了一定的简化, 而由此造成的偏差可通过材料实验数据反演获取的模型参数S和s得到一定的补偿, 故本文方法具有一定的合理性和普适性。
3) 完成了铝合金2024-T3光滑圆棒在6组应力水平下的疲劳实验, 基于实验数据获得了模型参数概率特性; 对某飞机搭接壁板进行了整体-细节有限元分析, 计算危险点应力应变。利用本文双尺度疲劳损伤演化模型对搭接壁板在常幅谱下的概率疲劳寿命进行预测。完成了相应的疲劳实验, 实验疲劳破坏位置与本文有限元分析预测的危险点位置一致, 实验寿命均落在模型预测寿命概率分布的95%置信度双侧分位数区间内, 验证了本文模型预测金属构件概率疲劳寿命的有效性。
| [1] | ASCE Committee on Fatigue and Fracture Reliability. Fatigue reliability[J]. ASCE Journal of the Structural Division, 1982, 108(1): 3-88. DOI:10.1061/JSDEAG.0005869 |
| [2] |
白恩军, 黄树涛, 谢里阳. 威布尔分布下小样本P-S-N曲线拟合方法[J]. 西安交通大学学报, 2019, 53(9): 96-101.
BAI Enjun, HUANG Shutao, XIE Liyang. P-S-N curves fitting of small sample with fatigue lives obeying Weibull distribution[J]. Journal of Xi'an Jiaotong University, 2019, 53(9): 96-101. (in Chinese) |
| [3] |
刘红彬, 陈伟, 刘林. 某轮盘螺栓孔处疲劳寿命预测方法与验证[J]. 西北工业大学学报, 2019, 37(4): 744-750.
LIU Hongbin, CHEN Wei, LIU Lin. Prediction method and verification of fatigue life for a turbine disk with bolt hole[J]. Journal of Northwestern Polytechnical University, 2019, 37(4): 744-750. (in Chinese) DOI:10.3969/j.issn.1000-2758.2019.04.014 |
| [4] |
罗媛, 颜东煌, 袁明. 随机车辆冲击作用下简支梁桥疲劳可靠度评估[J]. 计算力学学报, 2018, 35(4): 417-423.
LUO Yuan, YAN Donghuang, YUAN Ming. Fatigue reliability assessment of simply supported bridges subject to impact of stochastic vehicle loads[J]. Chinese Journal of Computational Mechanics, 2018, 35(4): 417-423. (in Chinese) |
| [5] |
马玉娥, 王博, 熊晓枫. 玻璃纤维铝合金层板(FMLs)的疲劳损伤特性及S-N曲线[J]. 西北工业大学学报, 2016, 34(2): 41-45.
MA Yu'e, WANG Bo, XIONG Xiaofeng. Experimental study of fatigue damage of glass-fiber reinforced aluminum laminates(FMLs)[J]. Journal of Northwestern Polytechnical University, 2016, 34(2): 41-45. (in Chinese) |
| [6] | DOUDARD C, CALLOCH S, CUGY P, et al. A probabilistic two-scale model for high-cycle fatigue life predictions[J]. Fatigue & Fracture of Engineering Materials & Structures, 2005, 28(3): 279-288. |
| [7] | SANTECCHIA E, HAMOUDA A M S, MUSHARAVATI F, et al. A review on fatigue life prediction methods for metals[J]. Advances in Materials Science and Engineering, 2016, 2016: 1-26. |
| [8] | TODINOV M T. Necessary and sufficient condition for additivity in the sense of Palmgren-Miner rule[J]. Comput Mater Sci, 2001, 21: 101-110. DOI:10.1016/S0927-0256(00)00221-4 |
| [9] | KRAJCINOVIC D, LEMAITRE J. Continuum damage mechanics: theory and applications[M]. New York: Springer-Verlag, 1987. |
| [10] | LEMAITRE J, SERMAGE J P, DESMORAT R. A two scale damage concept applied to fatigue[J]. International Journal of Fracture, 1999, 97: 67-81. DOI:10.1023/A:1018641414428 |
| [11] | ESHELBY J D. The determination of the elastic field of an ellipsoidal inclusion and related problems[J]. Proc Roy Soc London, 1957, A241: 376-196. |
| [12] | LEMAITRE J, DESMORAT R, SAUZAY M. Anisotropic damage law of evolution[J]. European Journal of Mechanics, 2000, 19: 187-208. DOI:10.1016/S0997-7538(00)00161-3 |
| [13] | STEPHENS M A. EDF statistics for goodness of fit and some comparisons[J]. Journal of the American Statistical Association, 1974, 69(347): 730-737. |
2. Research Institute of Unmanned Aerial Vehicles, Zhengzhou University of Aeronautics, Zhengzhou 450046, China;
3. School of Aeronautics, Northwestern Polytechnical University, Xi'an 710072, China
