CN114492110A

CN114492110A - 基于权函数的轮盘表面裂纹应力强度因子计算方法及系统

Info

Publication number: CN114492110A
Application number: CN202111674199.2A
Authority: CN
Inventors: 杨晓光; 王志芳; 范永升; 杨秦政; 谭龙
Original assignee: Beihang University
Current assignee: Beihang University
Priority date: 2021-12-31
Filing date: 2021-12-31
Publication date: 2022-05-13

Abstract

本发明涉及一种基于权函数的轮盘表面裂纹应力强度因子计算方法及系统。该方法包括：根据轮盘上任一表面裂纹所在平面的横截面的几何特征确定对应简化平板中心表面裂纹的几何尺寸，并代入Newman‑Raju经验公式中求出参考解，以此来确定权函数表达式中的待定系数；然后提取裂纹所在平面的横截面上的正应力分布数据并采用最小二乘法拟合成二元多项式，最后在横截面全域上对权函数与应力分布的乘积做二维积分运算即可得到裂纹尖端的应力强度因子。本发明能够高效灵活的计算轮盘这类复杂的三维结构上裂纹的应力强度因子，不需要在裂纹的形状大小位置发生变化时重新划分网格做有限元计算，且不存在需对裂纹尖端划分高质量的网格的问题。

Description

基于权函数的轮盘表面裂纹应力强度因子计算方法及系统

技术领域

本发明涉及航空发动机结构安全性和可靠性领域，特别是涉及一种基于权函数的轮盘表面裂纹应力强度因子计算方法及系统。

背景技术

航空发动机的涡轮盘是结构关键件，在服役过程中处于高温、高转速的条件下，低周疲劳是其最主要的失效原因。在启停循环的交变载荷作用下，材料缺陷、加工制造过程中产生的表面划痕或应力水平较大的位置上容易产生疲劳裂纹，其中，表面裂纹是轮盘上常见的一种裂纹形式，裂纹进一步扩展就会导致结构发生失效。

对于含裂纹结构的安全性评估通常采用断裂力学的方法。在线弹性断裂力学理论中，应力强度因子是最基础的参量，它与含裂纹结构的几何形状与尺寸及所承受的载荷相关。虽然有应力强度因子手册可供查阅，但是应用范围仅限于一些简单的几何结构和加载形式，而实际工程结构的几何形状和应力分布都较为复杂，无法直接通过解析解或手册解计算应力强度因子。因此，应力强度因子K的准确高效计算成为复杂工程结构损伤容限设计需要解决的关键问题。

由于轮盘几何结构及应力分布复杂，限制了常规解析法和经验公式的使用，目前轮盘上裂纹应力强度因子的计算以有限元法居多，但有限元法对网格划分的要求较高且计算成本大，一旦裂纹的形状尺寸位置发生改变就需要重新划分网格提交计算，存在费时费力的缺点。

发明内容

本发明的目的是提供一种基于权函数的轮盘表面裂纹应力强度因子计算方法及系统，能够高效灵活的计算轮盘这类复杂的三维结构上裂纹的应力强度因子，不需要在裂纹的形状大小位置发生变化时重新划分网格做有限元计算，且不存在需对裂纹尖端划分高质量的网格的问题。

为实现上述目的，本发明提供了如下方案：

一种基于权函数的轮盘表面裂纹应力强度因子计算方法，包括：

将轮盘沿表面裂纹所在平面进行剖开，获取表面裂纹所在平面的横截面；根据所述横截面的几何特征，确定平板表面中心裂纹结构模型的几何尺寸；

将所述几何尺寸代入Newman-Raju经验公式中，求解出表面裂纹在两端受均布拉伸载荷作用下的应力强度因子；

将所述表面裂纹在两端受均布拉伸载荷作用下的应力强度因子作为在参考载荷作用下的应力强度因子的参考解，反解出权函数中的待定系数；根据所述待定系数，确定权函数；

提取所述横截面上的正应力分布数据，并采用最小二乘法将提取的正应力分布数据拟合成二元多项式，得到拟合后的应力分布；

在所述横截面的全域上对所述权函数与所述拟合后的应力分布的乘积做二维积分运算，得到轮盘表面裂纹的应力强度因子。

可选地，所述根据所述横截面的几何特征，确定平板表面中心裂纹结构模型的几何尺寸，具体包括：

将所述横截面的裂纹中心点作为原点，以x轴沿表面裂纹所在边的切线、y轴沿表面裂纹的深度方向建立直角坐标系；

