一种基于应变能理论的横断型斜裂纹转子变刚度特性计算方法.pdf

本发明公开了一种基于应变能理论的横断型斜裂纹转子变刚度特性计算方法，其特征在于，该方法包括如下步骤：1)无裂纹转子轴单元弹性应变能计算；2)无裂纹转子轴单元的位移计算；3)三种类型裂纹转子轴应力密度因子计算；4)裂纹引起的转子轴单元外加弹性应变能计算；5)裂纹引起的转子轴单元外加位移计算；6)裂纹转子轴单元柔度系数计算；7)基于静平衡变换的变刚度特性求解。本发明的有益效果是，方法设计合理，计算精确。

1.  一种基于应变能理论的横断型斜裂纹转子变刚度特性计算方法，其特征在于，该方法包括如下步骤：
1)无裂纹转子轴单元弹性应变能计算；
2)无裂纹转子轴单元的位移计算；
3)三种类型裂纹转子轴应力密度因子计算；
4)裂纹引起的转子轴单元外加弹性应变能计算；
5)裂纹引起的转子轴单元外加位移计算；
6)裂纹转子轴单元柔度系数计算；
7)基于静平衡变换的变刚度特性求解。

2.  根据权利要求1所述的基于应变能理论的横断型斜裂纹转子变刚度特性计算方法，其特征在于，所述无裂纹转子轴单元弹性应变能计算式为：U0|E=12∫[M12EI+M22EI+T2GI0+F2AE]dx,]]>式中M₁、M₂为弯矩，T为扭矩，F为轴向力。G为刚性模量，E为杨氏弹性模量。/为裂纹截面惯性矩，/₀为裂纹截面极惯矩。

3.  根据权利要求1所述的基于应变能理论的横断型斜裂纹转子变刚度特性计算方法，其特征在于，所述无裂纹转子轴单元的位移计算式为：u10|E=∂U0∂P1=P1LAE,u20|E=∂U0∂P2=P2L33EI-P6L22EI,u30|E=∂U0∂P3=P3L33EI+P5L22EI,u40|E=∂U0∂P4=P4LGI0,u50|E=∂U0∂P5=P5LEI+P3L22EI,u60|E=∂U0∂P6=P6LEI-P2L22EI..]]>

4.  根据权利要求1所述的基于应变能理论的横断型斜裂纹转子变刚度特性计算方法，其特征在于，所述三种类型裂纹转子轴应力密度因子计算式为：
张开模式S/F：

滑移模式S/F：

剪切模式SIF：
K3III=kP3pR2sinqpaFIII(a/h),K4III=2P4pR4b sinqpaFIII(a/h),K1III=K2III=K5III=K6III=0.]]>
其中
F1=2hπαtan(πα2h)0.752+2.02(α/h)+0.37[1-sin(πα/2h)]3cos(πα/2h)F2=2hπαtan(πα2h)0.923+0.199[1-sin(πα/2h)]4cos(πα/2h)FII=1.122-0.561(α/h)+0.085(α/h)2+0.18(α/h)31-(α/h)FIII=2hπαtan(πα2h).]]>

5.  根据权利要求1所述的基于应变能理论的横断型斜裂纹转子变刚度特性计算方法，其特征在于，所述裂纹引起的转子轴单元外加弹性应变能计算式为：J(A)|E=1E′[(Σi=16KIi|E)2+(Σi=16KIIi|E)2+η(Σi=16KIIIi|E)2]·,]]>式中E′＝E/(1-v)且η＝1+v。v为泊松比，K_Ii|_E为裂纹位移张开模式的应力密度因子，K_IIi|_E为裂纹位移滑移模式的应力密度因子，K_IIIi|_E为裂纹位移剪切模式的应力密度因子，i＝1，...，6。

6.  根据权利要求1所述的基于应变能理论的横断型斜裂纹转子变刚度特性计算方法，其特征在于，裂纹引起的转子轴单元外加位移计 算式为：uic|E=∂Uc|E∂Pi,Uc|E=∫AJ(A)|EdA.,]]>式中J(A)|_E为采用欧拉梁单元建模时根据断裂力学概念给出的应变能密度函数，由权利要去5所述的式求得。

