当代战争要求武器装备在战备状态及应急突发事件时，能够以不惜牺牲部分发动机的热疲劳寿命来实现最短时间内的快速启动。而这种不计代价的快速启动方式，将使舰船发动机的涡轮转子无法经历暖机过程，而在最短的时间内快速达到运行工况，进入流道的气流温度与转子金属表面间存在较大的温度梯度，导致转子承受着巨大的热冲击。而由于涡轮盘常工作在复杂非均匀的热载荷环境下，极易产生热疲劳积累而导致转子总应力增加及材料性能恶化^[1]。因此，对舰船发动机涡轮转子的热分析应侧重热冲击及非对称热载荷分布下的瞬态热应力变化及温度分布特性研究。

国内外学者关于发动机转子的热应力分析侧重理论计算^[2-4]和在线监测^[5-8]两方面。而对于转子的瞬稳态温度场、热应力计算方面的研究主要是建立理论解析模型^[9-11]和有限元数值仿真^[12-14]，其中理论解析模型正是通过求解热传导方程和弹性理论控制方程得到热变形及热应力，并建立普遍形式的转子热应力求解方程，这样能够较为简便地给出温度分布下的热变形及热应力，多适用于在线监测中热疲劳寿命的快速预测分析。而有限元数值仿真是基于数值法解决基本传热方程和弹性理论，它是一种能够有效准确地分析各种复杂结构和边界条件的方法。回顾近年来研究人员在进行发动机转子的热分析来看，有限元法数值仿真结果与在线实测结果也更加接近^[15]。而其缺点就是计算工作量较大、需要验证仿真结果的精确性。

事实上，瞬态热冲击载荷是影响转子寿命的关键，而针对瞬态热冲击下发动机转子典型结构参数及热冲击载荷边界对热应力和温度场的影响研究却不多。同时，针对发动机转子在实际工作中受到的周向非对称热载荷的问题，如燃烧室火焰筒点火不均匀、熄火等故障，涡轮盘将承受局部热冲击载荷，极易导致转子发生热变形，最终严重影响转子热疲劳寿命并产生热致振动现象，这方面研究也较缺乏。因此，本文基于热弹耦合理论建立了转子热冲击理论计算及有限元分析模型，研究周向均布和非对称两种热冲击载荷下涡轮转子的最大热应力变化和温度场分布特性，为舰船发动机方案设计初期的热疲劳寿命预估提供有意义的借鉴。

1 转子主轴结构—热分析模型

转子主轴的瞬态热分析建模，将采用基于热传导方程和弹性理论推导出的理论方法。假设转子主轴外表面温度为T_b，初始温度分布为T₀，转子外径大小为b，假设温度分布沿轴向及周向均匀分布，即温度T函数仅是半径r的函数，如图 1所示。

1.1 主轴结构热传导分析

上述问题可以描述为一维径向导热问题，引入过余温度函数：

$\theta \left( r,\tau \right)=T\left( r,\tau \right)-{{T}_{b}}$

(1)

圆柱坐标系下无热源的热传导控制方程为

$\frac{\partial {{T}^{2}}}{\partial {{r}^{2}}}+\frac{1}{r}\frac{\partial T}{\partial r}+\frac{1}{{{r}^{2}}}\frac{{{\partial }^{2}}T}{\partial {{\theta }^{2}}}+\frac{{{\partial }^{2}}T}{\partial {{z}^{2}}}=\frac{1}{a}\frac{\partial T}{\partial t}$

(2)

式中:a为热扩散率，m²/s。引入过余温度函数的一维导热问题，式(2) 可以简化为

$\frac{1}{a}\frac{\partial \theta }{\partial \tau }=\frac{{{\partial }^{2}}\theta }{\partial {{r}^{2}}}+\frac{1}{r}\frac{\partial \theta }{\partial r}$

(3)

采用过余温度函数表示的第一类温度边界以及条件初始条件可以表示为

$\theta \left( b,\tau \right)=0,\theta \left( r,0 \right)={{T}_{0}}-{{T}_{b}}$

(4)