在所述直角坐标系的第一象限和第二象限内分别作角平分线，将所述横截面分为四个区域；

在每个区域内分别寻找所述横截面边界上距离原点最近的点，将靠近x轴区域内最近的点与原点之间的距离作为半宽尺寸，将靠近y轴区域内最近的点与原点之间的距离作为厚度尺寸；

分别选取每个区域内的半宽和厚度尺寸的最小值作为所述平板表面中心裂纹结构模型的几何尺寸。

可选地，所述在参考载荷作用下的应力强度因子的参考解的计算公式为：

其中，K_r(P')在参考载荷作用下的应力强度因子的参考解；s为裂纹横截面上点(x,y)到裂纹前缘的最短距离；ρ是点(x,y)到应力强度因子计算点P'的距离；φ代表P'在裂纹前缘的位置；r与

为裂纹面上点(x,y)的极坐标，R是点(x,y)所对应的裂纹前沿上的点的极坐标；σ_r(x,y)为参考载荷作用下裂纹横截面的应力分布；M(θ,α,β)是裂纹前沿点的位置参数θ的函数；α为裂纹的形状比；β为裂纹深度与板厚之比。

可选地，所述权函数的计算公式为：

其中，m(x,y；P')为权函数；s为裂纹横截面上点(x,y)到裂纹前缘的最短距离；ρ是点(x,y)到应力强度因子计算点P'的距离；φ代表P'在裂纹前缘的位置；r与

为裂纹面上点(x,y)的极坐标，R是点(x,y)所对应的裂纹前沿上的点的极坐标；M(θ,α,β)是裂纹前沿点的位置参数θ的函数，即待定系数；α为裂纹的形状比；β为裂纹深度与板厚之比。

可选地，所述提取所述横截面上的正应力分布数据，具体包括：

在所述横截面上做多条平行于x轴和y轴的直线，提取所述直线的交点处的所有正应力分布数据。

可选地，所述采用最小二乘法将提取的正应力分布数据拟合成二元多项式，得到拟合后的应力分布，具体包括：

采用最小二乘法将提取的正应力分布数据拟合成二元多项式：

令：

ψ_k(y)＝y^k

并使得：

其中，N和M为二元多项式的次数，c_dk为多项式系数，x^d、y^k、σ_ij为提取的正应力分布数据。

可选地，所述多项式系数的矩阵表达式为：

C＝[c_dk]_(M+1)×(N+1)＝(B^TB)^-1B^TUG(G^TG)^-1

其中，

G＝[ψ_k(y_j)]_(n+1)×(N+1)

U＝[σ_ij]_(m+1)×(n+1)

可选地，所述二元多项式中的次数N和M均取3。

可选地，所述轮盘表面裂纹的应力强度因子的计算公式为：

K＝∫∫m(x,y；P')·σ(x,y)dS

其中，K为轮盘表面裂纹的应力强度因子，m(x,y；P')为权函数，σ(x,y)为拟合后的应力分布，S为横截面的全域。

一种基于权函数的轮盘表面裂纹应力强度因子计算系统，包括：

几何尺寸获取模块，用于将轮盘沿表面裂纹所在平面进行剖开，获取表面裂纹所在平面的横截面，根据所述横截面的几何特征，确定平板表面中心裂纹结构模型的几何尺寸；

两端均布载荷应力强度因子计算模块，用于将所述几何尺寸代入Newman-Raju经验公式中，求解出表面裂纹在两端受均布拉伸载荷作用下的应力强度因子；

权函数确定模块，用于将所述表面裂纹在两端受均布拉伸载荷作用下的应力强度因子作为在参考载荷作用下的应力强度因子的参考解，反解出权函数中的待定系数，根据所述待定系数，确定权函数；

应力分布获取模块，用于提取所述横截面上的正应力分布数据，并采用最小二乘法将提取的正应力分布数据拟合成二元多项式，得到拟合后的应力分布；

表面裂纹应力强度因子确定模块，用于在所述横截面的全域上对所述权函数与所述拟合后的应力分布的乘积做二维积分运算，得到轮盘表面裂纹的应力强度因子。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明根据轮盘上任一表面裂纹所在平面的横截面的几何特征确定对应简化平板中心表面裂纹的几何尺寸，并代入Newman-Raju经验公式中求出参考解，以此来确定权函数表达式中的待定系数；然后提取裂纹所在平面的横截面上的正应力分布数据并采用最小二乘法拟合成二元多项式，最后在和横截面全域上对权函数与应力分布的乘积做二维积分运算即可得到裂纹尖端的应力强度因子。与有限元法相比，本发明能够高效灵活的计算轮盘这类复杂的三维结构上裂纹的应力强度因子，不需要在裂纹的形状大小位置发生变化时重新划分网格做有限元计算，且不存在需对裂纹尖端划分高质量的网格的问题。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例基于权函数的轮盘表面裂纹应力强度因子计算方法流程图；

图2为本发明实施例基于权函数的轮盘表面裂纹应力强度因子计算系统示意图；

图3为本发明实施例某发动机涡轮盘中心孔边1/24扇形模型的应力分布情况；

图4为本发明实施例平板表面中心裂纹结构示意图；

图5为本发明实施例根据横截面形状确定对应的简化平板结构的半宽和厚度尺寸示意图；

图6为本发明实施例横截面上应力分布数据提取点示意图；

图7为本发明实施例轮盘上裂纹的应力强度因子沿裂纹前沿的分布图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1为本发明实施例提供的基于权函数的轮盘表面裂纹应力强度因子计算方法流程图。

如图1所示，本实施例中的一种基于权函数的轮盘表面裂纹应力强度因子计算方法，包括：

步骤101：将轮盘沿表面裂纹所在平面进行剖开，获取表面裂纹所在平面的横截面；根据所述横截面的几何特征，确定平板表面中心裂纹结构模型的几何尺寸。

步骤101具体为：

轮盘的实际几何特征复杂，为便于确定权函数，可将轮盘沿表面裂纹所在平面进行剖开，然后对剖开后的横截面几何特征进行合理简化，采用平板表面中心裂纹结构对含有表面裂纹的轮盘做近似，确定对应的平板表面中心裂纹结构模型的宽度和厚度尺寸。

首先，将轮盘上表面裂纹所在平面的横截面上的裂纹中心点作为原点，以x轴沿表面裂纹所在边的切线、y轴沿裂纹的深度方向建立直角坐标系。其次，在直角坐标系的两个象限(第一象限和第二象限)内分别做角平分线，这样就将横截面划分为四个区域。然后，在每个区域内寻找横截面边界上距离原点最近的点，以此定义平板表面中心裂纹结构模型的半宽和厚度尺寸(将靠近x轴区域内最近的点与原点之间的距离作为半宽尺寸，将靠近y轴区域内最近的点与原点之间的距离作为厚度尺寸)。对这四个区域重复上述操作，分别选择平板半宽和厚度尺寸的最小值作为最终尺寸，即平板表面中心裂纹结构模型的几何尺寸。

对位于轮盘表面上每一位置和方向上的裂纹，均可采用本发明实施例步骤101的方法确定其所对应的平板表面中心裂纹结构模型的平板尺寸。

步骤102：将所述几何尺寸代入Newman-Raju经验公式中，求解出表面裂纹在两端受均布拉伸载荷作用下的应力强度因子。

步骤102具体为：

将平板表面中心裂纹结构模型的半宽和厚度尺寸均代入Newman-Raju经验公式中，得到两端受均布拉伸载荷的平板表面中心裂纹的应力强度因子K，其表达式为：

其中，S表示平板两端所受均布拉伸载荷的大小，a表示裂纹深度，c表示裂纹表面半宽，

表示裂纹前沿的点所处位置；t表示平板厚度，b表示平板半宽，t和b根据步骤101的方法确定；Q为形状因子，与裂纹的形状比有关，F为边界修正因子，与裂纹的形状比、裂纹深度与平板厚度的比值、裂纹宽度与平板宽度的比值以及裂纹前沿点的位置有关。

步骤103：将所述表面裂纹在两端受均布拉伸载荷作用下的应力强度因子作为在参考载荷作用下的应力强度因子的参考解，反解出权函数中的待定系数；根据所述待定系数，确定权函数。

步骤103具体为：

首先，将表面裂纹在两端受均布拉伸载荷作用下的应力强度因子K作为在参考载荷作用下的应力强度因子的参考解K_r(P')，反解出权函数中的待定系数M(θ,α,β)，其表达式为：

其中，在参考载荷作用下，应力强度因子等于权函数与应力分布在裂纹面全域上的二维积分。K_r(P')在参考载荷作用下的应力强度因子的参考解；s为裂纹横截面上点(x,y)到裂纹前缘的最短距离；ρ是点(x,y)到应力强度因子计算点P'的距离；φ代表P'在裂纹前缘的位置；r与

为裂纹面上点(x,y)的极坐标，R是点(x,y)所对应的裂纹前沿上的点的极坐标；σ_r(x,y)为参考载荷作用下裂纹横截面的应力分布；M(θ,α,β)是裂纹前沿点的位置参数θ的函数，α为裂纹的形状比，β为裂纹深度与板厚之比。

值得注意的是，M(θ,α,β)与裂纹的形状比α＝α/c，裂纹深度与板厚之比β＝a/t有关，且裂纹前沿各每个点的M(θ,α,β)的值都不同，M(θ,α,β)是裂纹前沿点的位置参数θ的函数。

其次，根据反解出的待定系数M(θ,α,β)，采用适用于平板中心表面裂纹的二维权函数的通用形式来确定权函数，其表达式为：

步骤104：提取所述横截面上的正应力分布数据，并采用最小二乘法将提取的正应力分布数据拟合成二元多项式，得到拟合后的应力分布。

步骤104具体为：

基于权函数的表达式，可得到任意载荷分布下平板表面中心裂纹前缘任意一点处的应力强度因子的计算表达式：

此时，还需要有裂纹面上应力分布的表达式σ(x,y)。轮盘上的应力分布可通过对受服役条件下载荷作用的不含裂纹的轮盘做一次有限元计算得到。但有限元计算得到的应力分布是一系列离散的数据，为了从这些离散的数据中提炼出应力分布的表达式，本发明提出根据最小二乘拟合的原则将裂纹所在平面的横截面上的应力分布数据拟合成二元多项式。

首先，在横截面上做一系列平行于x轴和y轴的直线，提取这些直线的交点处的应力数据(x_i,y_j,σ_ij)，其中，i＝0,1,…,m；j＝0,1,…,n。

其次，采用最小二乘法将提取的正应力分布数据拟合成二元多项式：

令：

ψ_k(y)＝y^k

并使得：

兼顾考虑拟合精度与计算效率，二元多项式的次数N和M均取3，多项式系数矩阵C＝[c_dk]_(M+1)×(N+1)可根据下面的公式进行计算：

C＝(B^TB)^-1B^TUG(G^TG)^-1

其中：

G＝[ψ_k(y_j)]_(n+1)×(N+1)

U＝[σ_ij]_(m+1)×(n+1)

步骤105：在所述横截面的全域上对所述权函数与所述拟合后的应力分布的乘积做二维积分运算，得到轮盘表面裂纹的应力强度因子。

其中，基于权函数法的裂纹应力强度因子计算公式如下：

式中，x为裂纹面坐标，a为裂纹深度，σ(x)为无裂纹体在假想裂纹处应力分布，m(x,a)为权函数，它表示作用于x点处的单位集中力载荷诱导的裂纹尖端应力强度因子。

对于双向变化应力场中的二维裂纹，应力强度因子可由权函数m(x,y；P')和应力分布σ(x,y)在裂纹面全域上的二维积分计算得到：

K＝∫∫m(x,y；P')·σ(x,y)dS

其中，权函数m(x,y；P')表示x,y)点处的单位集中力载荷在裂纹前缘P'点处诱导的应力强度因子，K为轮盘表面裂纹的应力强度因子，S为横截面的全域。

因此，本发明轮盘表面裂纹的应力强度因子的计算公式最终为：

K＝∫∫m(x,y；P')·σ(x,y)dS

图2为本发明实施例提供的基于权函数的轮盘表面裂纹应力强度因子计算系统示意图。

如图2所示，本实施例中的一种基于权函数的轮盘表面裂纹应力强度因子计算系统，包括：

几何尺寸获取模块，用于将轮盘沿表面裂纹所在平面进行剖开，获取表面裂纹所在平面的横截面，根据所述横截面的几何特征，确定平板表面中心裂纹结构模型的几何尺寸。

两端均布载荷应力强度因子计算模块，用于将所述几何尺寸代入Newman-Raju经验公式中，求解出表面裂纹在两端受均布拉伸载荷作用下的应力强度因子。

权函数确定模块，用于将所述表面裂纹在两端受均布拉伸载荷作用下的应力强度因子作为在参考载荷作用下的应力强度因子的参考解，反解出权函数中的待定系数，根据所述待定系数，确定权函数。

应力分布获取模块，用于提取所述横截面上的正应力分布数据，并采用最小二乘法将提取的正应力分布数据拟合成二元多项式，得到拟合后的应力分布。

实施例1

下面以某发动机涡轮盘上中心孔处表面裂纹的应力强度因子为例，对本发明基于函数法的轮盘表面裂纹应力强度因子计算方法予以详细说明，具体实施步骤如下：

1、考虑离心载荷、热载荷以及气动载荷等复杂因素对航空发动机涡轮盘扇形做服役应力状态有限元分析，分析过程可借助ABAQUS或ANSYS等商业有限元软件进行。其中，发动机涡轮盘中心孔边1/24扇形模型的应力分布情况如图3所示，圈出的位置为孔边应力较高处，假设该处存在一个径向的表面裂纹。

2、对轮盘几何特征进行合理简化，采用如图4所示的平板表面中心裂纹结构模型做近似，以便确定权函数。模型近似简化方法如下：

2.1、将轮盘的三维模型沿确定裂纹所在平面进行切割，获取如图5所示的横截面，其中灰色半椭圆表示裂纹面。

2.2、以裂纹中心点为原点，使x轴沿裂纹所在边的切线，y轴沿裂纹的深度方向建立坐标系。

2.3、利用x和y两个坐标轴构成直角坐标系，在两个象限内分别作角平分线，如图5中虚线所示，将轮盘横截面分为四个区域。

2.4、在每个区域内寻找横截面边界上距离原点最近的点，以此定义平板的半宽或厚度尺寸，尺寸定义遵循原则为：靠近x轴的区域定义半宽尺寸，靠近y轴的区域定义厚度尺寸。

2.5、对四个区域均进行搜索之后，分别选择平板半宽和厚度的最小尺寸作为最终简化模型尺寸。在此是实施例中y轴两侧对应区域搜索到的点到原点的距离相等。

3、将轮盘简化模型的平板半宽b和平板厚度t的尺寸数据代入Newman-Raju经验公式中计算出参考解，

然后根据参考解反算出二维权函数通用表达式中的待定系数M：

其中，表面裂纹前沿各点的应力强度因子是不同的，因此权函数表达式中的待定系数M也是与裂纹前沿的点的位置相关的。

求出待定系数M之后，权函数表达式也就完全确定下来：

4、在裂纹所在平面的横截面上做一系列平行于x轴和y轴的线，如图6所示，将这线的交点即图中红色的点作为数据提取点，提取这些点的应力数据：

(x_i,y_j,σ_ij)(i＝0,1,…,m；j＝0,1,…,n)

然后采用最小二乘法将这些应力分布数据拟合成二元多项式：

ψ_k(y)＝y^k

其中，多项式的次数N和M均取3。

5、在裂纹面全域上对权函数与应力分布表达式的乘积做二维积分，从而计算得到应力强度因子：

K＝∫∫m(x,y；P')·σ(x,y)dS

最终，计算得到的该轮盘上裂纹的应力强度因子沿裂纹前沿的分布如图7所示。应力强度因子表示了裂纹扩展的驱动力与裂纹扩展速率密切相关，是线弹性断裂力学中计算裂纹扩展寿命的基本参量，同时它还被用来表征材料抵抗裂纹扩展的能力，即断裂韧度。

综上所述，权函数法是求解复杂应力场中裂纹应力强度因子高效、可靠的方法。权函数法将影响应力强度因子的几何和载荷因素做出分离，权函数仅是结构几何特征的函数，一经确定，即可求解任意载荷条件下裂纹的应力强度因子。本发明只需对不含裂纹的轮盘做一次有限元计算得到轮盘的应力分布情况，之后对于轮盘上任一表面裂纹，首先根据裂纹所在平面的横截面的几何特征确定对应简化平板中心表面裂纹的几何尺寸，并代入Newman-Raju经验公式中求出参考解，以此来确定权函数表达式中的待定系数；然后提取裂纹所在平面的横截面上的正应力分布数据并采用最小二乘法拟合成二元多项式，最后在裂纹面全域上对权函数与应力分布的乘积做二维积分运算即可得到裂纹尖端的应力强度因子。与有限元法相比，本发明所提出的方法不需要在裂纹的形状大小位置发生变化时重新划分网格做有限元计算，且不存在需对裂纹尖端划分高质量的网格的问题，本发明为轮盘这类复杂的三维结构上裂纹的应力强度因子计算提供了一种高效、灵活的方法，能够得到有效的应力强度因子是对涡轮盘等复杂结构进行裂纹扩展行为分析、寿命预测和开展损伤容限分析的前提，对于航空发动机的设计和安全适航具有重要的意义。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。