7.  根据权利要求1所述的基于应变能理论的横断型斜裂纹转子变刚度特性计算方法，其特征在于，所述裂纹转子轴单元柔度系数g_ij|_E计算式为：


8.  根据权利要求1所述的基于应变能理论的横断型斜裂纹转子变刚度特性计算方法，其特征在于，所述基于静平衡变换的变刚度特性求解式为：[K]^c|_E＝[T]G|_E[T]^T，式中G|_E为柔度矩阵，由[G]|_E＝{g_ij|_E}，i，j＝1，...，6给出，柔度系数g_ij|_E可参照权利要求7所述内容计算得出，T为逆变矩阵，
[T]T=100000-1000000100000-1000l00100000-10-l0000100000-1000000100000-1000000100000-1,]]>[K]^c|_E为基于静平衡变换的变刚度矩阵。

一种基于应变能理论的横断型斜裂纹转子变刚度特性计算方法
技术领域
本发明基于应变能理论，提出一种横断型斜裂纹转子变刚度特性计算方法。转子轴分别采用欧拉梁与铁木辛科梁单元建模，并考虑了纵向、弯曲及扭转六个方向的自由度。根据柔度系数计算裂纹单元的刚度矩阵，而柔度系数依据裂纹单元的应变能理论求得。该计算方法为分析与揭示横断型斜裂纹转子系统振动的非线性特性奠定了理论基础。
背景技术
裂纹是一种重要的转子系统故障，如果不及时加以监测将导致灾难性后果。现代的大型旋转机械设备中，转子系统通常工作于严酷的热应力与机械应力工况下，裂纹往往频繁发生，因此有关裂纹建模与故障诊断问题得到广泛重视和研究。已有的裂纹模型可分为两大类，即常开裂纹^[1]与呼吸裂纹^[2]。常开裂纹将造成转子轴刚度的一个局部非对称常值削弱，适用于拉应力始终作用于裂纹面的情形。在呼吸裂纹情况下，拉应力与压应力随时间变化交替作用于裂纹面上^[3]，可以观察到裂纹呈现出不断张开、闭合的呼吸形态，从而引起转子轴刚度的周期变化。目前，有关裂纹建模与故障诊断研究大都针对最简单的横向表面裂纹^[4-8]，这限制了研究成果的实际应用。由于问题的复杂性，已知的斜裂纹研究比较少。Darpe建立了包含斜裂纹的转子有限元模型，并基于断裂力学推导出一个新的柔度系数矩阵，该矩阵中包含由裂纹方向引入的外加应力密度因子。进而，对比了斜裂纹与横向裂纹转子不同的刚度参数以及耦合振动响应特性^[9]。Sekhar等对比了斜裂纹与横向裂纹转子的模型特性，并对不同深度、不同位置的裂纹进行了诊断研究^[10]。
裂纹源有多种，比如穿晶裂纹、机械疲劳裂纹、汇流裂纹、晶界 裂纹、热疲劳裂纹、淬火裂纹、龟裂等，产生原因十分复杂。根据失效分析理论，不同的裂纹源以及不同的材料和形状其裂纹生成与扩展方向均不相同。Darpe与Sekhar研究的是垂直型斜裂纹，即裂纹面垂直于xoy平面，如图2所示。
实际上，另一种类型的斜裂纹对于失效分析也是十分重要的，即横断型斜裂纹，裂纹表面横断线垂直于转子轴轴线但与xoy面不垂直，如图3所示。
发明内容
本发明的目的是为了解决上述问题，设计了一种基于应变能理论的横断型斜裂纹转子变刚度特性计算方法。
实现上述目的本发明的技术方案为，一种基于应变能理论的横断型斜裂纹转子变刚度特性计算方法，其特征在于，该方法包括如下步骤：
1)无裂纹转子轴单元弹性应变能计算；
2)无裂纹转子轴单元的位移计算；
3)三种类型裂纹转子轴应力密度因子计算；
4)裂纹引起的转子轴单元外加弹性应变能计算；
5)裂纹引起的转子轴单元外加位移计算；
6)裂纹转子轴单元柔度系数计算；
7)基于静平衡变换的变刚度特性求解。
所述无裂纹转子轴单元弹性应变能计算式为：
式中M₁、M₂为弯矩，T为扭矩，F为轴向力。G为刚性模量，E为杨氏弹性模量。I为裂纹截面惯性矩，I₀为裂纹截面极惯矩。
所述无裂纹转子轴单元的位移计算式为：
u10|E=∂U0∂P1=P1LAE,u20|E=∂U0∂P2=P2L33EI-P6L22EI,u30|E=∂U0∂P3=P3L33EI+P5L22EI,u40|E=∂U0∂P4=P4LGI0,u50|E=∂U0∂P5=P5LEI+P3L22EI,u60|E=∂U0∂P6=P6LEI-P2L22EI..]]>
所述三种类型裂纹转子轴应力密度因子计算式为：
张开模式S/F：
K1I=P1pR2sin2qpaF1(a/h),K2IkP2pR2sin2qpaF1(a/h),K4I=2P4pR4R2-b2sin2qpaF1(a/h),K5I=4(P5+P3(x-(R2-b2-0.5a)ctgq))pR4bsin2qpaF2(a/h),K6I=4(P2(x-(R2-b2-0.5a)ctgq)-P6)pR4R2-b2sin2qpaF2(a/h),K3I=0.]]>
滑移模式S/F：
K1II=P1pR2sinqcosqpaFII(a/h),K2II=kP2pR2cos2qpaFII(a/h),K4II=2P4R2-b2pR4cos2qpaFII(a/h),]]>
K5II=4(P5+P3(x-(R2-b2-0.5a)ctgq))pR4bsinqcosqpaFII(a/h),K6II=4(P2(x-(R2-b2-0.5a)ctgq)-P6)pR4R2-b2sinqcosqpaFII(a/h),K3II=0]]>
剪切模式SIF：
K3III=kp3pR2sinqpaFIII(a/h),K4III=2P4pR4bsinqpaFIII(a/h),K1III=K2III=K5III=K6III=0.]]>
其中
F1=2hπαtan(πα2h)0.752+2.02(α/h)+0.37[1-sin(πα/2h)]3cos(πα/2h)F2=2hπαtan(πα2h)0.923+0.199[1-sin(πα/2h)]4cos(πα/2h)FII=1.122-0.561(α/h)+0.085(α/h)2+0.18(α/h)31-(α/h)FIII=2hπαtan(πα2h).]]>
所述裂纹引起的转子轴单元外加弹性应变能计算式为：
J(A)|E=1E'[(Σi=16KIi|E)2+(Σi=16KIIi|E)2+η(Σi=16KIIIi|E)2]·,]]>式中E′＝E/(1-v)且η＝1+v。v为泊松比，K_Ii|_E为裂纹位移张开模式的应力密度 因子，K_IIi|_E为裂纹位移滑移模式的应力密度因子，K_IIIi|_E为裂纹位移剪切模式的应力密度因子，i＝1，...，6。
裂纹引起的转子轴单元外加位移计算式为：
式中J(A)|_E为采用欧拉梁单元建模时根据断裂力学概念给出的应变能密度函数，由权利要去5所述的式求得。
所述裂纹转子轴单元柔度系数g_ij|_E计算式为：

