CN106446367A - 基于弧长法非线性有限元分析的轮盘破裂转速预测方法 - Google Patents

Abstract

本发明公开了一种基于弧长法非线性有限元分析的轮盘破裂转速预测方法，步骤为：对轮盘材料取样，通过试验获取材料本构模型参数；采用传统方法近似估算轮盘的破裂转速；在有限元软件中定义轮盘的弹塑性本构模型；建立轮盘材料的三维有限元模型并进行网格划分；针对需要计算的轮盘转子有限元模型，施加约束及转速载荷；采用弧长法进行大变形非线性有限元计算；计算完成后在时间后处理器中获得转速或角速度‑轮盘径向位移曲线，从而可得到轮盘或转子的极限临界转速，即为破裂转速。本发明采用了相对更为合理基于弧长法非线性有限元分析和全局塑性失稳的轮盘破裂准则，因此能更精确地预测出轮盘或转子的极限转速。

Description

基于弧长法非线性有限元分析的轮盘破裂转速预测方法

技术领域

本发明属于材料试验方法技术领域，具体涉及一种基于弧长法非线性有限元分析和整体塑性失稳准则的轮盘破裂转速预测方法。

背景技术

在航空发动机中，轮盘是最重要的转子部件之一。对于涡轮盘，还处在较高的工作温度下，由于温度分布不均引起的热应力及高速旋转产生的离心载荷导致涡轮盘处在及其恶劣的工作环境中。为了确保发动机和飞机的正常服役，国内外的燃气涡轮发动机设计准则或规范对轮盘的强度即破裂转速性能做出了严格的要求和规定。

因此，破裂转速分析成为现代航空发动机轮盘强度设计最重要内容之一。随着发动机整体性能水平的不断提升，先进高强度合金材料在轮盘中大量应用，传统破裂转速预测方法不能满足先进材料轮盘设计分析的需要。

实际发动机轮盘几何形状复杂，存在榫槽、倒角、螺栓孔以及偏心孔等几何不连续情况，加上旋转圆盘内部本身主要存在的径向和周向双向应力，导致旋转轮盘内部处于复杂的多轴应力状态。材料在多轴应力状态下的屈服、塑性变形和失效行为都变得较为复杂。工程实践表明传统方法(如平均应力法)所采用的轮盘破裂准则与轮盘达到极限破裂转速时的实际状态并不相符，破裂转速误差较大，不能满足工程设计需要。有必要研究和建立与轮盘极限破裂转速状态吻合的破裂准则(即极限转速确定准则)，然后在此基础上发展和验证与之相适应的轮盘破裂转速计算方法。

发明内容

为弥补现有预测方法在复杂几何构型轮盘破裂转速预测精度方面的不足，本发明的目的是提供一种基于弧长法非线性有限元分析的轮盘破裂转速预测方法，以提高轮盘破裂转速判断的准确性。

为实现上述目的，本发明采用以下技术方案：

一种基于弧长法非线性有限元分析的轮盘破裂转速预测方法，包括以下步骤：

(1)对轮盘材料取样，通过试验获取材料本构模型参数；

(2)采用传统方法近似估算轮盘的破裂转速；

(3)在有限元软件中定义轮盘的弹塑性本构模型；

(4)建立轮盘材料的三维有限元模型并进行网格划分；

(5)针对需要计算的轮盘转子有限元模型，施加约束及转速载荷；

(6)采用弧长法进行大变形非线性有限元计算；

(7)计算完成后在时间后处理器中获得转速或角速度-轮盘径向位移曲线，从而可得到轮盘或转子的极限临界转速，即为破裂转速。

进一步的，所述步骤(1)中，试验获取轮盘材料的本构模型参数的步骤为：通过光滑试件拉伸试验获得材料的工程应力-应变曲线(σ_E-ε_E)，选取最高点前的曲线，由如下变换公式转化为真应力-应变曲线(σ_T-ε_T)；

上式中，σ_E表示工程应力，ε_E表示工程应变，σ_T表示真应力，ε_T表示真应变；

进一步将真应力-应变曲线有如下变换公式转换为真应力-塑性应变曲线(σ_T-ε_p)；

上式中，ε_p表示塑性应变，E表示弹性模量，该弹性模量由试验的弹性段斜率而获得；

由真应力-塑性应变曲线拟合获得材料本构硬化模型σ_T＝σ₀+r(ε_p)，式中σ₀表示初始屈服应力，r(ε_p)表示用塑性应变表达的硬化项；拟合曲线σ_T＝σ₀+r(ε_p)要求最高点，即最后一点满足即最高点斜率等于纵坐标值。

进一步的，所述步骤(2)中，由传统方法预估破裂转速ω_es，或者采用试算法大致确定破裂转速所在的范围。

