一种齿轮裂纹扩展模拟的小波扩展有限元仿真分析方法

技术领域

本发明涉及针对裂纹扩展的有限元分析法，是一种采用了小波扩展这一新的有限分析单元的数值分析方法。 

背景技术

由于齿轮裂纹尖端区域位移与应力都含有 奇异性(其中r为裂纹尖端附近某一点的极坐标矢径)，在裂纹中奇异点附近的解梯度大，还会发生突变，因此给传统有限元方法求数值解造成困难。 

传统有限元方法采用连续函数作为形函数，要求在单元内部形状函数连续且材料性能不能跳跃，对于处理像裂纹这样的不连续问题时，需要将裂纹面设置为单元的边，裂尖点设置为单元的结点，而且要在裂尖点附近不连续体的奇异场内进行高密度网格划分以及在模拟裂纹扩展时需要不断的进行网格的重新划分，这使得有限元程序设计相当复杂，且效率极低。在有限元框架内发展起来的扩展有限元法(XFEM)，以解决不连续问题为着眼点，在处理裂纹问题是虽然不需要对结构内存在的几何或者物理界面进行网格剖分，但是对于工程中裂纹等奇异性问题的处理还有困难。将小波分析理论引入有限元分析的小波有限元方法因其具有的多尺度、多分辨的特性，可以解决因裂纹尖端区域奇异性带来的困难，但却不能反映裂纹的生长情况。 

发明内容

本发明为了方便有效地解决齿根裂纹尖端应力应变奇异性给传统有限元带来的困难，同时得到裂纹的实时生长情况，提出了一种采取小波扩展这一 新型有限分析单元的数值分析方法。本发明将两种有限元方法的优势结合，借鉴小波形函数更加准确的逼近特性和扩展有限元无需重新划分网格的特点，形成一种新的有限元数值方法。本发明不仅可以跟踪裂纹的生长状况，还能解决应力高度集中带来的困难，具有较高的计算精度和较高的计算效率。此外采用此方法模拟齿根裂纹扩展能够为机械设备故障诊断的研究带来方便。 

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

一种齿轮裂纹扩展模拟的小波扩展有限元仿真分析方法，该方法基于有限元软件ABAQUS和计算机辅助制图以及数学辅助计算软件构建的平台，对齿轮裂纹扩展进行分析，其特征在于包括以下步骤： 

1)建立含齿根裂纹的齿根裂纹模型： 

应用计算机辅助制图软件绘制给定齿轮模数和齿数的齿轮模型，将所建立的齿轮模型导入到有限元软件ABAQUS中自动划分网格生成含有H个单元的齿轮模型，H为正整数，并假设在齿根区域内L个单元内存在裂纹，L为自然数，且L＜H； 

2)获得齿根裂纹模型内各个单元的节点坐标： 

根据步骤1)中建立的齿根裂纹模型，应用有限元软件ABAQUS命令流输出功能输出齿根裂纹模型中H个单元的各个节点坐标； 

3)根据所建立的齿根裂纹模型计算出H个单元的小波单元刚度矩阵： 

将步骤2)齿根裂纹模型中所含有的H个单元的各个节点坐标数据导出，应用平面问题的势能泛函数得到小波单元刚度矩阵； 

小波单元刚度矩阵计算式： 

$k=\left(\begin{array}{cc}{K}^{e,1}& K^{e,2}\\ K^{e,3}& K^{e,4}\end{array}\right)---\left(1\right)$

分块子矩阵为： 

$K^{e,1}=\mathrm{Et}/(1-{\mu }^2)\left(\begin{array}{c}\left(\frac{1}{b-a}{\int }_0^1\frac{{\mathrm{d\Phi }}_1^T}{\mathrm{d\xi }}\frac{d{\Phi }_1}{\mathrm{d\xi }}\mathrm{d\xi }\right)\otimes \left((d-c){\int }_0^1{\Phi }_2^T{\Phi }_2\mathrm{d\eta }\right)+\\ \frac{(1-\mu )}{2}\left((b-a){\int }_0^1{\Phi }_1^T{\Phi }_1\mathrm{d\xi }\right)\otimes \left(\frac{1}{d-c}{\int }_0^1\frac{d{\Phi }_2^T}{\mathrm{d\eta }}\frac{d{\Phi }_2}{\mathrm{d\eta }}\mathrm{d\eta }\right)\end{array}\right)$