所述基于静平衡变换的变刚度特性求解式为：[K]^c|_E＝[T]G|_E[T]^T，式中G|_E为柔度矩阵，由[G]|_E＝{g_ij|_E}，i，j＝1，…，6给出，柔度系数g_ij|_E可参照权利要求7所述内容计算得出，T为逆变矩阵，
[T]T=100000-1000000100000-1000l00100000-10-l0000100000-1000000100000-1000000100000-1,]]>[K]^c|_E为基于静平衡变换的变刚度矩阵。
利用本发明的技术方案制作的基于应变能理论的横断型斜裂纹转子变刚度特性计算方法，设计合理，计算结果准确度高，可很好的计算出基于应变能理论的横断型斜裂纹转子变刚度特性。
附图说明
图1是本发明所述基于应变能理论的横断型斜裂纹转子变刚度特性计算方法的流程示意图：
图2垂直型斜裂纹转子受力坐标系统及裂纹面；
图3横向型斜裂纹转子受力坐标系统及裂纹面。
具体实施方式
下面结合附图对本发明进行具体描述，如图1是本发明所述基于应变能理论的横断型斜裂纹转子变刚度特性计算方法的流程示意图，如图所示，转子轴分别采用欧拉梁单元与铁木辛科梁单元建模两种，并考虑了纵向、弯曲及扭转六个方向的自由度如图3为一段包含横断型斜裂纹的转子轴，转盘质量为m，转轴直径为D，长度为L。裂纹位于左端x处，裂纹深度为a。单元受力情况如下：剪力P₂、P₃与 P₈、P₉，弯矩P₅、P₆与P₁₁、P₁₂，轴向力P₁与P₇，扭矩P₄与P₁₀。
本技术方案以欧拉梁单元建模转子轴为例子阐述求解基于静平衡变换的变刚度特性的过程，其基本求解原理为：根据柔度系数计算裂纹单元的刚度矩阵，而柔度系数依据裂纹单元的Castingliano应变能理论求得，即
ui=∂U∂Pi=∂U0∂Pi+∂Uc∂Pi=ui0+uic]]>
式中U⁰为无裂纹转子轴的应变能，U^c为裂纹导致的附加应变能。u_i与P_i分别为沿着第i个坐标方向的位移和力。
实施例1
基于欧拉梁单元建模的横断型斜裂纹转子变刚度特性计算
采用欧拉梁单元建模时，忽略了裂纹面上剪切变形，只考虑轴向力、扭矩及弯矩作用。这种情况下无裂纹转子轴单元的弹性应变能U₀|_E可表达为
U0|E=12∫[M12EI+M22EI+T2GI0+F2AE]dx]]>
式中M₁、M₂为弯矩，T为扭矩，F为轴向力。G为刚性模量，E为杨氏弹性模量。I为裂纹截面惯性矩，I₀为裂纹截面极惯矩。
由图2，有M₁＝P₂x-P₆，M₂＝P₃x+P₅，T＝P₄以及F＝P₁。从而式(2)变为
U0|E=12[P12LAE+P22L33EI+P32L33EI+P42LGI0+P52LEI+P62LEI-P2P6L2EI+P3P5L2EI]]]>
无裂纹欧拉梁单元位移为
u10|E=∂U0∂P1=P1LAE,u20|E=∂U0∂P2=P2L33EI-P6L22EI,u30|E=∂U0∂P3=P3L33EI+P5L22EI,u40|E=∂U0∂P4=P4LGI0,u50|E=∂U0∂P5=P5LEI+P3L22EI,u60|E=∂U0∂P6=P6LEI-P2L22EI.]]>
裂纹引起的外加位移为
U^c|_E＝∫_AJ(A)|_EdA.
式中J(A)|_E为采用欧拉梁单元建模时根据断裂力学概念给出的应变能密度函数，由下式给出
J(A)|E=1E′[(Σi=16KIi|E)2+(Σi=16KIIi|E)2+η(Σi=16KIIIi|E)2].]]>
式中E′＝E/(1-v)且η＝1+v。v为泊松比，J_Ii|_E为裂纹位移张开模式的应力密度因子，K_IIi|_E为裂纹位移滑移模式的应力密度因子，K_IIIi|_E为裂纹位移剪切模式的应力密度因子，i＝1，...，6。
这些应力密度因子(SIF)由下式给出：
(1)张开模式SIF：
K1I=P1pR2sin2qpaF1(a/h),K2I=kP2pR2sin2qapF1(a/h),K4I=2P4pR4R2-b2sin2qpaF1(a/h),K5I=4(P5+P3(x-(R2-b2-0.5a)ctgq))pR4bsin2qpaF2(a/h),K6I=4(P2(x-(R2-b2-0.5a)ctgq)-P6)pR4R2-b2sin2qpaF2(a/h),K3I=0.]]>
(2)滑移模式SIF：
K1II=P1pR2sinq cos qpaFII(a/h),K2II=kP2pR2cos2qpaFII(a/h),]]>
K4II2P4R2-b2pR4cos2qpaFII(a/h),K5II=4(P5+P3(x-(R2-b2-0.5a)ctgq))pR4bsinqcosqpaFII(a/h),K6II=4(P2(x-(R2-b2-0.5a)ctgq)-P6)pR4R2-b2sinq cos qpaFII(a/h),K3II=0]]>
(3)剪切模式SIF：
K3III=kP3pR2sin qpaFIII(a/h),K4III=2P4pR4b sin qpaFIII(a/h),K1III=K2III=K5III=K6III=0.]]>
其中
F1=2hπαtan(πα2h)0.752+2.02(α/h)+0.37[1-sin(πα/2h)]3cos(πα/2h)F2=2hπαtan(πα2h)0.923+0.199[1-sin(πα/2h)]4cos(πα/2h)FII=1.122-0.561(α/h)+0.085(α/h)2+0.18(α/h)31-(α/h)FIII=2hπαtan(πα2h)]]>
应用这些应力密度因子(SIF)公式以及式(6)的应变能J(A)|_E表达式，求得裂纹引起的外加位移u_i^c|_E，i＝1，…，6。进而，整体位移u_i|_E可以通过加上无裂纹转子单元位移u_i⁰|_E求得，可以写成矩阵形式u_i|_E＝u_i⁰|_E+u_i^c|_E.＝[G]|_EP_i，其中G|_E为柔度矩阵，由[G]|_E＝{g_ij|_E}，i，j＝1，…，6给出。柔度系数g_ij|_E为

