基于虚拟应变能的金属材料多轴疲劳寿命预测方法与系统

1.本发明涉及一种基于虚拟应变能的金属材料多轴疲劳寿命预测方法与系统，属于增材制造金属材料寿命预测领域，可用于航空系统技术领域。


背景技术：

2.飞机部件在服役过程中承受循环载荷。疲劳损伤是飞机结构失效的主要原因之一。航空公司为节约运营成本，对轻量化和长寿命的设计方案要求越发迫切。增材制造(am)技术适用于复杂结构近净成形,可以满足结构拓扑优化的要求，同时具有快速成型和材料利用率高的优点而受到航空工业的广泛关注。
3.研究表明，增材制造加工的零件的静力性能与常规加工方式加工出的零件相差不大，但疲劳性能差距较大，表现为疲劳寿命较低、疲劳分散性较大。因此，增材制造技术尚未应用于飞机主承力结构的制造中。由于增材制造加工工艺的特点，在加工过程中不可避免地会产生大量的缺陷，对疲劳裂纹的萌生提供了较好的环境。目前，科学界针对增材制造中的疲劳行为进行了大量的研究，可以分为微观和宏观两个研究方向。如从微观角度分析，往往通过ct扫描，显微镜拍摄等手段，通过缺陷、表面粗糙度、裂纹扩展速率等方面进行分析；宏观的分析方法分为建立针对增材制造材料的应力寿命曲线或采用孔洞附近的局部参数代替全局参数等方法。微观分析有助于深入理解疲劳破坏的原理，而宏观分析在工程应用中更为便捷。总的来说，科学界对于增材制造的研究方向，主要集中在单轴状况的预测分析，多轴载荷下的疲劳问题尚未得到很好的解决。
4.飞机的结构在服役中承受复杂的多轴载荷。例如，机翼整体受到剪切力、弯矩和扭矩的作用。机翼上的蒙皮将同时承受正应力和剪应力。飞机起落架结构在起飞和着陆过程中承受着来自多个方向的载荷。多轴疲劳失效通常是这些结构失效的原因。目前，多轴疲劳寿命的预测主要通过将多轴应力状态等效到单轴应力状态，再根据单轴应力寿命曲线预测疲劳寿命。这些方法可分为应力准则、应变准则和能量准则。能量准则认为，疲劳危险点附近能量的累积是疲劳破坏的原因，较为符合疲劳破坏的内在机理。同时能量是标量，计算过程也较为简便。然而，各种方法适用于特定的材料或载荷条件，目前还没有公认的符合工程结构设计要求的对增材制造金属材料多轴疲劳寿命预测方法。


技术实现要素：