$K^{e,2}=\mathrm{Et}/(1-{\mu }^2)\left(\begin{array}{c}\left(\mu {\int }_0^1\frac{{\mathrm{d\Phi }}_1^T}{\mathrm{d\xi }}{\Phi }_1\mathrm{d\xi }\right)\otimes \left({\int }_0^1{\Phi }_2^T\frac{{\mathrm{d\Phi }}_2}{\mathrm{d\eta }}\mathrm{d\eta }\right)+\\ \frac{(1-\mu )}{2}\left({\int }_0^1{\Phi }_1^T\frac{d{\Phi }_1}{\mathrm{d\xi }}\mathrm{d\xi }\right)\otimes \left({\int }_0^1\frac{d{\Phi }_2^T}{\mathrm{d\eta }}{\Phi }_2\mathrm{d\eta }\right)\end{array}\right)$

K^e，3＝(K^e，2)^T

$K^{e,4}=\mathrm{Et}/(1-{\mu }^2)\left(\begin{array}{c}\left((b-a){\int }_0^1{\Phi }_1^T{\Phi }_1\mathrm{d\xi }\right)\otimes \left(\frac{1}{d-c}{\int }_0^1\frac{d{\Phi }_2^T}{\mathrm{d\eta }}\frac{d{\Phi }_2}{\mathrm{d\eta }}\mathrm{d\eta }\right)+\\ \frac{(1-\mu )}{2}(\frac{1}{b-a}{\int }_0^1\frac{d{\Phi }_1^T}{\mathrm{d\xi }}\frac{{\mathrm{d\Phi }}_1}{\mathrm{d\xi }}\mathrm{d\xi }+)\otimes \left((d-c){\int }_0^1{\Phi }_2^T{\Phi }_2\mathrm{d\eta }\right)\end{array}\right)$

其中a、b、c、d分别为平面单元求解域Ω_e＝{x，y|x∈[a，b]，y∈[c，d]}内x与y坐标取值的上下区间；x与y坐标值为步骤2)中所输出H个单元的各个节点坐标值；ξ，η为自然坐标，ξ＝(x-a)/(b-a)，η＝(y-c)/(d-c)，则标准求解域Ω_s＝{ξ，η|ξ，η∈[0，1]}；E为弹性模量；μ为泊松比；t为厚度；以上各分块子矩阵中Φ₁和Φ₂分别为m阶j尺度下的一维区间B样条小波尺度函数向量(详细请参看何正嘉陈雪峰著《小波有限元理论及其工程应用》科学出版社第36、171页)； 

4)应用扩展有限元对齿根裂纹模型中齿根区域内含有齿根裂纹的L个单元的小波单元位移场函数 进行改进，计算公式如下： 

$\left(\begin{array}{c}u(\xi ,\eta )\\ v(\xi ,\eta )\end{array}\right)=\underset{s}{\Sigma }{N}_s(\xi ,\eta )\{\left(\begin{array}{c}{u}_s\\ v_s\end{array}\right)+H_s(\xi ,\eta )\left(\begin{array}{c}{a}_{s1}\\ a_{s2}\end{array}\right)\}---\left(2\right)$

其中ξ，η为自然坐标，ξ＝(x-a)/(b-a)，η＝(y-c)/(d-c)，ξ，η取 值范围为齿根裂纹模型中含有齿根裂纹的L个单元内部区域的自然坐标值；s代表含有齿根裂纹的L个单元内各个节点编号；N_s(ξ，η)代表含有齿根裂纹的L个单元内部区域对应的小波形函数，N_s(ξ，η)＝Φ(R^e)^-1，其中 Φ₁和Φ₂与步骤3)中相同，分别为m阶j尺度下的一维区间B样条小波尺度函数向量；R^e为单位矩阵；u_s、v_s表示含有齿根裂纹的L个单元内各个节点的常规自由度；a_s1和a_s2为含有齿根裂纹的L个单元内每一个节点的附加自由度；H_s(ξ，η)为反映含有齿根裂纹的L个单元的位移不连续性函数，计算公式如下： 

sign(x)为符号函数； 为符号距离函数用于描述含有齿根裂纹的L个单元内任意点与齿根裂纹之间的垂直距离，计算公式如下； 

为含有齿根裂纹的L个单元内各个节点处的 值，(ξ_s，η_s)代表含有齿根裂纹的L个单元内各个节点的自然坐标值； 

