一种基于有限元结果的多轴疲劳寿命矩阵化计算方法与流程

1.本发明涉及疲劳寿命计算技术领域，具体涉及一种基于有限元结果的多轴疲劳寿命矩阵化计算方法。


背景技术：

2.在工程实践中，疲劳寿命的评估通常涉及多轴应力状态下的复杂几何体，这使得疲劳寿命的计算变得非常复杂。大量实验表明，裂纹总是在某些特定的材料平面上形成和生长，因此基于临界平面的疲劳准则得到了广泛的应用，这些准则的算法都包含寻找具有最大疲劳损伤的材料临界平面方向。但是，临界平面的计算成本很高，现有算法需要基于每个点的应力应变历史对整个取向空间进行遍历计算。
3.对于真实的工程案例，结构中的载荷历史通常通过有限元软件计算得到。对于大型结构，有限元模型通常在空间上有几十万或上百万个节点，在时间上存在几百或上千个离散点。对于这种情况，现有临界平面疲劳评估方法可能变得过于昂贵，特别是在早期设计迭代阶段，单次的疲劳寿命需要消耗几十甚至上百小时。


技术实现要素：

4.本发明的目的在于提供一种基于有限元结果的多轴疲劳寿命矩阵化计算方法，以解决现有技术中临界平面疲劳评估方法可能变得过于昂贵，特别是在早期设计迭代阶段，单次的疲劳寿命需要消耗几十甚至上百小时的技术问题。
5.为解决上述技术问题，本发明具体提供下述技术方案：
6.一种基于有限元结果的多轴疲劳寿命矩阵化计算方法，包括以下步骤：
7.步骤100、基于全局三维坐标轴对应的应力/应变张量，利用坐标变换和矩阵乘法计算斜面的三维坐标轴对应的应力/应变张量，将应力张量中的正应力/应变数据和剪应力/应变数据进行分类；
8.步骤200、设定斜面的应力/应变张量的整个应力/应变数据的存储结构，将随时间变化的应力/应变存储为三维数组，建立应力/应变存储数组；
9.步骤300、将斜面三维坐标轴的空间离散矩阵存储为三维数组，建立空间离散坐标变换数组；
10.步骤400、将空间离散坐标变换数组的三维数据与应力/应变存储数组的三维数据进行矩阵乘法运算，得到空间离散后的所有斜面应力/应变的空间变换应力数组；
11.步骤500、根据不同的临界面准则得到正应变临界面的计算公式以及剪应变临界面计算公式，且分别提取所述空间变换应力/应变数组中的正应力/应变数据和剪应力/应变数据，计算所有节点对应的最大正应变幅值和最大剪应变幅值；
12.步骤600、基于最大正应变幅值和最大剪应变幅值，以及swt疲劳寿命计算公式以及fatemi-socie疲劳寿命计算公式，求解所有节点对应的疲劳寿命。
13.作为本发明的一种优选方案，在步骤100中，斜面的三维坐标轴对应的应力张量具
有对称性，采用voigt方法获取斜面应力/应变的计算公式为：
14.{σ
″
}＝[t]{σ}
[0015]
{ε
″
}＝[t]{ε}
[0016]
其中，σ为任意应力状态的全局坐标xyz下的应力张量，ε为任意应力状态的全局坐标xyz下的应变张量；
[0017]
[t]为6
×
6的坐标变换矩阵；
[0018]
σ
″
为斜面三维坐标下的应力张量，且{σ
″
}和{σ}为采用voigt记法6
×
1的列向量。
[0019]
ε
″
为斜面三维坐标下的应变张量，且{ε
″
}和{ε}为采用voigt记法6
×
1的列向量。
[0020]
作为本发明的一种优选方案，利用voigt方法获取斜面应力计算方式，根据矩阵乘法计算所述斜面上的正应力/应变和剪应力/应变，将{σ
″
}中的正应力和剪应力进行分类，将{ε
″
}中的正应变和剪应变进行分类。
[0021]
作为本发明的一种优选方案，在步骤200中，所述应力应变存储数组的第一维度代表节点数目，长度为i；
[0022]
第二维度代表时间离散数目，长度为j；
[0023]
第三维度代表应力/应变向量，长度为6；
[0024]
所述应力/应变存储数组的表达式为(i，j，6)。
[0025]
作为本发明的一种优选方案，在步骤300中，所述空间离散坐标变换数组为将所述坐标变换矩阵进行三维空间离散，所述空间离散坐标变换数组的第一维度表示任意应力/应变状态的三维空间离散的坐标系数目，长度为k；
[0026]
第二维度和第三维度表示所述坐标变换矩阵的维度，长度均为6；
[0027]
所述空间离散坐标变换数组的表达式为(k，6，6)。
[0028]
作为本发明的一种优选方案，在步骤400中，空间离散后的所有斜面应力/应变的空间变换应力/应变数组为四维数组，所述空间变换应力/应变数组的大小为(i，j，k，6)。
[0029]
作为本发明的一种优选方案，在步骤500中，首先，最大正应变幅值作为临界面准则，计算所述空间变换应变数组，从所述空间变换应变数组中提取正应变数据，计算斜面最大正应力幅值作为临界面并记录其所对应的坐标系，具体的实现步骤为：
[0030]
基于ε
″
中的正应变分类，取出空间变换应变数组(i，j，k，6)中的正应力部分，建立正应变子集(i，j，k，1)；
[0031]
对该子集的第2维j做最大值减去最小值的运算，得到斜面上的正应变幅值，在第三维k方向求最大值，即可得到所有节点对应的最大正应变幅值和临界面对应的斜面坐标系。
[0032]
作为本发明的一种优选方案，在步骤600中，首先，以最大剪应变幅值作为临界面准则，结合所述空间变换应变数组计算最大剪应变幅值临界面，从所述空间变换应变数组中提出剪应变数据，来计算斜面最大剪应变临界面及对应的坐标系，具体的实现步骤为：
[0033]
基于ε
″
中的剪应变分类，取出空间变换应变数组(i，j，k，6)中的剪应变部分，建立剪应变子集(i，j，k，2)；
[0034]
对该子集的第2维j和第4维做最小内接矩形运算，得到斜面上的剪应变幅值，在第三维k方向求最大值，即可得到所有节点对应的最大剪应变幅值和临界面对应的斜面坐标系。
[0035]
作为本发明的一种优选方案，将计算的最大剪应变幅值以及最大正应变幅值，分别代入计算疲劳寿命nf的两个公式内，通过牛顿法对非线性方程进行求解，即可得到所有节点所对应的疲劳寿命nf；其中，计算疲劳寿命nf分别为：
[0036][0037]
式中σ
n,max
代表临界面上的最大正应力，代表正应变幅值，σ
′f为疲劳强度系数，e为杨氏模量，nf为疲劳寿命，b为疲劳强度指数，ε
′f为疲劳延展系数，c为疲劳延展指数。
[0038][0039]
式中σ
n,max
代表临界面上的最大正应力，k为材料参数，σy为材料的屈服应力，代表剪应变幅值，τ
′f为剪切疲劳强度系数，g为剪切模量，nf为疲劳寿命，b
γ
为剪切疲劳强度指数，γ
′f为剪切疲劳延展系数，c
γ
为剪切疲劳延展指数。
[0040]
本发明与现有技术相比较具有如下有益效果：
[0041]
本发明适用于各种临界面疲劳模型的求解和计算，基于最一般的三维复杂载荷情况，适用于任何几何形状，适用于任意复杂的应力应变输入，且矩阵乘法非常适合并行计算，有丰富和成熟的代码库，并行效率极高。只要通过提高计算机硬件的配置，就可以实现极短时间内复杂结构的寿命评估。
附图说明
[0042]
为了更清楚地说明本发明的实施方式或现有技术中的技术方案，下面将对实施方式或现有技术描述中所需要使用的附图作简单地介绍。显而易见地，下面描述中的附图仅仅是示例性的，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图引伸获得其它的实施附图。
[0043]
图1为本发明实施例提供的临界面任意应力状态的三维坐标图；
[0044]
图2为本发明实施例提供的斜面的应力图；
[0045]
图3为本发明实施例提供的斜面应力的计算示意图；
[0046]
图4为本发明实施例提供的正应力的矩阵计算示意图；
[0047]
图5为本发明实施例提供的剪切力的矩阵计算示意图；
[0048]
图6为本发明实施例提供的应力/应变存储数组的三维图；
[0049]
图7为本发明实施例提供的空间离散坐标变换数组的三维图；
[0050]
图8为本发明实施例提供的空间变换应力数组的三维图；
[0051]
图9为本发明实施例提供的多轴疲劳寿命矩阵化计算流程图。
具体实施方式
[0052]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他
实施例，都属于本发明保护的范围。
[0053]
如图9所示，本发明提供了一种基于有限元结果的多轴疲劳寿命矩阵化计算方法，传统临界面搜索方法需要基于每个点的应力应变历史对整个取向空间进行遍历计算，如图1和图2所示，计算斜面最大正应力的临界面的计算流程如下：
[0054]
已知任意应力状态：
[0055][0056]
则斜面δ上的正应力σn(t)和剪应力τn(t)的计算公式为
[0057]nx
＝sin(θ)
·
cos(φ)
[0058]
ny＝sin(θ)
·
sin(φ)
[0059]
nz＝cos(θ)
[0060]ax
＝-sin(φ)
[0061]ay
＝cos(φ)
[0062]az
＝0
[0063]bx
＝-cos(θ)
·
cos(φ)
[0064]by
＝-cos(θ)
·
sin(φ)
[0065]bz
＝sin(θ)
[0066]
t
x
(t)＝σ
x
(t)
·nx
τ
xy
(t)
·
ny τ
xz
(t)
·
nz[0067]
ty(t)＝τ
xy
(t)
·nx
σy(t)
·
ny τ
yz
(t)
·
nz[0068]
tz(t)＝τ
xz
(t)
·nx
τ
yz
(t)
·
ny σz(t)
·
nz[0069]
σn(t)＝t
x
(t)
·nx
ty(t)
·
ny tz(t)
·
nz[0070]
τ
na
(t)＝t
x
(t)
·ax
ty(t)
·ay
tz(t)
·az
[0071]
τ
nb
(t)＝t
x
(t)
·bx
ty(t)
·by
tz(t)
·bz
[0072][0073]
此时如果我们想检索最大正应力的临界面，则算法如下：
[0074]
①
循环遍历有限元模型中的节点；
[0075]
②
循环遍历theta角；
[0076]
③
循环遍历phi角；
[0077]
④
循环遍历载荷时间t；
[0078]
⑤
计算当前角度和时间下的σn(t)；
[0079]
⑥
计算当前theta，phi面下的δσn；
[0080]
⑦
所有theta，phi面下的最大的δσn；
[0081]
则该theta，phi面为当前节点的临界面，假设有限元中节点数为n，theta角数量为i，phi角数量为j，载荷时间数量为m。则该算法共包含n*i*j*m次循环，该算法效率低，不利于并行化计算。
[0082]
而本实施方式基于有限元计算结果，设计了临界面疲劳模型的矩阵化计算方法，便于实现大型结构的并行化计算，大大提高了结构寿命计算的效率。
[0083]
本实施方式包括以下步骤，如图9所示：
[0084]
步骤100、基于全局三维坐标轴对应的应力/应变张量，利用坐标变换和矩阵乘法计算斜面的三维坐标轴对应的应力/应变张量，将应力张量中的正应力/应变数据和剪应力/应变数据进行分类。
[0085]
结合附图1和图2所示，σ为全局坐标xyz下的应力张量，σ
″
为斜面坐标nba下的应力张量，水平面为xyz坐标系的xy平面法向为z。阴影面为nba坐标系的ab平面，法向为n。全局坐标系是人为定义的一个参考坐标系。则传统的斜面坐标nba下的应力张量的计算公式为：σ
″
＝mσm
t
；
[0086]
其中，
[0087][0088]
即
[0089]
则有，斜面上的正应力σn(t)＝σ
″
11
；
[0090]
斜面上的剪应力τ
na
(t)＝σ
″
；τ
nb
(t)＝σ
″
13
；
[0091]
在步骤100中，斜面的三维坐标轴对应的应力张量具有对称性，采用voigt方法获取斜面应力/应变的计算公式为：
[0092]
{σ
″
}＝[t]{σ}
[0093]
{ε
″
}＝[t]{ε}
[0094]
其中，σ为任意应力状态的全局坐标xyz下的应力张量，ε为任意应力状态的全局坐标xyz下的应变张量；
[0095]
[t]为6
×
6的坐标变换矩阵；
[0096]
σ
″
为斜面三维坐标下的应力张量，且{σ
″
}和{σ}为采用voigt记法6
×
1的列向量。
[0097]
ε
″
为斜面三维坐标下的应变张量，且{ε
″
}和{ε}为采用voigt记法6
×
1的列向量。
[0098]
其中，
[0099]
[t]矩阵与m矩阵的关系如下，“**”代表乘方：
[0100]
t11＝m11**2；t12＝m12**2；t13＝m13**2；t14＝2*m11*m12；t15＝2*m11*m13；t16＝2*m12*m13；
[0101]
t21＝m21**2；t22＝m22**2；t21＝m23**2；t24＝2*m21*m22；t25＝2*m21*m23；t26＝2*m22*m23；
[0102]
t31＝m31**2；t32＝m32**2；t33＝m33**2；t34＝2*m31*m32；t35＝2*m31*m33；t36＝2*m32*m33；
[0103]
t41＝m11*m21；t42＝m12*m22；t43＝m13*m23；t44＝m12*m21 m11*m22；t45＝m13*m21 m11*m23；t46＝m13*m22 m12*m23；
[0104]
t51＝m11*m31；t52＝m12*m32；t53＝m13*m33；t54＝m12*m31 m11*m32；t55＝m13*m31 m11*m33；t56＝m13*m32 m12*m33；
[0105]
t61＝m21*m31；t62＝m22*m32；t63＝m23*m33；t64＝m22*m31 m21*m32；t65＝m23*m31 m21*m33；t66＝m23*m32 m22*m33。
[0106]
利用voigt方法获取斜面应力计算方式，根据矩阵乘法计算所述斜面上的正应力/应变和剪应力/应变，将{σ
″’
}中的正应力和剪应力进行分类，将{ε
″
}中的正应变和剪应变进行分类。如图4和图5所示，基于正应力和剪应力的定义。1方向为平面法向因此为正应力，2,3方向为平面的切向因此为剪应力，同样的斜面上的正应力σn(t)＝σ
″
；斜面上的剪应力τ
na
(t)＝σ
″
12
；τ
nb
(t)＝σ
″
13
。
[0107]
步骤200、设定斜面的应力/应变张量的整个应力/应变数据的存储结构，将随时间变化的应力/应变存储为三维数组，建立应力/应变存储数组。
[0108]
如图6所示，所述应力/应变存储数组的第一维度代表节点数目，长度为i；
[0109]
第二维度代表时间离散数目，长度为j；
[0110]
第三维度代表应力/应变向量，长度为6；
[0111]
所述应力/应变存储数组的表达式为(i，j，6)。
[0112]
步骤300、将斜面三维坐标轴的空间离散矩阵存储为三维数组，建立空间离散坐标变换数组；
[0113]
如图7所示，在步骤300中，所述空间离散坐标变换数组为将所述坐标变换矩阵进行三维空间离散，所述空间离散坐标变换数组的第一维度表示任意应力/应变状态的三维空间离散的坐标系数目，长度为k；
[0114]
第二维度和第三维度表示所述坐标变换矩阵的维度，长度均为6；
[0115]
所述空间离散坐标变换数组的表达式为(k，6，6)，，如果对3维空间按180
°×
180
°
进行离散，则该矩阵的大小为(32400，6，6)。
[0116]
步骤400、将空间离散坐标变换数组的三维数据与应力/应变存储数组的三维数据进行矩阵乘法运算，得到空间离散后的所有斜面应力/应变的空间变换应力数组。
[0117]
在步骤400中，空间离散后的所有斜面应力/应变的空间变换应力/应变数组为四维数组，所述空间变换应力/应变数组的大小为(i，j，k，6)
[0118]
如图8所示，空间离散后的所有斜面应力的空间变换应力数组为四维数组，所述空间变换应力数组的大小为(i，j，k，6)，本实施方式通过一次矩阵乘法就实现了上述传统算法中
①②③④
步的所有运算，即采用矩阵运算代替了传统算法中所有的循环、判断和标量代数运算。
[0119]
步骤500、根据不同的临界面准则得到正应变临界面的计算公式以及剪应变临界面计算公式，且分别提取所述空间变换应力/应变数组中的正应力/应变数据和剪应力/应变数据，计算所有节点对应的最大正应变幅值和最大剪应变幅值。
[0120]
步骤600、基于最大正应变幅值和最大剪应变幅值，以及swt疲劳寿命计算公式以及fatemi-socie疲劳寿命计算公式，求解所有节点对应的疲劳寿命。
[0121]
在步骤500中，首先，最大正应变幅值作为临界面准则，计算所述空间变换应变数
组，从所述空间变换应变数组中提取正应变数据，计算斜面最大正应力幅值作为临界面并记录其所对应的坐标系，具体的实现步骤为：
[0122]
基于ε
″
中的正应变分类，取出空间变换应变数组(i，j，k，6)中的正应力部分，建立正应变子集(i，j，k，1)；
[0123]
对该子集的第2维j做最大值减去最小值的运算，得到斜面上的正应变幅值，在第三维k方向求最大值，即可得到所有节点对应的最大正应变幅值和临界面对应的斜面坐标系。
[0124]
在步骤600中，首先，以最大剪应变幅值作为临界面准则，结合所述空间变换应变数组计算最大剪应变幅值临界面，从所述空间变换应变数组中提出剪应变数据，来计算斜面最大剪应变临界面及对应的坐标系，具体的实现步骤为：
[0125]
基于ε
″
中的剪应变分类，取出空间变换应变数组(i，j，k，6)中的剪应变部分，建立剪应变子集(i，j，k，2)；
[0126]
对该子集的第2维j和第4维做最小内接矩形运算，得到斜面上的剪应变幅值，在第三维k方向求最大值，即可得到所有节点对应的最大剪应变幅值和临界面对应的斜面坐标系。
[0127]
将计算的最大剪应变幅值以及最大正应变幅值，分别代入计算疲劳寿命nf的两个公式内，通过牛顿法对非线性方程进行求解，即可得到所有节点所对应的疲劳寿命nf；其中，计算疲劳寿命nf分别为：，
[0128]
其中，计算疲劳寿命nf分别为：
[0129][0130]
式中σ
n,max
代表临界面上的最大正应力，代表正应变幅值，σ
′f为疲劳强度系数，e为杨氏模量，nf为疲劳寿命，b为疲劳强度指数，ε
′f为疲劳延展系数，c为疲劳延展指数。
[0131][0132]
式中σ
n,max
代表临界面上的最大正应力，k为材料参数，σy为材料的屈服应力，代表剪应变幅值，τ
′f为剪切疲劳强度系数，g为剪切模量，nf为疲劳寿命，b
γ
为剪切疲劳强度指数，γ
′f为剪切疲劳延展系数，c
γ
为剪切疲劳延展指数。
[0133]
将最大正应变(i,1)和最大剪正应变(i,1)分别代入上述两个公式(1)(2)，疲劳寿命nf为带求变量，其余为材料参数可以通过试验或手册获取。通过牛顿法对非线性方程进行求解，即可得到所有节点所对应的疲劳寿命nf。
[0134]
因此本实施方式适用于各种临界面模型的求解和计算，基于最一般的三维复杂载荷情况，适用于任何几何形状，适用于任意复杂的应力应变输入，且矩阵乘法非常适合并行计算，有丰富和成熟的代码库，并行效率极高。只要通过提高计算机硬件的配置，就可以实现极短时间内复杂结构的寿命评估。
[0135]
以上实施例仅为本技术的示例性实施例，不用于限制本技术，本技术的保护范围由权利要求书限定。本领域技术人员可以在本技术的实质和保护范围内，对本技术做出各
种修改或等同替换，这种修改或等同替换也应视为落在本技术的保护范围内。