本发明涉及确定转子变刚度的技术领域，特别是涉及一种基于空间任意斜裂纹的转子变刚度确定方法及系统。
背景技术：
现代转子系统长时间工作于严酷的机械、热应力环境中，疲劳裂纹时有发生，如果得不到及时监测将严重威胁系统的运行安全。裂纹动力学是转子裂纹监测的基础，涉及裂纹的动力学建模与分析等内容。
现有的研究针对横向裂纹(裂纹面同时垂直于转子轴线和基面)的情况较多，如图3所示。近年来，研究方向开始关注于斜裂纹。相比于横向裂纹，斜裂纹会导致刚度矩阵中不同参数的更多耦合，进而造成转子弯曲、扭转甚至是纵向振动的耦合。bachschmid等指出，通常情况下转子裂纹面垂直于转轴轴线方向(即横向)，但在大扭矩和强弯矩载荷作用下裂纹会沿着螺旋方向扩展，从而形成螺旋裂纹或斜裂纹。它们以精细的三维非线性模型为基础，建立起更为简单的裂纹转子一维有限元模型，据此对横裂纹与斜裂纹的呼吸机理、弯扭联合作用下裂纹的扩展行为进行数值仿真与实验验证。han等认为齿轮啮合力导致的大扭矩往往引发裂纹的斜向扩展，因此对含有斜裂纹的齿轮-转子-轴承系统的动力学特性进行了数值仿真研究。国内刘长利等对横裂纹与斜裂纹转子的变刚度特性进行了对比，并进一步研究了双盘、双呼吸型裂纹转子系统的典型非线性动力学特性。夏恒恒等分析了斜裂纹转子在裂纹全部张开状态下刚度的影响因素。但是，这些确定转子变刚度特性的绝大部分针对的是裂纹面倾斜于转轴轴线垂直于转轴基面的横-斜或直-斜裂纹，如图4所示。而对于如图5-6所示的空间任意斜向裂纹(即裂纹面与转轴轴线和转轴基面均不垂直)，其几何模型相比于横向、横-斜向的几何模型有着本质的区别，因此它们确定转子变刚度特性方法也有很大差别，另外，由于其推导过程极端复杂，所以迄今为止并没有针对空间任意斜裂纹的转子变刚度特性进行研究。基于上述问题，如何确定基于空间任意斜裂纹的转子变刚度特性成为本领域亟需解决的技术问题。
技术实现要素：
本发明的目的是提供一种基于空间任意斜裂纹转子的变刚度确定方法及系统，以实现对空间任意斜裂纹转子的变刚度的确定，便于后期根据变刚度设计转子强度以及对转子进行故障诊断。
为实现上述目的，本发明提供一种基于空间任意斜裂纹转子的变刚度确定方法，所述方法包括：
确定空间任意斜裂纹转子的总应变能模型；
根据所述总应变能模型确定空间任意斜裂纹转子的柔度参数；
根据所述柔度参数确定空间任意斜裂纹转子的变刚度模型；
将待检测的空间任意斜裂纹转子的参数值输入所述变刚度模型确定空间任意斜裂纹转子的变刚度；所述待检测的空间任意斜裂纹转子的参数值为空间任意斜裂纹转子的长度l、直径d、裂纹深度a以及裂纹面空间方位角度θ1和θ2中至少一者。
可选的，所述确定空间任意斜裂纹转子的总应变能模型，具体包括：
获取无裂纹转子的应变能；
确定裂纹面张开、滑移与剪开三种位移模式的应力密度因子；
根据裂纹面张开、滑移与剪开三种位移模式的应力密度因子确定外加应变能；
根据无裂纹转子的应变能和外加应变能确定空间任意斜裂纹转子的总应变能模型。
可选的，所述根据裂纹面张开、滑移与剪开三种位移模式的应力密度因子确定外加应变能，具体公式为：
其中，e'＝e/(1-v)，ms＝1+v，v为泊松比，kiⅰ、kiⅱ与kiⅲ分别为描述裂纹面张开、滑移与剪开三种位移模式的应力密度因子，uc为外加应变能，e为杨氏弹性模量，a为转子轴横截面积，da为积分运算中对a求导数。
可选的，所述根据所述总应变能模型确定空间任意斜裂纹转子的柔度参数，具体公式为：
其中，u0+uc为总应变能模型，u0为无裂纹转子的应变能，uc为外加应变能，ui为沿着第i个坐标方向的节点位移与节点力，pi为沿着第i个坐标方向的节点力，pj为沿着第j个坐标方向的节点力。
可选的，所述根据所述柔度参数确定空间任意斜裂纹转子的变刚度模型，具体公式为：
kc＝tg-1tt；
其中，g＝{gij},i＝1、2、…、6，j＝1、2、…、6，gij为空间任意斜裂纹转子的柔度参数，t为变换矩阵。
本发明还提供一种基于空间任意斜裂纹转子的变刚度确定系统，所述系统包括：
总应变能模型确定模块，用于确定空间任意斜裂纹转子的总应变能模型；
柔度参数确定模块，用于根据所述总应变能模型确定空间任意斜裂纹转子的柔度参数；
变刚度模型确定模块，用于根据所述柔度参数确定空间任意斜裂纹转子的变刚度模型；
变刚度确定模块，用于将待检测的空间任意斜裂纹转子的参数值输入所述变刚度模型确定空间任意斜裂纹转子的变刚度；所述待检测的空间任意斜裂纹转子的参数值为空间任意斜裂纹转子的长度l、直径d、裂纹深度a以及裂纹面空间方位角度θ1和θ2中至少一者。
可选的，所述总应变能模型确定模块，具体包括：
获取单元，用于获取无裂纹转子的应变能；
应力密度因子确定单元，用于确定裂纹面张开、滑移与剪开三种位移模式的应力密度因子；
外加应变能确定单元，用于根据裂纹面张开、滑移与剪开三种位移模式的应力密度因子确定外加应变能；
总应变能模型确定单元，用于根据无裂纹转子的应变能和外加应变能确定空间任意斜裂纹转子的总应变能模型。
可选的，所述根据裂纹面张开、滑移与剪开三种位移模式的应力密度因子确定外加应变能，具体公式为：
其中，e'＝e/(1-v)，ms＝1+v，v为泊松比，kiⅰ、kiⅱ与kiⅲ分别为描述裂纹面张开、滑移与剪开三种位移模式的应力密度因子，uc为外加应变能，e为杨氏弹性模量，a为转子轴横截面积，da为积分运算中对a求导数。
可选的，所述根据所述总应变能模型确定空间任意斜裂纹转子的柔度参数，具体公式为：
其中，u0+uc为总应变能模型，u0为无裂纹转子的应变能，uc为外加应变能，ui为沿着第i个坐标方向的节点位移与节点力，pi为沿着第i个坐标方向的节点力，pj为沿着第j个坐标方向的节点力。
可选的，所述根据所述柔度参数确定空间任意斜裂纹转子的变刚度模型，具体公式为：
kc＝tg-1tt；
其中，g＝{gij},i＝1、2、…、6，j＝1、2、…、6，gij为空间任意斜裂纹转子的柔度参数，t为变换矩阵。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明先确定空间任意斜裂纹转子的总应变能模型；然后根据所述总应变能模型确定空间任意斜裂纹转子的柔度参数；根据所述柔度参数确定空间任意斜裂纹转子的变刚度模型；最后将待检测的空间任意斜裂纹转子的参数值输入所述变刚度模型确定空间任意斜裂纹转子的变刚度；以实现对空间任意斜裂纹转子的变刚度的确定，便于后期根据变刚度设计转子强度以及对转子进行故障诊断。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例基于空间任意斜裂纹转子的变刚度确定方法流程图；
图2为本发明实施例基于空间任意斜裂纹转子的变刚度确定系统结构图；
图3为本发明实施例横向裂纹转子受力图；
图4为本发明实施例横-斜裂纹转子受力图；
图5为本发明实施例空间任意斜裂纹转子受力图；
图6为本发明实施例三类裂纹转子受力综合对比图；
图7为本发明实施例三类裂纹转子的变刚度对比图；
图8为本发明实施例空间任意斜裂纹转子的变刚度特性随着裂纹面空间方位角的变化图；
图9为本发明实施例空间任意斜裂纹转子的变刚度特性随着裂纹深度的变化图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明的目的是提供一种基于空间任意斜裂纹转子的变刚度确定方法及系统，以实现对空间任意斜裂纹转子的变刚度的确定，便于后期根据变刚度设计转子强度以及对转子进行故障诊断。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
图1为本发明实施例基于空间任意斜裂纹转子的变刚度确定方法流程图；如图1所示，本发明提供一种基于空间任意斜裂纹转子的变刚度确定方法，所述方法包括：
步骤10：确定空间任意斜裂纹转子的总应变能模型.
步骤20：根据所述总应变能模型确定空间任意斜裂纹转子的柔度参数。
步骤30：根据所述柔度参数确定空间任意斜裂纹转子的变刚度模型。
步骤40：将待检测的空间任意斜裂纹转子的参数值输入所述变刚度模型确定空间任意斜裂纹转子的变刚度。所述待检测的空间任意斜裂纹转子的参数值为空间任意斜裂纹转子的长度l、直径d、裂纹深度a以及裂纹面空间方位角度θ1和θ2中至少一者。
下面对各个步骤进行详细分析：
步骤10：确定空间任意斜裂纹转子的总应变能模型，具体包括：
步骤101：获取无裂纹转子的应变能u0。
步骤102：确定裂纹面张开、滑移与剪开三种位移模式的应力密度因子；
(1)张开模式的应力密度因子公式为：
(2)滑移模式的应力密度因子公式为：:
(3)剪开模式的应力密度因子公式为：
其中，f1、f2分别是张开(i型)模式的广义力，fii为滑移(ii型)模式的广义力，fiii为剪开(iii型)模式的广义力，p1～p6分别为六个自由度方向上的节点力，r为转子轴横截面半径，a为裂纹深度，θ1和θ2分别为裂纹面空间方位角度，rb为在斜裂纹面上的投影半径，β为裂纹宽度方向上的积分变量，k为剪力修正系数，通常取k＝10/9，xs＝x+βcosθ1，x为裂纹面到转子轴左端面的距离变量，rb2＝(r2-β2sin2θ1)/sin2θ2。oiθ1θ2、siθ1θ2和tiθ1θ2分别为对应裂纹面张开(i型)、滑移(ii型)与剪开(iii型)三种位移模式的应力，其上标(·)θ1θ2表示裂纹面方位角θ1与θ2的共同作用，下标(·)i对应于单元节点力序号，以区别不同类型裂纹面上的应力作用。
步骤103：根据裂纹面张开、滑移与剪开三种位移模式的应力密度因子确定外加应变能；具体公式为：
其中，i＝1,2,3,4,5,6，e'＝e/(1-v)，ms＝1+v，v为泊松比，kiⅰ、kiⅱ与kiⅲ分别为描述裂纹面张开、滑移与剪开三种位移模式的应力密度因子，uc为外加应变能，e为杨氏弹性模量，a为转子轴横截面积，da为积分运算中对a求导数。
步骤104：根据无裂纹转子的应变能和外加应变能确定空间任意斜裂纹转子的总应变能模型，具体公式为：
u＝u0+uc；
其中，u0为无裂纹转子的应变能，uc为外加应变。
步骤20：根据所述总应变能模型确定空间任意斜裂纹转子的柔度参数，具体公式为：
其中，u0+uc为总应变能模型，u0为无裂纹转子的应变能，uc为外加应变能，ui为沿着第i个坐标方向的节点位移与节点力，pi为沿着第i个坐标方向的节点力，pj为沿着第j个坐标方向的节点力。
其中，f1、f2分别是张开(i型)模式的广义力，fii为滑移(ii型)模式的广义力，fiii为剪开(iii型)模式的广义力，r为转子轴横截面半径，a为裂纹深度，θ1和θ2分别为裂纹面空间方位角度，ms＝1+v为基于泊松比v所定义的新变量，rb为在斜裂纹面上的投影半径，β为裂纹宽度方向上的积分变量，k为剪力修正系数，通常取k＝10/9，xs＝x+βcosθ1，x为裂纹面到转子轴左端面的距离，l为转子轴长度，αs为剪切变形参数，g为刚度模量，i为横截面积惯性矩，e'＝e/(1-v)，ms＝1+v，v为泊松比，kiⅰ、kiⅱ与kiⅲ分别为描述裂纹面张开、滑移与剪开三种位移模式的应力密度因子，uc为外加应变能，e为杨氏弹性模量，a为转子轴横截面积，da为积分运算中对a求导数，rb2＝(r2-β2sin2θ1)/sin2θ2，i0为横截面积极惯矩。
步骤30：根据所述柔度参数确定空间任意斜裂纹转子的变刚度模型，具体公式为：
kc＝tg-1tt；
其中，g＝{gij},i＝1、2、…、6，j＝1、2、…、6，gij为空间任意斜裂纹转子的柔度参数，t为变换矩阵。
图2为本发明实施例基于空间任意斜裂纹转子的变刚度确定系统结构图，如图2所示，本发明还提供一种基于空间任意斜裂纹转子的变刚度确定系统，所述系统包括：
总应变能模型确定模块1，用于确定空间任意斜裂纹转子的总应变能模型；
柔度参数确定模块2，用于根据所述总应变能模型确定空间任意斜裂纹转子的柔度参数；
变刚度模型确定模块3，用于根据所述柔度参数确定空间任意斜裂纹转子的变刚度模型；
变刚度确定模块4，用于将待检测的空间任意斜裂纹转子的参数值输入所述变刚度模型确定空间任意斜裂纹转子的变刚度；所述待检测的空间任意斜裂纹转子的参数值为空间任意斜裂纹转子的长度l、直径d、裂纹深度a以及裂纹面空间方位角度θ1和θ2中至少一者。
下面对各个模块进行详细论述：
所述总应变能模型确定模块1，具体包括：
获取单元，用于获取无裂纹转子的应变能；
应力密度因子确定单元，用于确定裂纹面张开、滑移与剪开三种位移模式的应力密度因子；
外加应变能确定单元，用于根据裂纹面张开、滑移与剪开三种位移模式的应力密度因子确定外加应变能；
总应变能模型确定单元，用于根据无裂纹转子的应变能和外加应变能确定空间任意斜裂纹转子的总应变能模型。
外加应变能确定单元，用于根据裂纹面张开、滑移与剪开三种位移模式的应力密度因子确定外加应变能，具体公式为：
其中，e'＝e/(1-v)，ms＝1+v，v为泊松比，kiⅰ、kiⅱ与kiⅲ分别为描述裂纹面张开、滑移与剪开三种位移模式的应力密度因子，uc为外加应变能，e为杨氏弹性模量，a为转子轴横截面积，da为积分运算中对a求导数。
柔度参数确定模块2，用于根据所述总应变能模型确定空间任意斜裂纹转子的柔度参数，具体公式为：
其中，u0+uc为总应变能模型，u0为无裂纹转子的应变能，uc为外加应变能，ui为沿着第i个坐标方向的节点位移与节点力，pi为沿着第i个坐标方向的节点力，pj为沿着第j个坐标方向的节点力。
变刚度确定模块4，用于根据所述柔度参数确定空间任意斜裂纹转子的变刚度模型，具体公式为：
kc＝tg-1tt；
其中，g＝{gij},i＝1、2、…、6，j＝1、2、…、6，gij为空间任意斜裂纹转子的柔度参数，t为变换矩阵，
本发明获得空间任意斜裂纹转子的变刚度模型后，即可对不同类型裂纹的刚度变化特性进行对比，找出裂纹类型对裂纹转子刚度特性的影响规律。此外，还可以通过研究具有不同空间方位角度以及不同裂纹深度的空间任意斜裂纹的变刚度特性，明确这种特殊裂纹转子的刚度变化机理，并进一步通过动力学建模对具有空间任意斜裂纹转子的多自由度耦合振动特征进行分析，所得到的一系列分析结果可以为转子强度设计以及转子系统的故障诊断等应用奠定基础。
具体举例
转子轴分别采用铁木辛科梁单元建模，并考虑了纵向、弯曲及扭转六个方向的自由度。如图6为一段包含空间任意斜裂纹的转子轴，转盘质量为m，转轴直径为d，长度为l。裂纹位于左端x处，裂纹深度为α。单元受力情况如下：剪力p2、p3与p8、p9，弯矩p5、p6与p11、p12，轴向力p1与p7，扭矩p4与p10。
随着转子轴旋转，裂纹发生交替性开/合，导致裂纹转子刚度的周期性变化。为了仿真裂纹转子的变刚度特性，在一个旋转周期内(360°)将裂纹边均分为100个点，每个点确定一个ccl裂纹面闭合线位置，由此确定不同的积分限。从而，可将刚度参数变化视为是ccl位置的函数。本发明以铁木辛科梁单元建模转子轴为例子阐述求解空间任意斜裂纹转子变刚度特性的过程，其确定的空间任意斜裂纹转子的变刚度模型不仅取决于裂纹的类型，还与空间任意斜裂纹转子的长度l、直径d、裂纹深度a以及裂纹面空间方位角度θ1和θ2中有关。
实施例一
本发明令空间任意斜裂纹转子的有限元长度l＝0.7/14＝0.05m，裂纹深度比ā＝a/d＝0.3，a为裂纹深度，d为直径，裂纹面方位角θ1＝45°，θ2＝60°，通过改变裂纹的类型进行变刚度特性分析，即t-typecrack、ts-typecrack和as-typecrack分别为横向、横-斜向以及空间任意斜向裂纹，具体如图7所示，图7中t-typecrack型、ts-typecrack型和as-typecrack型裂纹的刚度参数分别以“-----”、“---”和“—”表示。
由图7可见，k11为纵向刚度值、k22为水平剪切方向刚度值、k33为垂直剪切方向刚度值、k44为扭转方向刚度值、k55为垂直弯曲方向刚度值、k66为水平弯曲方向刚度值。在转子整个旋转周期中，空间任意斜裂纹(as)转子相比于横裂纹和横-斜裂纹(t与ts)转子具有最大的轴向与垂直剪切刚度(k11与k33)。在旋转初期，as型裂纹转子的扭转刚度(k44)略小于ts型裂纹转子；在旋转末期，as型裂纹转子的弯曲刚度(k55与k66)略小于t型裂纹转子。不过，从整体上来看其刚度差不多也是最大的。但是，在转子旋转的大部分过程中as型裂纹转子的水平剪切刚度(k22)明显小于其他两类裂纹转子，虽然在旋转末期其刚度有明显提高。此外，as型与ts型裂纹转子的扭转刚度(k44)很接近，弯曲刚度(k55与k66)的变化趋势相近，表明在扭转和弯曲方向上裂纹面方位角θ1是影响转子刚度特性的主要因素。而在其他方向上，方位角θ2则发挥了重要的作用，导致as型裂纹转子具有与其他两类裂纹转子明显不同的、更加复杂的变刚度特性，尤其在水平剪切方向上(k22)表现得更加明显，刚度曲线的对称性被严重破坏，刚度值明显下降，整体小于其他两类裂纹转子。
实施例二
保持裂纹面方位角θ1＝45°不变，θ2在30°到90°之间变化。在不同方位角的共同作用下，as型裂纹转子的刚度特性曲线如图8所示。图8中分别以线型“-----”、“---”、“—”、“-+-”以及“-o-”按照θ2递增的顺序加以对比描述。需要注意的是，当θ2＝90°时as型裂纹退化为ts型裂纹。
由图8可见，随着θ2的增大，纵向刚度值(k11)、垂直剪切方向刚度值(k33)、扭转方向刚度值(k44)以及垂直弯曲方向刚度值(k55)均呈现单调下降趋势。但是，在水平弯曲方向刚度值(k66)以及水平剪切方向刚度值(k22)不存在单调下降趋势，刚度变化曲线彼此出现了明显交叉，表明as型裂纹转子的裂纹面方位角θ1与θ2之间的交互作用效应。显然，在转子整个旋转周期的不同转角位置，这种效应对转子刚度的影响是不同的，因此是动态的，也是非线性的。
实施例三
进一步地，保持as型裂纹转子的裂纹面方位角θ1＝45°，θ2＝60°不变，裂纹深度比ā从0.1变化到0.5。对应于不同裂纹深度的刚度特性变化曲线如图9所示，图9中分别以线型“-----”、“---”、“—”、“-+-”以及“-o-”按照ā递增的顺序进行描述。
由图9可见，随着ā的增大，as型裂纹转子各方向的刚度参数kii，i＝1～6均单调下降。在转子整个旋转周期的不同转角位置，刚度值是不断变化的，但是对应于不同裂纹深度的转子变刚度特性曲线几乎具有相同的走向和变化趋势。这表明，裂纹深度参数对转子刚度的影响虽然也是动态的，但是具有明显的线性特征。
对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。