5)应用扩展有限元改进齿根裂纹模型中齿根区域内含有齿根裂纹的L个单元的的小波单元位移场函数之后，计算齿根区域内含有齿根裂纹的L个单元的小波扩展单元刚度矩阵(详细请参看：卡坦著韩来彬译《MATLAB有限元分析与应用》清华大学出版社第十三章)，计算公式如下： 

$K^e=\underset{{\Omega }_s}{\int }t{B}^{\mathrm{eT}}D{B}^e\mathrm{Jd\xi d\eta }---\left(5\right)$

t为厚度；B^e为应变矩阵(参考文献杨珺平面裂纹扩展分析的扩展有限 元法[D].南京航空航天大学航空宇航学院，2007.第31页)；D为弹性矩阵， $D=\frac{E}{1-{\mu }^2}\left(\begin{array}{ccc}1& \mu & 0\\ \mu & 1& 0\\ 0& 0& \frac{1-\mu }{2}\end{array}\right);$E为弹性模量；μ为泊松比；J为雅可比矩阵(详细请参看冷纪桐赵军著《有限元技术基础》化学工业出版社第84页)；Ω_s为标准求解域Ω_s＝{ξ，η|ξ，η∈[0，1]}； 

6)计算齿根裂纹模型的结构刚度矩阵； 

根据步骤3)、4)、5)对齿根裂纹模型中H个小波单元刚度矩阵的计算以及齿根区域内含有齿根裂纹的L个单元的小波扩展单元刚度矩阵的计算，按照齿轮裂纹模型中H个单元的各个节点次序集成结构刚度矩阵K(详细请参看胡于进王璋奇著《有限元分析及应用》清华大学出版社第59页)。 

7)建立有限元方程求解节点常规自由度和附加自由度的数值解： 

基于步骤6)中齿根裂纹模型结构刚度矩阵K的计算，建立以下有限元求解方程： 

$K\left(\begin{array}{c}u\\ a\end{array}\right)=\left(\begin{array}{c}{P}_u\\ P_a\end{array}\right)---\left(6\right)$

式中u、a分别为齿根裂纹模型中各个节点的常规自由度与附加自由度向量；P_u、P_a分别为与模型中常规自由度向量和附加自由度向量相对应的载荷向量，解方程可求出各个节点常规自由度和附加自由度的数值解； 

8)计算齿根裂纹模型中裂纹的张开位移和裂纹开裂角度： 

基于步骤7)求解有限元方程得到的齿根裂纹模型中附加自由度的数值解，可求齿根裂纹模型中裂纹张开位移，计算式如下： 

$w=\stackrel{v}{n}·2\underset{s\in n_s}{\Sigma }{N}_s(\xi ,\eta )a_s---\left(7\right)$

w代表裂纹张开位移值； 为裂纹上任一点的法向量；s代表含有齿根裂纹的L个单元内各个节点编号；n_s为齿根区域内含有齿根裂纹的L个单元包含的所有改进节点的数量；N_s(ξ，η)代表含有齿根裂纹的L个单元内部区域对应的小波形函数，计算裂纹张开位移时N_s(ξ，η)自变量取值为齿根裂纹面上，即当 时对应的数值；a_s为含有齿根裂纹的L个单元内每一个节点的附加自由度向量，每个该向量中含有步骤4)中所述a_s1和a_s2两个附加自由度值； 

裂纹的裂纹开裂角度θ₁可由下式计算： 

$\left(\begin{array}{cc}{\theta }_1={\tan }^{-1}\frac{{\sigma }_1-{\sigma }_x}{{\tau }_{\mathrm{xy}}}+\pi /2& {\tau }_{\mathrm{xy}}\ne 0\\ {\theta }_1={\tan }^{-1}\frac{{\tau }_{\mathrm{xy}}}{{\sigma }_1-{\sigma }_x}+\pi /2& {\sigma }_1-{\sigma }_x\ne 0\end{array}\right)---\left(8\right)$

式中σ₁为单元的最大主应力，计算公式如下： 

${\sigma }_1=\frac{{\sigma }_x+{\sigma }_y}{2}+\sqrt{\frac{{\sigma }_x-{\sigma }_y}{2}+{{\tau }_{\mathrm{xy}}}^2}---\left(9\right)$

以上两式中σ_x、σ_y、τ_xy分别表示单元应力在横向、纵向以及切向的数值分量(详细请参看冷纪桐赵军著《有限元技术基础》化学工业出版社第53页)，根据步骤7)中求解有限元方程得到的齿根裂纹模型中各个节点常规自由度数值解计算； 