运用数学物理方法中分离变量的方法，可将θ(r,τ)分解成关于时间项及空间项特征函数：

$\theta \left( r,\tau \right)=R\left( r \right)T\left( \tau \right)$

(5)

其中，时间项特征函数特解形式为

$T\left( t \right)={{\text{e}}^{\text{-}a{{\lambda }^{2}}\tau }}$

(6)

空间项特征函数解的表达式为

$R\left( {{\lambda }_{n}},r \right)=\sum\limits_{n=1}^{\infty }{{{A}_{\text{n}}}{{\text{J}}_{\text{0}}}}\left( {{\lambda }_{n}},r \right){{\text{e}}^{\text{-}a{{\lambda }^{2}}\tau }}$

(7)

因此过余温度函数可以整理成：

$\theta \left( r,t \right)=\sum\limits_{n=1}^{\infty }{{{A}_{n}}}{{\text{J}}_{\text{0}}}\left( {{\lambda }_{n}}r \right){{\text{e}}^{\text{-a}{{\lambda }^{2}}t}}$

(8)

联立式(4) 、(8) ，可以得到关于特征值λ_n以及系数A_n的方程：

${{\text{J}}_{\text{0}}}\left( {{\lambda }_{n}},b \right)=0,\text{ }{{A}_{\text{n}}}=\frac{2\left( {{T}_{0}}-{{T}_{b}} \right)}{\left( {{\lambda }_{n}},b \right){{\text{J}}_{1}}\left( {{\lambda }_{n}},b \right)}$

(9)

故转子主轴结构瞬态温度场函数的表达式为

$\frac{T\left( r,t \right)-{{T}_{b}}}{{{T}_{0}}-{{T}_{b}}}=2\sum\limits_{n=1}^{\infty }{\frac{{{\text{e}}^{\text{-}a{{\lambda }^{2}}\tau }}{{\text{J}}_{\text{0}}}\left( {{\lambda }_{n}},r \right)}{\left( {{\lambda }_{n}},b \right){{\text{J}}_{1}}\left( {{\lambda }_{n}},b \right)}}$

(10)

式中：J_nλ_nr为贝塞尔函数，λ_n为贝塞尔函数特征值。

1.2 主轴结构热弹耦合分析

由于转子的轴向尺寸比另外两向要大得多，故可等效为平面应变问题。根据平面应变理论有

圆柱坐标系下的运动平衡微分方程：

$\frac{d{{\sigma }_{rr}}}{\text{d}r}+\frac{{{\sigma }_{rr}}-{{\sigma }_{00}}}{r}=0$

(11)

几何方程：

${{\sigma }_{rr}}=\frac{du}{\text{d}r},{{\varepsilon }_{00}}=\frac{u}{r}$

(12)

物理方程：