进一步的，所述步骤(3)中，通过二次开发方法将材料本构硬化模型σ_T＝σ₀+r(ε_p)植入到有限元子程序中；并将步骤(1)中获得的参数采用命令流或GUI方式输入到有限元程序中。

进一步的，所述步骤(4)中，将轮盘简化为扇形模型进行建模。

进一步的，所述步骤(4)中，对已建好的三维有限元模型进行网格划分时，针对螺栓孔、榫槽、圆角这些关键部位进行局部网格细化。

进一步的，所述步骤(5)中，对步骤(4)中建立的三维有限元模型周向两端面施加周向固定约束，并对整体模型施加一个大于轮盘破裂转速估算值的转速载荷。

进一步的，所述步骤(6)中，求解器设置采用弧长法进行大变形非线性有限元计算，并设置计算的停止条件。

进一步的，所述步骤(6)中，设置计算结束时Time值为1,关闭自动控制时间步选项，根据结构的复杂程度设置子步数和并记录结果的频次。

进一步的，所述步骤(7)中，计算完成后在时间后处理器中获得转速或角速度-轮盘径向位移曲线，其曲线极值即为轮盘或转子的极限临界转速。

有益效果：本发明所述的轮盘破裂转速预测方法，采用了相对更为合理的基于弧长法非线性有限元分析和全局塑性失稳的轮盘破裂准则，有效地解决了传统数值方法无法对复杂轮盘进行准确分析以及传统有限元方法无法收敛的问题，同时使用了更贴近实际轮盘物理行为的全局塑性失稳破裂准则，因而能更精确地预测出轮盘或转子的极限转速。

附图说明

图1a是模拟轮盘试样的图纸，图1b是模拟轮盘试样的实物图片；

图2a是光滑试件拉伸试验试样图纸，图2b是光滑试件拉伸试验试样实物图片；

图3是力学性能测试取样方位图；

图4是拉伸应力-应变曲线试验结果曲线；

图5是拉伸工程应力-应变曲线和真应力-应变曲线；

图6是有限元模型及网格划分示意图；

图7是转速或角速度-轮盘径向位移曲线。

具体实施方式

下面结合附图对本发明做更进一步的解释。

本发明的一种基于弧长法非线性有限元分析的轮盘破裂转速预测方法，包括以下步骤：

(1)对轮盘材料取样，通过试验获取材料本构模型参数，具体步骤为：通过光滑试件拉伸试验获得材料的工程应力-应变曲线(σ_E-ε_E)，选取最高点前的曲线，由如下变换公式转化为真应力-应变曲线(σ_T-ε_T)；

上式中，σ_E表示工程应力，ε_E表示工程应变，σ_T表示真应力，ε_T表示真应变；

进一步将真应力-应变曲线有如下变换公式转换为真应力-塑性应变曲线(σ_T-ε_p)；

上式中，ε_p表示塑性应变，E表示弹性模量，该弹性模量由试验的弹性段斜率而获得；

由真应力-塑性应变曲线拟合获得材料本构硬化模型σ_T＝σ₀+r(ε_p)，式中σ₀表示初始屈服应力，r(ε_p)表示用塑性应变表达的硬化项；拟合曲线σ_T＝σ₀+r(ε_p)要求最高点，即最后一点满足即最高点斜率等于纵坐标值。

(2)采用传统方法近似估算轮盘的破裂转速；可由传统方法预估破裂转速ω_es，或者采用试算法大致确定破裂转速所在的范围，此步骤只需确定破裂转速的大致下限。

(3)在有限元软件中定义轮盘的弹塑性本构模型；通过二次开发方法将材料本构硬化模型σ_T＝σ₀+r(ε_p)植入到有限元子程序中；并将步骤(1)中获得的参数采用命令流或GUI方式输入到有限元程序中。

(4)建立轮盘材料的三维有限元模型并进行网格划分；对轮盘进行有限元建模，这里为节省计算时间，可将轮盘根据实际情况简化为扇形进行建模；对已建好的三维模型进行网格划分，对于螺栓孔、榫槽、圆角等关键部位应进行局部网格细化。

(5)针对需要计算的轮盘转子有限元模型，施加约束及转速载荷；对三维有限元模型周向两端面施加周向固定约束，并对整体模型施加一个大于轮盘破裂转速估算值的转速载荷，如1.5ω_es或2ω_es。

(6)采用弧长法进行大变形非线性有限元计算；求解器设置采用弧长法进行大变形非线性有限元计算，并设置计算的停止条件，例如可设置为轮盘径向伸长一定程度作为迭代计算的停止条件，设置的径向伸长量要足够大；设置计算结束时Time值为1,关闭自动控制时间步选项，根据结构的复杂程度设置子步数和并记录结果的频次。