以上所述为应用小波扩展有限元方法模拟齿根裂纹扩展的仿真方法，本发明的流程图如图1所示。 

本发明的有益效果是：利用计算机辅助制图比如Pro/E在齿轮建模方面的优势，建立精确的齿轮啮合模型，避免了有限元软件在建模，特别是其在有 复杂裂纹故障的非正常齿轮建模时的缺憾；再结合ABAQUS强大的公式计算和数据输出功能，获得所划分网格的坐标数据；应用数学辅助计算软件如MATLAB编制程序计算出小波扩展这种单元的结构刚度矩阵，最后引入边界条件和载荷求解有限元方程，得到齿根裂纹扩展的模拟结果，还可根据需要改变差值函数的阶次获得不同的求解精度，这为齿轮系统故障的机理研究提供了准确、可靠的理论基础。 

附图说明

图1是本发明的流程图； 

图2是本实施例工作流程示意图； 

图3是本实施例齿轮副二维实体模型示意图； 

图4是本实施例含齿根裂纹轮齿二维实体模型示意图； 

图5是本实施例在ABAQUS中划分网格后的齿轮有限元模型示意图； 

图6是本实施例要改进的单元1与单元10的示意图； 

图7是本实施例所中单元10对应的数值解； 

图8是本实施例裂纹张开位移随裂纹长度变化的曲线； 

图9是是本实施例模拟的裂纹扩展效果图； 

具体实施方式

下面具体结合附图与实例对本发明作进一步的说明。 

如图2所示，是本实施例工作流程示意图。主要由3大部分组成，详细的步骤如下。 

一、利用Pro/E绘图软件绘制模数为2，齿数为19的齿轮与另一大齿轮啮合模型如图3所示。再根据需求可在轮齿上添加任意形式的裂纹故障，可 以改变裂纹的位置，大小，形态等特征，仿真出更符合实际情况的复杂裂纹故障。图4所示为假设一个在齿根处有图示裂纹故障的轮齿。 

二、利用ABAQUS划分有限元网格并输出网格节点坐标： 

1)将用Pro/E绘制的齿轮模型导入ABAQUS并划分生成含有22个有限单元的齿轮模型；假设裂纹故障存在于单元1与单元10中，为方便计算，实例中仅用到了含有裂纹故障的一个齿而且采用了较为粗糙的网格划分，利用有限元软件的命令显示各个单元编号和节点编号，如图5所示为划分网格后显示单元和节点编号的齿轮有限元模型； 

2)在进行步骤1)的同时，ABAQUS会根据菜单操作的步骤自动生成对应的命令流文件，此命令流文件包含了菜单操作时的所有信息。本发明中，利用ABAQUS的计算功能和结果数据写入文本命令，即可得到整个结构中各个节点的几何坐标数据。 

三、应用MATLAB计算所有22个单元的单元刚度矩阵以及齿根裂纹模型的结构刚度矩阵并建立有限元方程： 

1)首先将所得节点坐标数据导入到小波单元刚度矩阵计算式 $k=\left(\begin{array}{cc}{K}^{e,1}& K^{e,2}\\ K^{e,3}& K^{e,4}\end{array}\right),$应用数学软件MATLAB编程计算得到全部22个小波单元刚度矩阵； 

2)对含有齿根裂纹的两个小波单元位移场即单元1与单元10进行改进，如图6所示为需要改进的单元1与单元10的示意图，图中斜线表示裂纹，方形符号表示经过改进的裂纹贯穿单元节点，方程η＝0.8ξ+0.1为裂纹在自然坐标下的表达式，(0，0.3)、(1，0.9)为裂纹与单元10的边界交点坐标，改进后 的单元位移场表达式为： 

$\left(\begin{array}{c}u(\xi ,\eta )\\ v(\xi ,\eta )\end{array}\right)=\underset{s}{\Sigma }{N}_s(\xi ,\eta )\{\left(\begin{array}{c}{u}_s\\ v_s\end{array}\right)+H_s(\xi ,\eta )\left(\begin{array}{c}{a}_{s1}\\ a_{s2}\end{array}\right)\}$

其中(ξ，η)代表单元1与单元10内部区域的自然坐标值；s代表单元1与单元10内各个节点编号；H_s(ξ，η)反映齿根裂纹的位移不连续性函数；改进的主要内容是在每一个原有节点常规自由度的基础上加入了一个附加自由度；含有齿根裂纹的单元1与单元10内每一个节点的附加自由度为a_s1和a_s2； 