5.发明目的：针对于上述现有技术的不足，本发明的目的在于提供一种多轴虚拟应变能处理方法，以解决增材制造金属材料多轴疲劳寿命预测问题。
6.技术方案：为达到上述目的，本发明采用的技术方案如下：
7.一种基于虚拟应变能的金属材料多轴疲劳寿命预测方法，包括如下步骤：
8.(1)对金属构件进行有限元分析，确定危险点并求出危险点处应力应变载荷时间历程；
9.(2)计算经过危险点的每个平面上的剪应变能，以剪应变能最大的平面为临界面，
得到临界面上的剪应变能以及多轴应力比；其中剪应变能定义为剪应力变幅与剪应变变幅的乘积，多轴应力比定义为临界面上正应力最大值与剪应力幅值的比值；
10.(3)计算危险点处的虚拟应变能并采用扭转载荷的剪应变能寿命曲线进行寿命预测；所述虚拟应变能为对单轴拉压载荷、扭转载荷下剪应变能寿命曲线中剪应变能进行归一化处理后的应变能。
11.作为优选，所述步骤(1)包括：
12.对于给定外载荷，采用三维建模软件绘制构件几何模型，并进行网格划分，对疲劳危险点附近进行网格细化；将材料属性赋予有限元模型，并对有限元模型添加边界条件，模拟构件在真实环境中的载荷状况；通过有限元分析得出危险点处应力应变载荷时间历程。
13.作为优选，所述步骤(2)中经过危险点的每个平面上的剪应变能计算方法包括：
14.设金属构件在时间[0,t]内承受载荷，在疲劳危险点o处建立局部笛卡尔坐标系oxyz，经过危险点o的任意平面δ可以用在其平面上的两个方向向量和一个平面的法向量描述。具体地，向量a为平面δ与坐标系oxyz中平面oxy的交线，向量b为在平面δ上与a轴垂直的向量，n为平面δ的法向量。向量n,a,b与坐标系oxyz的位置关系可以用三个角度表示，其中，为平面δ的方向法向量n与x轴的夹角；θ为平面δ的方向法向量n与z轴的夹角,α为平面δ上任意向量q与向量a的夹角。危险点o处的应力应变载荷时间历程用如下矩阵表示：
[0015][0016][0017]
式中，σ
x
(t),σy(t),σz(t)为正应力分量，τ
xy
(t),τ
yz
(t),τ
xz
(t)为剪应力分量，ε
x
(t),εy(t),εz(t)为正应变分量，γ
xy
(t),γ
yz
(t),γ
xz
(t)剪应变分量，t为[0,t]区间内的任一时刻；
[0018]
描述平面δ的三个单位向量n,a,b通过坐标变换表示为：
[0019][0020][0021][0022]
平面δ上任意单位向量q写作：
[0023][0024]
任意角度的正应力σn(t)、剪应力τq(t)、剪应变γq(t)写作：
[0025]
σn(t)＝n
t
[σ(t)]n
[0026]
τq(t)＝q
t
[σ(t)]n
[0027]
γq(t)＝2q
t
[ε(t)]n
[0028]
任意平面上的正应力均值σ
n,m
定义为[0,t]区间内正应力积分的平均值，正应力幅值σ
n,a
通过正应力的方差定义
[0029][0030][0031]
式中，
[0032][0033]
同理，可以定义剪应力、剪应变的均值和幅值τm,τa,γm,γa:
[0034][0035][0036][0037][0038]
式中，
[0039][0040][0041]
剪应变能δτδγ定义为构件在稳定循环下包围迟滞回线的最小矩形，即为剪应力变幅δτ与剪应变变幅δγ的乘积，变幅δτ＝2τa，δγ＝2γa，于是δτδγ＝4τa×
γa多轴应力比ρ定义为正应力最大值σ
n,max
与剪应力幅值τa的比值：
[0042][0043]
作为优选，所述步骤(3)中单轴拉压载荷、扭转载荷下剪应变能寿命曲线根据如下方法确定：
[0044]
分别计算单轴拉压、扭转载荷下每个疲劳数据点临界面上的剪应变能，首尾相连
绘制剪应变能-寿命曲线，以幂函数分别拟合求解两种载荷状况下，剪应变能-寿命曲线的各项参数a1,b1,a3,b3；两种载荷下的剪应变能-寿命曲线函数表示为：
[0045]
单轴拉压载荷
[0046]
扭转载荷
[0047]
其中，δτδγ为剪应变能，nf为寿命。
[0048]
作为优选，所述步骤(3)中归一化后的虚拟应变能wv表示为：
[0049][0050]
其中，为归一化系数，ρ为多轴应力比。
[0051]
基于相同的发明构思，本发明提供的一种基于虚拟应变能的金属材料多轴疲劳寿命预测系统，包括输入模块，处理模块以及输出模块，所述输入模块用于输入金属构件几何模型以及外载荷；所述输出模块用于显示预测的金属构件在给定外载荷下的寿命；所述处理模块包括：
[0052]
应力应变计算单元，用于对金属构件进行有限元分析，确定危险点并求出危险点处应力应变载荷时间历程；
[0053]
临界面计算单元，用于计算经过危险点的每个平面的剪应变能，以剪应变能最大的平面为临界面，得到临界面上的剪应变能以及多轴应力比；其中剪应变能定义为剪应力变幅与剪应变变幅的乘积，多轴应力比定义为临界面上正应力最大值与剪应力幅值的比值；
[0054]
以及预测单元，用于计算危险点处的虚拟应变能并采用扭转载荷的剪应变能寿命曲线进行寿命预测；所述虚拟应变能为对单轴拉压载荷、扭转载荷下剪应变能寿命曲线中剪应变能进行归一化处理后的应变能。
[0055]
作为优选，所述单轴拉压载荷、扭转载荷下剪应变能寿命曲线由幂函数拟合单元生成，所述幂函数拟合单元，用于分别计算单轴拉压、扭转载荷下每个疲劳数据点临界面上的剪应变能，首尾相连绘制剪应变能-寿命曲线，以幂函数分别拟合求解两种载荷状况下，剪应变能-寿命曲线的各项参数a1,b1,a3,b3。
[0056]
作为优选，所述预测单元中归一化后的虚拟应变能wv表示为：
[0057][0058]
基于相同的发明构思，本发明提供的一种计算机系统，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述计算机程序被加载至处理器时实现所述的基于虚拟应变能的金属材料多轴疲劳寿命预测方法。
[0059]
有益效果：本发明分析了典型韧性金属材料单轴拉压、扭转状态下的剪应变能寿命曲线，发现两种载荷状况下的剪应变能与寿命具有较好的双对数线性关系，基于如下两点考虑将剪应变能进行归一化处理得到虚拟应变能，采用扭转载荷的剪应变能寿命曲线进行寿命预测，从而能够用于不同载荷路径的寿命分析问题：一是韧性金属材料的疲劳裂纹行为体现为剪破坏，采用扭转载荷的剪应变能寿命曲线更为直观；二是扭转载荷的剪应变能寿命曲线往往位于所有曲线的最上方，拟合到该曲线下更为直观，且归一化系数形式较为简单。此外，传统应力准则往往用于塑性应变较小的高周疲劳，应变准则往往用于塑性应
变较大的低周疲劳，本发明采用应力与应变相乘的能量准则可以兼顾低周、高周疲劳。
附图说明
[0060]
图1为本发明实施例的总体流程图。
[0061]
图2为本发明实施例中所涉及的危险点处局部坐标系定义示意图。
[0062]
图3为本发明实施例中疲劳寿命详细计算流程图。
[0063]
图4为本发明实验例的316l钢单轴拉压、扭转载荷下剪应变能随寿命变化曲线图。
[0064]
图5为本发明实验例的有限元建模网格划分示意图。
[0065]
图6为本发明实验例的预测结果与试验结果比较图。
[0066]
图7为本发明实施例的模块结构示意图。
具体实施方式
[0067]
为了便于本领域技术人员的理解，下面结合实施例与附图对本发明作进一步的说明，实施方式提及的内容并非对本发明的限定。
[0068]
如图1所示，本发明实施例公开的一种基于虚拟应变能的金属材料多轴疲劳寿命预测方法，主要步骤如下：
[0069]
(1)对金属构件进行有限元分析，确定危险点并求出危险点处应力应变载荷时间历程。
[0070]
主要包括：
[0071]
对于给定外载荷，采用三维建模软件建立构件的几何模型，导入有限元软件；
[0072]
对几何模型进行网格划分，并对疲劳危险点附近进行网格细化，直到计算结果收敛；
[0073]
将材料属性赋予有限元模型，根据材料在疲劳试验过程有无塑性应变，选择线弹性或弹塑性本构模型；
[0074]
对有限元模型添加边界条件，模拟构件在真实环境的载荷状况；
[0075]
最后有限元分析得出危险点o处应力应变载荷时间历程以备后续计算。
[0076]
(2)计算经过危险点的每个平面上的剪应变能，以剪应变能最大的平面为临界面，得到临界面上的剪应变能以及多轴应力比。主要包括：
[0077]
设金属构件在时间[0,t]内承受载荷，在疲劳危险点o处建立局部笛卡尔坐标系oxyz，经过危险点o的任意平面δ可以用在其平面上的两个方向向量和一个平面的法向量描述。具体地，向量a为平面δ与坐标系oxyz中平面oxy的交线，向量b为在平面δ上与a轴垂直的向量，n为平面δ的法向量。向量n,a,b与坐标系oxyz的位置关系可以用三个角度表示，其中，为平面δ的方向法向量n与x轴的夹角；θ为平面δ的方向法向量n与z轴的夹角,α为平面δ上任意向量q与向量a的夹角。危险点o处的应力应变载荷时间历程可以用如下矩阵表示：
[0078]
[0079][0080]
式中，σ
x
(t),σy(t),σz(t)为正应力分量，τ
xy
(t),τ
yz
(t),τ
xz
(t)为剪应力分量，ε
x
(t),εy(t),εz(t)为正应变分量，γ
xy
(t),γ
yz
(t),γ
xz
(t)剪应变分量，t为[0,t]区间内的任一时刻。
[0081]
描述平面δ的三个单位向量n,a,b通过坐标变换表示为：
[0082][0083][0084][0085]
平面δ上任意单位向量q可以写作：
[0086][0087]
任意角度的正应力σn(t)、剪应力τq(t)、剪应变γq(t)可以写作：
[0088]
σn(t)＝n
t
[σ(t)]n
[0089]
τq(t)＝q
t
[σ(t)]n
[0090]
γq(t)＝2q
t
[ε(t)]n
[0091]
任意平面上的正应力均值σ
n,m
定义为[0,t]区间内正应力积分的平均值，正应力幅值σ
n,a
通过正应力的方差定义
[0092][0093][0094]
式中，
[0095][0096]
同理，可以定义剪应力、剪应变的均值和幅值τm,τa,γm,γa:
[0097]
[0098][0099][0100][0101]
式中，
[0102][0103][0104]
剪应变能δτδγ定义为构件在稳定循环下包围迟滞回线的最小矩形，即为剪应力变幅δτ与剪应变变幅δγ的乘积，变幅δτ＝2τa，δγ＝2γa，于是
[0105]
δτδγ＝4τa×
γa[0106]
多轴应力比ρ定义为正应力最大值σ
n,max
与剪应力幅值τa的比值，即：
[0107][0108]
(3)计算危险点处的虚拟应变能并采用扭转载荷的剪应变能寿命曲线进行寿命预测。这里虚拟应变能为对单轴拉压载荷、扭转载荷下剪应变能寿命曲线中剪应变能进行归一化处理后的应变能。单轴拉压载荷、扭转载荷下剪应变能寿命曲线预先根据材料单轴拉压、扭转疲劳试验数据确定。具体为：
[0109]
分别计算单轴拉压、扭转载荷下每个疲劳数据点临界面上的剪应变能，首尾相连绘制剪应变能-寿命曲线，拟合求解两种载荷状况下，剪应变能-寿命曲线的各项参数a1,b1,a3,b3。
[0110]
分析典型韧性金属材料单轴拉压、扭转状态下的剪应变能-寿命曲线。可以发现两种载荷状况下的剪应变能与寿命具有较好的双对数线性关系，分别对两条曲线进行幂函数拟合，即可用以描述对应寿命的剪应变能。
[0111]
单轴拉压载荷
[0112]
扭转载荷
[0113]
虚拟应变能通过归一化不同载荷路径的剪应变能，从而计算疲劳寿命：
[0114][0115][0116]
式中，多轴应力比ρ来反映不同载荷路径的影响。
[0117]
基于虚拟应变能的扭转载荷下寿命曲线计算疲劳寿命的详细过程可见图3。
[0118]
下面以一具体实验例对本发明方案做进一步验证。本例采用的材料是由选区激光熔化生产的316l不锈钢，进行了多种载荷路径的疲劳试验，即单轴拉压、扭转、比例和90
°
非
比例多轴载荷。首先对试验件进行有限元分析，确定危险点并求出危险点处应力应变载荷时间历程[σ(t)]，[ε(t)]；然后通过坐标变换求解临界面上剪应变能δτδγ和多轴应力比ρ；再根据材料单轴拉压、扭转疲劳试验数据确定模型的参数a1,b1,a3,b3；最后计算待求载荷下危险点处虚拟应变能wv并求解寿命nf。详细过程如下：
[0119]
s1：本例采用的试验件材料为增材制造316l钢，标距段长度为38mm，直径为12mm的实心圆棒。采用mts809拉扭试验机进行不含平均应力的单轴拉压、扭转、多轴比例载荷、90
°
非比例载荷和包含平均应力的单轴拉压、多轴比例载荷、90
°
非比例载荷试验。单轴拉压、扭转试验件疲劳数据用以确定本发明模型的参数，其他载荷路径的疲劳数据用以验证本发明模型的正确性，具体的试验数据如表1所示，其中载荷比r为应力最小值与最大值之比，用以描述平均应力的大小，相位角φ为正应力与剪应力的相位差，用以描述非比例程度。
[0120]
表1 316l钢疲劳试验结果
[0121][0122][0123]
本例采用有限元软件patran&nastran进行有限元分析。以试验件质心为原点，试验件轴线方向为x轴，建立整体坐标系。由于光滑试件标距段表面每个点应力应变状态都相同，裂纹可能从标距段表面任意位置萌生，不妨假定疲劳危险点为坐标(0,6,0)的点。
[0124]
网格划分采取先建立二维壳单元，再沿轴线旋转扫掠生成三维实体单元的方法，
并对疲劳危险点局部加密，1mm划分20个单元，有限元建模结果如图5所示。
[0125]
在疲劳试验过程中，未发现塑性应变或塑性应变可以忽略，因此采用线弹性本构模型，弹性模量e＝190.8gpa,泊松比v＝0.3由静力试验测得。
[0126]
为模拟试验件在试验机上的受载状况，对试验件一端固支，另一端建立mpc并施加载荷。
[0127]
最后进行有限元计算，并导出应力应变载荷时间历程[σ(t)]，[ε(t)]。
[0128]
s2：采用matlab软件进行编程，对三个角度从0
°
到180
°
，步长为1，分别求解每个平面δ的剪应变能，找出其中剪应变能最大的平面定位临界面，并求解本模型需要的参数，列于表2的3-7列。
[0129]
表2 316l钢采用本模型计算的数据
[0130][0131]
s3：绘制单轴拉压(编号p-03至p-12)，扭转(编号p-21至p-28)状况下剪应变能δτδγ-寿命nf曲线，如图4所示，采用matlab拟合程序或excel拟合功能对曲线进行幂函数拟合，求出a1＝5.42,b1＝-0.13,a3＝3.03,b3＝-0.06。
[0132]
s4：根据模型
[0133]
[0134][0135]
求解各个载荷状况下的虚拟应变能值wv，并计算寿命n
f,e
，分别列于表2的第8、9列。为了便于比较结果，计算结果与试验结果绘制与图6。由图可见，79％的寿命预测点处于2倍误差带内，93％的寿命预测点处于3倍误差带内，对各种载荷路径的疲劳寿命预测均取得了较好的预测效果。
[0136]
基于相同的发明构思，本发明实施例提供的一种基于虚拟应变能的金属材料多轴疲劳寿命预测系统，如图7所示，包括输入模块，处理模块以及输出模块，所述输入模块用于输入金属构件几何模型以及外载荷；所述输出模块用于显示预测的金属构件在给定外载荷下的寿命；所述处理模块包括：应力应变计算单元，用于对金属构件进行有限元分析，确定危险点并求出危险点处应力应变载荷时间历程；临界面计算单元，用于计算经过危险点的每个平面的剪应变能，以剪应变能最大的平面为临界面，得到临界面上的剪应变能以及多轴应力比；其中剪应变能定义为剪应力变幅与剪应变变幅的乘积，多轴应力比定义为临界面上正应力最大值与剪应力幅值的比值；以及预测单元，用于计算危险点处的虚拟应变能并采用扭转载荷的剪应变能寿命曲线进行寿命预测；所述虚拟应变能为对单轴拉压载荷、扭转载荷下剪应变能寿命曲线中剪应变能进行归一化处理后的应变能。单轴拉压载荷、扭转载荷下剪应变能寿命曲线由幂函数拟合单元生成，所述幂函数拟合单元，用于分别计算单轴拉压、扭转载荷下每个疲劳数据点临界面上的剪应变能，首尾相连绘制剪应变能-寿命曲线，以幂函数分别拟合求解两种载荷状况下，剪应变能-寿命曲线的各项参数。
[0137]
本领域技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的各单元的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。所述单元的划分仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元可以结合或者可以集成到另一个系统。
[0138]
基于相同的发明构思，本发明实施例提供的一种计算机系统，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述计算机程序被加载至处理器时实现所述的基于虚拟应变能的金属材料多轴疲劳寿命预测方法。
[0139]
本领域技术人员可以理解的是，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明实施例所述方法的全部或部分步骤。存储介质包括：u盘、移动硬盘、只读存储器rom、随机存取存储器ram、磁碟或者光盘等各种可以存储计算机程序的介质。
[0140]
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。