(7)计算完成后在时间后处理器中获得转速或角速度-轮盘径向位移曲线，其曲线极值即为轮盘或转子的极限临界转速，即为破裂转速。

下面结合具体实施例对本发明做进一步说明。

实施例：

本发明以镍基高温合金GH4169模拟轮盘试样为例，在室温条件下对该模拟轮盘的破裂转速进行预测。试样图纸及实物图片如图1a和b所示。

步骤一，为准确预测GH4169模拟盘的破裂转速，需要从模拟盘锻件上取样测试其力学性能和拉伸应力-应变曲线，即材料的工程应力应变曲线，试样图纸及照片如图2a和b所示,图中尺寸为：h1’＝100mm，h2’＝28.39mm，h3’＝30mm，Φ1＝10±0.02mm，Φ2＝5±0.02mm，R＝10mm，Ra＝0.8。

GH4169模拟盘采用自由锻圆饼加工制造，粗加工后采用直接时效热处理制度进行热处理，然后按图纸精加工至最终尺寸。力学性能测试试样从同一块自由锻圆饼上取样，其中，6根径向取样，6根弦向取样，取样方位如图3所示，随粗加工后的毛坯盘同炉热处理。

模拟盘锻件取样的拉伸应力-应变曲线试验结果如图4所示。从图中可以看出，GH4169模拟盘锻件圆饼径向和弦向的力学性能差异较小，拉伸应力-应变曲线分散性较小，故在确定其本构模型参数时将选取所有拉伸应力-应变曲线的平均曲线，作为建立材料弹塑性本构模型的基本数据，即材料的工程应力-应变曲线(σ_E-ε_E)。选取最高点前的曲线，由如下变换公式转化为真应力-应变曲线(σ_T-ε_T)，

上式中，σ_E表示工程应力，ε_E表示工程应变，σ_T表示真应力，ε_T表示真应变。

上述步骤获取的GH4169室温下的拉伸工程应力-应变曲线和真应力-应变曲线如图5所示。

进一步将真应力-应变曲线有如下变换公式转换为真应力-塑性应变曲线(σ_T-ε_p)，

上式中，ε_p表示塑性应变，E表示弹性模量，该弹性模量可由试验的弹性段斜率而获得。

由真应力-塑性应变曲线拟合获得材料本构硬化模型σ_T＝σ₀+r(ε_p)，式中σ₀表示初始屈服应力，r(ε_p)表示用塑性应变表达的硬化项。拟合曲线σ_T＝σ₀+r(ε_p)要求最高点，即最后一点满足即最高点斜率等于纵坐标值，满足大变形失稳条件。

这里采用如下非线性各项同性硬化模型规律来描述其单向拉伸真应力-应变曲线的塑性变形阶段：

式中，σ₀为初始屈服应力，第二项为线性硬化项，第三项、第四项和第五项为非线性硬化项，r₁、r₂、r₃、r₄、b₁、b₂、b₃和b₄为各向同性硬化参数。

拟合结果如下表1：

表1 GH4169弹塑性本构模型材料参数拟合结果

步骤二，根据经验数据初步估计破裂转速ω_es≥25000r/min。

步骤三，将步骤一中试验获得的材料属性参数采用命令流或GUI方式输入到已植入硬化模型的Ansys子程序中。

步骤四，为简化计算过程，本例对轮盘试样模型化简为轮盘模型，即圆形角为30°的扇形模型后，进行三维有限元建模，并对其划分网格，对于螺栓孔、榫槽、圆角等关键部位应进行局部网格细化，如图6所示，本例采用Solid 185单元。

步骤五，对步骤四中建立的扇形有限元模型周向两端面施加周向固定约束，对整体模型施加一个大于轮盘破裂转速估算值的转速载荷，本例取2ω_es即50000r/min。

步骤六，求解器设置采用弧长法进行大变形非线性有限元计算，指定载荷步末的时间为1，关闭自动确定载荷步选项，设置载荷步数为2000步，每10步记录一次数据，计算时设置为轮盘外缘径向伸长15mm作为迭代计算的停止条件。

步骤七，求解，获得有限元计算结果。

步骤八，在时间后处理器中获取轮盘外缘每一载荷步的径向位移和转速，做出转速或角速度-轮盘径向位移曲线，如图7，曲线的极值36750r/min即为轮盘的破裂转速。

将本预测方法计算结果与实验结果对比如表2，从误差可以看出本发明一种能够有效预测轮盘破裂转速的方法。

表2结果对比