3)然后将扩展有限元对单元位移场的改进引入到单元刚度矩阵的积分方程 中得到相应的两个小波扩展单元刚度矩阵，t为厚度；B^e为改进单元位移场后所得到的应变矩阵(参考文献杨珺平面裂纹扩展分析的扩展有限元法[D].南京航空航天大学航空宇航学院，2007.第31页)，弹性矩阵 $D=\frac{E}{1-{\mu }^2}\left(\begin{array}{ccc}1& \mu & 0\\ \mu & 1& 0\\ 0& 0& \frac{1-\mu }{2}\end{array}\right);$E为弹性模量；μ为泊松比；J为雅可比矩阵；Ω_s＝{ξ，η|ξ，η∈[0，1]}为标准求解域。 

4)计算完成所有22个单元的刚度矩阵之后，根据结构的单元和节点排列次序构造出整体的结构刚度矩阵K(详细请参看胡于进王璋奇著《有限元分析及应用》清华大学出版社第59页)。 

5)引入边界约束条件和载荷，单元1和单元10含有初始裂纹，限制齿轮两侧及底边节点的自由度，并在节点4施加载荷求解。基于本实施例步骤4)中齿根裂纹模型结构刚度矩阵K的计算，建立以下有限元求解方程： 

$K\left(\begin{array}{c}u\\ a\end{array}\right)=\left(\begin{array}{c}{P}_u\\ P_a\end{array}\right)$

式中u、a分别为齿根裂纹模型中各个节点的常规自由度与附加自由度向量：u＝[u₁ v₁ L u_r v_r]^T  a＝[a₁₁ a₁₂ L a_r1 a_r2]^T，模型节点总数为r＝33；u₁ v₁、u_r v_r分别代表模型中第一个节点常规自由度值和最后一个节点常规自由度值；a₁₁ a₁₂、a_r1 a_r2分别代表模型中第一个节点的附加自由度值和最后一个节点的附加自由度值；P_u、P_a为与模型中常规自由度向量和附加自由度向量相对应的载荷向量，利用MATLAB解以上方程可得到齿根裂纹模型中所有节点常规自由度和附加自由度的数值解； 

6)根据5)中所求节点附加自由度数值解计算单元1和单元10中裂纹的张开位移，以单元10为例， 为自然坐标下裂纹上的任意法向量，节点29、21、3、25的小波形函数值分别为N₁＝0.75、N₂＝0.05、N₃＝0.85、N₄＝0.25，各节点的附加自由度值分别为(-0.1714，0.0375)、(-0.2013，-0.1112)、(-0.0566，0.0030)、(-0.0656，-0.0826)，将以上数据代入公式 中即可得到如图7所示扩展位移数值。 

由此例可知，附加自由度是与常规自由度类似的节点位移，在每一个常规自由度上增加一个附加自由度，从而可以反映裂纹的张开位移，这就是对含有齿根裂纹的单元1和单元10的小波单元位移场改进之后所产生的效果；图7为根据齿根单元10中常规自由度和附加自由度求得的裂纹张开位移，图8为裂纹张开位移随裂纹长度的变化曲线； 

7)根据5)中所求节点常规自由度数值解计算单元22中的单元应力在横向、纵向及切向上的数值分量σ_x、σ_y、τ_xy(详细请参看冷纪桐赵军著《有限元技术基础》化学工业出版社第53页)；然后根据下列公式计算裂纹在单元22中的扩展角： 

$\left(\begin{array}{cc}{\theta }_1={\tan }^{-1}\frac{{\sigma }_1-{\sigma }_x}{{\tau }_{\mathrm{xy}}}+\pi /2& {\tau }_{\mathrm{xy}}\ne 0\\ {\theta }_1={\tan }^{-1}\frac{{\tau }_{\mathrm{xy}}}{{\sigma }_1-{\sigma }_x}+\pi /2& {\sigma }_1-{\sigma }_x\ne 0\end{array}\right)---\left(8\right)$

式中σ₁为单元的最大主应力，计算公式如下： 

${\sigma }_1=\frac{{\sigma }_x+{\sigma }_y}{2}+\sqrt{\frac{{\sigma }_x-{\sigma }_y}{2}+{{\tau }_{\mathrm{xy}}}^2}$

得到裂纹的扩展角度为θ₁＝32.46度；裂纹在单元22中的扩展方向以及裂纹扩展的整个过程如图9所示。