$\left\{ \begin{matrix} {{\varepsilon }_{rr}}=\frac{1}{E}\left[ {{\sigma }_{rr}}-\mu \left( {{\sigma }_{00}}+{{\sigma }_{zz}} \right) \right]+\alpha \left( T-{{T}_{0}} \right) \\ {{\varepsilon }_{00}}=\frac{1}{E}\left[ {{\sigma }_{00}}-\mu \left( {{\sigma }_{rr}}+{{\sigma }_{zz}} \right) \right]+\alpha \left( T-{{T}_{0}} \right) \\ {{\sigma }_{zz}}=\mu \left( {{\sigma }_{rr}}+{{\sigma }_{00}} \right)-E\alpha \left( T-{{T}_{0}} \right) \\ \end{matrix} \right.$

(13)

联立式(10) 、(13) ，得到转子柱坐标系下瞬态热应力分量如下

$\begin{align} & {{\sigma }_{rr}}=\frac{E\alpha }{1-v} \\ & \left\{ \frac{2\left( {{T}_{0}}-{{T}_{b}} \right)}{{{b}^{2}}}\sum\limits_{n=1}^{\infty }{\frac{1}{{{\lambda }_{n}}^{2}}{{\text{e}}^{\text{-}a\lambda _{n}^{2}t}}-\frac{2\left( {{T}_{0}}-{{T}_{b}} \right)}{rb}\sum\limits_{n=1}^{\infty }{\frac{1}{{{\lambda }_{n}}^{2}}\frac{{{\text{J}}_{\text{1}}}\left( {{\lambda }_{n}},r \right)}{{{\text{J}}_{\text{1}}}\left( {{\lambda }_{n}},b \right)}{{\text{e}}^{\text{-}a\lambda _{n}^{2}t}}}} \right\} \\ \end{align}$

(14)

$\begin{align} & \left\{ \frac{2\left( {{T}_{0}}-{{T}_{b}} \right)}{{{b}^{2}}} \right.\sum\limits_{n=1}^{\infty }{\frac{{{\text{e}}^{\text{-}a\lambda _{n}^{2}\tau }}}{{{\lambda }_{n}}^{2}}+\frac{2\left( {{T}_{0}}-{{T}_{b}} \right)}{rb}\cdot \sum\limits_{n=1}^{\infty }{\frac{{{\text{e}}^{\text{-}a\lambda _{n}^{2}\tau }}}{{{\lambda }_{n}}^{2}}\frac{{{\text{J}}_{\text{1}}}\left( {{\lambda }_{n}},r \right)}{{{\text{J}}_{\text{1}}}\left( {{\lambda }_{n}},b \right)}-}} \\ & \left. \sum\limits_{n=1}^{\infty }{\frac{2\left( {{T}_{0}}-{{T}_{b}} \right){{\text{e}}^{\text{-}a\lambda _{n}^{2}\tau }}{{\text{J}}_{0}}\left( {{\lambda }_{n}},r \right)}{\left( {{\lambda }_{n}},b \right){{\text{J}}_{\text{1}}}\left( {{\lambda }_{n}},b \right)}} \right\} \\ \end{align}$

(15)

$\begin{align} & {{\sigma }_{zz}}=\frac{E\alpha }{1-v} \\ & \left\{ \frac{4v}{{{b}^{2}}} \right.\left( {{T}_{0}}-{{T}_{b}} \right)\sum\limits_{n=1}^{\infty }{\frac{{{\text{e}}^{\text{-}a\lambda _{n}^{2}\tau }}}{{{\lambda }_{n}}^{2}}+\left( {{T}_{b}}-{{T}_{0}} \right)\left( v-1 \right)-} \\ & \frac{2\left( {{T}_{0}}-{{T}_{b}} \right)}{{{b}^{2}}}\sum\limits_{n=1}^{\infty }{\frac{1}{{{\lambda }_{n}}}\frac{{{\text{J}}_{0}}\left( {{\lambda }_{n}},r \right)}{{{\text{J}}_{\text{1}}}\left( {{\lambda }_{n}},b \right)}}{{\text{e}}^{\text{-}a\lambda _{n}^{2}\tau }} \\ \end{align}$

(16)

而米塞斯应力常作为衡量透平机械强度的准则，可通过下式求解：

$\sigma =\sqrt{{{\left( {{\sigma }_{rr}}-{{\sigma }_{\varphi \varphi }} \right)}^{2}}+{{\left( {{\sigma }_{rr}}-{{\sigma }_{zz}} \right)}^{2}}+{{\left( {{\sigma }_{\varphi \varphi }}-{{\sigma }_{zz}} \right)}^{2}}}$

(17)

1.3 有限元仿真分析

为验证转子瞬态温度场及热应力场的理论模型，对比了有限元数值仿真结果。主轴模型共划分400个温度—位移耦合单元，并完成了网格密度及时间步长的无关性验证，二者之间的误差在5%以内，证明了本文轴单元热分析模型的准确性。

图 2为瞬态热冲击载荷下圆柱体主轴的瞬态径向、周向和轴向热应力，以及瞬态温度场的变化情况，同时对比了有限元计算结果。

图 2(a)所示径向热应力在初始阶段急剧增大，最大径向热应力随半径发生变化，随着半径的增加，最大径向热应力逐渐减少。在达到最大值后，径向热应力随时间的增加而减小到0。如图所示，对比在r=0，r=0.5b，r=0.75b和r=b处的径向应力两者之间的平均径向应力误差在1.03%，式(14) 的理论结果与有限元结果是一致的。由于式(14) 中半径在分母上，r=0理论上是径向热应力分析的奇点线，但有限元可以给出结果。在r=b处，式(14) 预测的结果比有限元结果小，其他位置公式二者拟合的很好。

图 2(b)给出了周向热应力分量。在中心附近位置是拉应力，而在轴外表面附近变成了压应力。式(15) 预测周向热应力要大于有限元法，两者的误差不大于2.35%。半径处由于开始阶段外表面温度梯度较大，在初始阶段对外表面的网格质量要求高，而有限元模型结果受限于轴表面处的网格尺寸，而理论方法却没有这个问题。随着轴表面单元的径向尺寸减小，有限元结果更接近于理论解。

轴向热应力分布如图 2(c)所示。轴向热应力的分布规律与周向相似。轴向热应力随时间逐渐接近一个定常的压应力，周向热应力逐渐变成0。对比式(16) 的理论结果和有限元结果，两种方法的轴向热应力结果一致，大部分位置误差不大于0.5%。式(16) 和有限元法的差别发生在初始阶段的外表面。误差在于轴表面单元的径向尺寸。

图 2(d)比较了两种方法的瞬态温度分布结果。径向温度分布开始迅速变化，这是由内外径巨大的温差引起的。径向温度分布随时间变化逐渐减慢。几分钟后，温度分布恒定。式(10) 的结果与有限元结果保持一致，两者之间的误差不大于1%，证实了理论模型的正确性。

图 2(e)比较了式(17) 和有限元法计算的米塞斯应力，由于在位置r=0和r=b处的应力分量的不同，米塞斯应力的差别在这两处比较明显。对于理论分析来说，轴线处是奇异的，理论方法无法求解。但轴表面的理论结果要比有限元法更精准，特别在瞬态传热径向温度梯度很明显的初始阶段。为了仿真瞬态热应力更准确，有限元法需要在轴表面处更优良的网格质量。

通过比较理论解与有限元仿真结果可以发现，圆柱体主轴的瞬态温度场分布及三维瞬态热应力变化规律基本一致，验证了理论解的准确性。

2 结构参数对轮盘根部倒圆结构的影响热分析

轮盘根部倒角是由主轴和轮盘构成的，是涡轮转子最普遍的一种结构形式。一般来讲，倒角处的热应力要远远大于其他位置，因此将采用有限单元模型对倒角处温度分布及热应力进行数值仿真析。

倒角结构如图 3所示，倒圆角是由轴半径b和盘半径D形成的，倒角的半径为R，倒角结构的初始温度为T₀，外表面的边界温度是T_b。理论上，由于倒圆角结构的奇异性，很难通过理论分析获得其热应力。故该处热应力将通过有限元法进行研究。材料属性与之前的研究保持一致。由于倒角处的热应力明显大于其他部位，因此该处的网格应更优良。分析前验证了网格密度和时间步长的无关性，二者差别不大于1.8%，证明了分析模型的准确性。倒角网格局部细化模型如图 4所示。

倒圆半径R=1 mm的模型，轮盘根部倒圆角处的最大瞬态热应力值随时间变化趋势如图 5。可以看到在受到瞬态热冲击后，热应力在初始阶段快速增加，转子的最大瞬态热应力发生在44 s左右，峰值过后热应力逐渐减小并趋于稳定值。

倒圆半径明显的影响最大热应力，取一系列同轴径、盘径，但不同倒圆角半径(1 mm~30 mm)的模型进行研究和比较，如图 6(a)所示。不同倒圆半径之间的热应力相差比较明显，最大热应力从2 338 MPa(R=1 mm)减小到731.5 MPa(R=30 mm)。增加倒圆半径，是减小倒角处热应力的有效方法。

将轴径与盘径之比简化成β=D/b，研究不同倒圆角R下的热应力随β的变化规律。结果如图 6(b)，对相同倒圆半径下的模型，在盘轴径比β＜1.4时，热应力随β而增加。当盘轴径比β>1.4时，热应力值逐渐降低到常值。不同模型的热应力在不同倒角时具有同样的变化趋势。

当盘轴径比在β=1.5时，倒圆半径为R=1 mm,3mm,5mm,7mm,10mm热应力随轴径的变化规律如图 6(c)。随轴径的增加，热应力几乎不变。不同倒圆半径具有相同的结果。

图 6(d)给出了不同轮盘厚度下的热应力,当轴径b=50mm，盘轴径比β=1.5，倒圆半径R=5mm时，研究轮盘厚度对倒角最大热应力的影响。倒角最大热应力随轮盘厚度的增加而增大。

图 6(e)给出了不同盘轴径之比下最大热应力发生时间随倒圆半径的变化。可以看到发生时间随盘轴径之比增加而增加，同时增加变得的越来越慢。在同样的盘轴径之比下，发生时间随小β值增加而缓慢增加。然而在大β值下增加的更加明显。在盘轴径之比不变下，最大热应力发生时间随倒圆半径的变化。随倒圆半径的增加，发生时间急剧增加。倒圆半径越大，出现时间越久。

综上，倒圆角半径对转子最大热应力的影响较大，随倒圆角的增加最大热应力将会大幅度降低；轮盘与轴的半径比对最大热应力有一定的影响；主轴的半径对最大热应力的影响不大；轮盘厚度的增加将会增大最大热应力。

3 涡轮转子非对称热冲击研究

燃烧筒点火不均匀及运行时出现故障，会导致燃机涡轮盘承受瞬态非对称热载荷。这类问题属于周向非对称边界问题，即在转子某一扇区(1/6，2/6，…，5/6) 加载瞬态热冲击载荷，其与周围扇区的温差骤升，在此情况下对转子进行瞬态热冲击分析。研究对象取经过适当简化后的某型燃机涡轮端转子。初始条件及边界条件用数学模型可表示为

$\left\{ \begin{array}{*{35}{l}} T\left( r,0 \right)={{T}_{0}}={{0}^{\circ }}\text{C} \\ T\left( r,\theta ,\tau \right)\left| _{r={{R}_{2}}}={{T}_{b}}={{200}^{\circ }}\text{C},{{0}^{\circ }}\le \theta \le {{300}^{\circ }} \right. \\ T\left( r,\theta ,\tau \right)\left| _{r={{R}_{2}}}={{T}_{0}}={{0}^{\circ }}\text{C},{{300}^{\circ }} ＜\theta \le {{360}^{\circ }} \right. \\ U(\left. x,y,z \right|=0 \\ \end{array} \right.$

(18)

通过对涡轮转子表面不同扇区大小进行等温差热冲击，发现转子最大热应力发生位置位于涡轮盘热冲击扇区表面与轴交界的倒角处，在一分钟内轮盘最大热应力水平随时间的变化规律如表 1。

如图 7，轮盘热应力值在前十秒快速增长并达到极值，而后缓慢降低，其中5/6扇区的升高速率最快且热应力值也最大，而后50 s最大应力逐渐降低，热冲击扇区的角度越大，对涡轮的强度越不利。

图 8表明：热冲击扇区越大，温度场、应力场分布越趋于均匀，且平均应力水平逐渐增大，对转子安全性的影响增大；同时，热应力值较大的危险区域主要集中在轮盘与轴结合的倒角的一定范围内。

4 结论

本文基于导热控制方程以及弹性力学对涡轮转子两种典型结构的瞬态热应力进行了研究，主要结论包括：

1) 瞬态最大热应力水平远远高于稳态结果，倒圆半径越大，瞬态热应力降低越明显，轮盘厚度及盘轴比会影响最大热应力的大小，而主轴半径对其几乎不影响。

2) 瞬态热冲击下，涡轮转子的最大热应力出现时间主要与主轴半径有关，盘轴比对其略微有影响。

3) 涡轮转子瞬态在非对称热冲击下的最大热应力开始快速增大并达到极值，而后缓慢减小，非对称热冲击的扇区越大，应力极值越大。

热冲击下的瞬态热应力分析结果对涡轮转子结构设计初期评估热疲劳寿命提供有意义的参考。