虽然本发明已以较佳实施例公开如上，但他们并不是用来限定本发明的，任何熟习此技艺者，在不脱离本发明之精神和范围内，自当可做各种变化或润饰，因此本发明的保护范围应当以本申请的专利保护范围所界定的为准。本发明未详尽描述的均为常规技术内容。

Claims (10)

1.一种基于弧长法非线性有限元分析的轮盘破裂转速预测方法，其特征在于：包括以下步骤：

(1)对轮盘材料取样，通过试验获取材料本构模型参数；

(2)采用传统方法近似估算轮盘的破裂转速；

(3)在有限元软件中定义轮盘的弹塑性本构模型；

(4)建立轮盘材料的三维有限元模型并进行网格划分；

(5)针对需要计算的轮盘转子有限元模型，施加约束及转速载荷；

(6)采用弧长法进行大变形非线性有限元计算；

(7)计算完成后在时间后处理器中获得转速或角速度-轮盘径向位移曲线，从而可得到轮盘或转子的极限临界转速，即为破裂转速。

2.根据权利要求1所述的基于弧长法非线性有限元分析的轮盘破裂转速预测方法，其特征在于：所述步骤(1)中，试验获取轮盘材料的本构模型参数的步骤为：通过光滑试件拉伸试验获得材料的工程应力-应变曲线(σ_E-ε_E)，选取最高点前的曲线，由如下变换公式转化为真应力-应变曲线(σ_T-ε_T)；

$\left\{\begin{array}{c}{\mathrm{\ϵ}}_T=ln({\mathrm{\ϵ}}_E+1)\\ {\mathrm{\σ}}_T={\mathrm{\σ}}_E\·({\mathrm{\ϵ}}_E+1)\end{array}\right.$

上式中，σ_E表示工程应力，ε_E表示工程应变，σ_T表示真应力，ε_T表示真应变；

进一步将真应力-应变曲线有如下变换公式转换为真应力-塑性应变曲线(σ_T-ε_p)；

${\mathrm{\ϵ}}_p={\mathrm{\ϵ}}_T-\frac{{\mathrm{\σ}}_T}{E}$

上式中，ε_p表示塑性应变，E表示弹性模量，该弹性模量由试验的弹性段斜率而获得；

由真应力-塑性应变曲线拟合获得材料本构硬化模型σ_T＝σ₀+r(ε_p)，式中σ₀表示初始屈服应力，r(ε_p)表示用塑性应变表达的硬化项；拟合曲线σ_T＝σ₀+r(ε_p)要求最高点，即最后一点满足即最高点斜率等于纵坐标值。

3.根据权利要求1所述的基于弧长法非线性有限元分析的轮盘破裂转速预测方法，其特征在于：所述步骤(2)中，由传统方法预估破裂转速ω_es，或者采用试算法大致确定破裂转速所在的范围。

4.根据权利要求1所述的基于弧长法非线性有限元分析的轮盘破裂转速预测方法，其特征在于：所述步骤(3)中，通过二次开发方法将材料本构硬化模型σ_T＝σ₀+r(ε_p)植入到有限元子程序中；并将步骤(1)中获得的参数采用命令流或GUI方式输入到有限元程序中。

5.根据权利要求1所述的基于弧长法非线性有限元分析的轮盘破裂转速预测方法，其特征在于：所述步骤(4)中，将轮盘简化为扇形模型进行建模。

6.根据权利要求1或5所述的基于弧长法非线性有限元分析的轮盘破裂转速预测方法，其特征在于：所述步骤(4)中，对已建好的三维有限元模型进行网格划分时，针对螺栓孔、榫槽、圆角这些关键部位进行局部网格细化。

7.根据权利要求1所述的基于弧长法非线性有限元分析的轮盘破裂转速预测方法，其特征在于：所述步骤(5)中，对步骤(4)中建立的三维有限元模型周向两端面施加周向固定约束，并对整体模型施加一个大于轮盘破裂转速估算值的转速载荷。

8.根据权利要求1所述的基于弧长法非线性有限元分析的轮盘破裂转速预测方法，其特征在于：所述步骤(6)中，求解器设置采用弧长法进行大变形非线性有限元计算，并设置计算的停止条件。

9.根据权利要求8所述的基于弧长法非线性有限元分析的轮盘破裂转速预测方法，其特征在于：所述步骤(6)中，设置计算结束时Time值为1,关闭自动控制时间步选项，根据结构的复杂程度设置子步数和并记录结果的频次。

10.根据权利要求1所述的基于弧长法非线性有限元分析的轮盘破裂转速预测方法，其特征在于：所述步骤(7)中，计算完成后在时间后处理器中获得转速或角速度-轮盘径向位移曲线，其曲线极值即为轮盘或转子的极限临界转速。