考虑有限元静平衡条件，采用下式由柔度矩阵通过变换T求解刚度矩阵：
{q_1-12}^T＝[T]{q_1-6}^T，
式中变换矩阵为
[T]T=100000-1000000100000-1000l00100000-10-l0000100000-1000000100000-1000000100000-1]]>
从而横断式斜裂纹转子欧拉梁单元刚度矩阵为
[K]^c|_E＝[T]G|_E[T]^T
实施例2
基于铁木辛科梁单元建模的横断型斜裂纹转子变刚度特性计算
对于铁木辛科梁单元，裂纹面上的剪切变形必须加以考虑，从而无裂纹转子轴单元的弹性应变能U⁰|_T可写为
U0|T=12∫[αsV12GA+αsV22GA+M12EI+M22EI+T2GI0+F2AE]dx]]>
式中V₁、V₂为剪切力，且V₁＝P₂，V₂＝P₃。α_S为剪切变形因子。
无裂纹欧拉梁单元位移u_i⁰|_T为
u10|T=u10|E,u40|T=u40|E,u50|T=u50|E,u60|T=u60|E,u20|T=∂U0|T∂P2=(αsLGA+L33EI)P2-P6L22EI,u30|T=∂U0|T∂P3=(αsLGA+L33EI)P3+P5L22EI.]]>
当分别采用欧拉梁与铁木辛科梁单元进行建模时，应力密度因子(SIF)有所不同，不同之处由下式给出：
K_II2|_E＝K_III3|_E＝0.           欧拉梁单元
KII2|T=kP2PR2paFII(a/h),KIII3|T=kP3pR2paFIII(a/h).]]>       铁木辛科梁单元
式中k为剪力分布修正因子。其他的应力密度因子是一样的，可参考前面的公式。
最终，横断型斜裂纹下铁木辛科梁单元与欧拉梁单元在柔度系数上的差异由下式给出：
g22|E=L33EI+∫A[2k2αFII2πER4+8x2h2αF22πER8]dA,g33|E=L33EI+∫A[32x2αβ2F12πER8+2ηk2αFIII2πER4]dA,g24|E=g42|E=0,g34|E=g43|E=0.]]>  欧拉梁单元
g22|T=αsLGA+L33EI+∫A[2k2αFII2πER4+8x2h2αF22πER8]dA,g33|T=αsLGA+L33EI+∫A[32x2αβ2F12πER8+2ηk2αFIII2πER4]dA,g24|T=g42|T∫A4kαβFII2πER6dA,g34|T=g43|T=∫A2ηkhαFIII2πER6dA.]]>  铁木辛科梁单元
其他的柔度系数是一样的，可参考前面的公式。计算出柔度系数g_ij|_T后，即可参照式实施例1中的公式计算横断型斜裂纹转子铁木辛科梁单元的刚度矩阵[K]^c|_T。
上述技术方案仅体现了本发明技术方案的优选技术方案，本技术领域的技术人员对其中某些部分所可能做出的一些变动均体现了本发明的原理，属于本发明的保护范围之内。