2018, Vol. 39
叶轮作为离心压缩机的核心部件，在实际运行中因受迫振动导致的疲劳破坏可能造成严重后果。为展开离心叶轮受迫响应的计算分析，应当获得叶轮部件的完整力学参数，其中气动阻尼是较难得到的关键数据。同时，气动阻尼亦是气弹稳定性预测的关键参数，在先进叶轮机械的设计阶段，获得气弹稳定性裕度是受到关注的重要问题。对于不带冠离心叶轮，因叶片和轮盘是一体加工的，机械阻尼可被忽略，且有研究表明，叶片材料阻尼较小，对动力响应影响不大，于是气动阻尼就成为影响叶轮阻尼特性的支配因素。气动阻尼表征叶片与周围气体之间的能量传递，其值受诸多气流参数影响，目前有关气动阻尼的文献数据有限，且研究主要针对轴流叶片展开。Crawley[1]较早进行了该方面的工作，对一跨声速压气机展开了转子叶片气动阻尼的测量，叶片由上游进气扰动装置激振，测量得到了叶片不同叶间相角下的气动阻尼。Kielb和Abhari[2]通过实验测量了涡轮叶片的气动和结构阻尼，将在一定进气压力条件下测得的叶片总阻尼，减去空载条件下测得的阻尼，即得到叶片的气动阻尼。Kammerer和Abhari[3]针对某高速离心压缩机叶轮测量了其在不同工况条件下的阻尼数据，着重分析了进气压力和工况点对叶片模态阻尼的影响。Kammerer[4]开展了离心压缩机叶片受迫振动的实验研究，测量了叶轮的阻尼数据，实测结果表明，气动阻尼随着进气压力的提高而增加；当进气压力为0.1MPa时，气动阻尼值可达材料阻尼的10倍。
应用计算流体力学技术，亦有研究采用数值分析手段获得叶片的气动阻尼。Parthasarathy[5]利用数值模拟得到了某跨声速风扇叶片的气动阻尼，并研究了叶间相角对气动阻尼特性的影响。Srivastava等[6]用能量法计算了叶片气动阻尼，分析了出口反压对气动阻尼的影响，以及不同反压下激波位置变化对气动阻尼的影响。Seinturier等[7]分别用Euler方程和线化N-S方程计算了某压气机叶片二弯模态的气动阻尼，并将其用于气流激振的受迫响应计算。Elder等[8]采用傅式转换法求解周期性流场，只使用两个通道计算，得到的NASA Rotor 67叶片不同叶间相角时的气动阻尼与全通道计算结果符合很好。Im和Zha[9]采用气固耦合方法研究了一跨声速风扇转子在前行波振动下的颤振特性。根据计算得到的叶片位移响应，可求得振动的气动阻尼值。结果表明随工况点向失速边界移动，气动阻尼逐渐由正变负。郭雪莲等[10]通过数值计算分析了亚声速和跨声速转子叶片的模态气动阻尼比，结果表明这两类叶片的气动阻尼特性基本一致。杨晓东等[11]基于能量法建立了跨声速实验转子的全周气动阻尼计算模型，数值分析了转子叶片频率失谐对其气弹稳定性的影响。失谐的存在弱化了叶间相角对叶片气动阻尼的影响，失谐转子中不同叶片的气动阻尼表现出显著差异。
上述有关气动阻尼的计算研究均以轴流叶片为对象，目前少有针对离心叶轮气动阻尼展开计算分析的文献。由于转子叶盘的结构特点，其叶片往往呈现行波型振动，特别对于轴流叶片一类高而短的结构，可假设所有叶片均固支于其根部振动，各叶片以同样的频率和振幅做简谐振动，相邻叶片的振动只相差一个叶间相角，该方法被广泛应用于轴流叶片气动阻尼的计算中。但对于离心叶轮结构，因其叶片展向短而弦向长，叶-盘振动必须整体分析，不同节径的振动不仅叶间相角不同，频率亦发生变化，故离心叶轮气动阻尼的分析方法与轴流叶片有所不同。因气动阻尼的计算涉及网格变形、流固数据传递和运动边界流场分析等多项数值技术，已有计算研究多采用商业软件实现。而目前商业软件对于大规模网格的变形和壁面距离计算，缺乏高效的解决方案，因而限制了离心叶轮气动阻尼问题的计算研究。
本文为实现离心叶轮在行波振动下的气动阻尼计算，独立开发了网格变形程序和流动分析程序，分别利用紧支撑径向基函数和二叉树技术，大大提高了网格变形和壁面距离的计算效率。以半开式叶轮为研究对象，针对盘振动模态在前、后行波振动条件下引起的叶轮气动响应和气动功率密度展开了计算分析，并得到了叶轮的气动阻尼数据。
2 数值方法 2.1 控制方程考虑运动网格，将相对旋转坐标系下的可压缩流雷诺时均Navier-Stokes方程转换到任意曲线坐标系(ξ, η, ζ)下，有
| $ \frac{{\partial \mathit{\boldsymbol{Q}}}}{{\partial t}} + \frac{{\partial \left( {\mathit{\boldsymbol{F}} - {\mathit{\boldsymbol{F}}_{\rm{v}}}} \right)}}{{\partial \xi }} + \frac{{\partial \left( {\mathit{\boldsymbol{G}} - {\mathit{\boldsymbol{G}}_{\rm{v}}}} \right)}}{{\partial \eta }} + \frac{{\partial \left( {\mathit{\boldsymbol{H}} - {\mathit{\boldsymbol{H}}_{\rm{v}}}} \right)}}{{\partial \zeta }}\mathit{\boldsymbol{S}} $ | (1) |
式中Q为守恒变量，Fv，Gv，Hv为粘性通量，无粘通量F，G，H和源项S的表达式为
| $ \begin{array}{l} \mathit{\boldsymbol{F}} = \frac{1}{\mathit{\boldsymbol{J}}}\left( {\begin{array}{*{20}{l}} {\rho U}\\ {\rho Uu + {\xi _x}p}\\ {\rho Uv + {\xi _y}p}\\ {\rho Uw + {\xi _z}p}\\ {\left( {E + p} \right)U - {\xi _t}p} \end{array}} \right)\;\;\;\;\;\;\mathit{\boldsymbol{G}} = \frac{1}{\mathit{\boldsymbol{J}}}\left( {\begin{array}{*{20}{l}} {\rho V}\\ {\rho Vu + {\eta _x}p}\\ {\rho Vv + {\eta _y}p}\\ {\rho Vw + {\eta _z}p}\\ {\left( {E + p} \right)V - {\eta _t}p} \end{array}} \right)\\ \mathit{\boldsymbol{H}} = \frac{1}{\mathit{\boldsymbol{J}}}\left( {\begin{array}{*{20}{l}} {\rho W}\\ {\rho Wu + {\zeta _x}p}\\ {\rho Wv + {\zeta _y}p}\\ {\rho Ww + {\zeta _z}p}\\ {\left( {E + p} \right)W - {\zeta _t}p} \end{array}} \right)\;\;\;\;\;\;\mathit{\boldsymbol{S}} = \frac{1}{\mathit{\boldsymbol{J}}}\left[ {\begin{array}{*{20}{l}} 0\\ {\rho \omega v}\\ { - \rho \omega u}\\ 0\\ 0 \end{array}} \right] \end{array} $ | (2) |
式中J为雅可比矩阵，ω为绕z轴旋转的角速度，E为总能量，p为压力，ρ为密度，U，V，W为逆变速度，其表达式为
| $ \begin{array}{*{20}{l}} {U = {\xi _x}\left( {u + \omega y - {x_t}} \right) + {\xi _y}\left( {v - \omega x - {y_t}} \right) + {\xi _z}\left( {w - {z_t}} \right)}\\ {V = {\eta _x}\left( {u + \omega y - {x_t}} \right) + {\eta _y}\left( {v - \omega x - {y_t}} \right) + {\eta _z}\left( {w - {z_t}} \right)}\\ {W = {\zeta _x}\left( {u + \omega y - {x_t}} \right) + {\zeta _y}\left( {v - \omega x - {y_t}} \right) + {\zeta _z}\left( {w - {z_t}} \right)} \end{array} $ | (3) |
| $ \begin{array}{l} {\xi _t} = - {\xi _x}\left( { - \omega y + {x_t}} \right) - {\xi _y}\left( {\omega x + {y_t}} \right) - {\xi _z}{z_t}\\ {\eta _t} = - {\eta _x}\left( { - \omega y + {x_t}} \right) - {\eta _y}\left( {\omega x + {y_t}} \right) - {\eta _z}{z_t}\\ {\zeta _t} = - {\zeta _x}\left( { - \omega y + {x_t}} \right) - {\zeta _y}\left( {\omega x + {y_t}} \right) - {\zeta _z}{z_t} \end{array} $ | (4) |
式中u，v，w为绝对速度分量，xt，yt，zt为网格运动速度分量。
程序中无粘项离散可使用二阶中心格式和迎风格式，粘性项离散使用中心格式。考虑到目前PC机和工作站已普遍使用多核多线程处理器，程序引入OPENMP共享内存并行编译技术，通过向程序中添加少量伪代码实现了单机多线程并行计算，显著减少了时间成本。
2.2 多块网格的块间数据交换因非结构网格空间填充效率不高，为适应叶轮流道的复杂几何外形，程序基于多块结构化网格开发。对于多块网格，块间对接面为非物理边界，边界条件的施加通过块间数据交换实现。目前实现的是完全匹配对接，即点对点对接。参见图 1，在网格块间对接处，为构造空间离散格式需要边界面元两侧的单元变量值，此时可将对接网格块边界内侧的两层网格作为当前网格块的哑元网格，即当前网格块对接面外侧哑元网格的变量值可直接取对接块内侧网格处的值。
|
Fig. 1 Interblock data exchange |
对粘性项采用中心差分离散，二阶导数按网格面上一阶导数的差分处理。以i向网格面为例，粘性项按下式离散
| $ {\left( {{\delta _\xi }{F_{\rm{v}}}} \right)_{i, j, k}} = {\left( {{F_{\rm{v}}}} \right)_{i + 1/2, j, k}} - {\left( {{F_{\rm{v}}}} \right)_{i - 1/2, j, k}} $ | (5) |
对于有限体积方法，当多块结构化网格呈现图 2所示的拓扑结构时，采用高斯定理在偏移控制体积上积分计算网格面上的速度导数和温度梯度，块间几何数据交换将会出现歧义。一般叶片前缘的叶顶间隙网格就具有这种拓扑结构，此时可采用雅可比变换法来计算面元上的导数，同样以i向网格面为例
| $ {\left( {\frac{{\partial \phi }}{{\partial x}}} \right)_{i + 1/2, j, k}} = {\left( {\frac{{\partial \phi }}{{\partial \xi }}{\xi _x} + \frac{{\partial \phi }}{{\partial \eta }}{\eta _x} + \frac{{\partial \phi }}{{\partial \zeta }}{\zeta _x}} \right)_{i + 1/2, j, k}} $ | (6) |
| $ \begin{array}{*{20}{l}} {{{\left( {\frac{{\partial \phi }}{{\partial \xi }}} \right)}_{i + 1/2, j, k}} = {\phi _{i + 1, j, k}} - {\phi _{i, j, k}}}\\ {{{\left( {\frac{{\partial \phi }}{{\partial \eta }}} \right)}_{i + 1/2, j, k}} = {\phi _{i + 1/2, j + 1/2, k}} - {\phi _{i + 1/2, j - 1/2, k}}}\\ {{{\left( {\frac{{\partial \phi }}{{\partial \zeta }}} \right)}_{i + 1/2, j, k}} = {\phi _{i + 1/2, j, k + 1/2}} - {\phi _{i + 1/2, j, k - 1/2}}} \end{array} $ | (7) |
| $ \begin{array}{*{20}{l}} {{{\left. {{\xi _x}} \right|}_{{\rm{}}i + 1/2, j, k}} = {{\left. {J\left( {{y_\eta }{z_\zeta } - {y_\zeta }{z_\eta }} \right)} \right|}_{{\rm{}}i + 1/2, j, k}}}\\ {{{\left. {{\eta _x}} \right|}_{{\rm{}}i + 1/2, j, k}} = {{\left. {J\left( {{y_\zeta }{z_\xi } - {y_\xi }{z_\zeta }} \right)} \right|}_{{\rm{}}i + 1/2, j, k}}}\\ {{{\left. {{\zeta _x}} \right|}_{{\rm{}}i + 1/2, j, k}} = {{\left. {J\left( {{y_\xi }{z_\eta } - {y_\eta }{z_\xi }} \right)} \right|}_{{\rm{}}i + 1/2, j, k}}} \end{array} $ | (8) |
|
Fig. 2 Ambiguous topology for interblock data exchange |
式中ϕ表示速度、温度等物理量，ϕ在半点处的值可以用周围格心处的值插值得到，如
| $ {\phi _{i + 1/2, j + 1/2, k}} = \left( {{\phi _{i, j, k}} + {\phi _{i + 1, j, k}} + {\phi _{i, j + 1, k}} + {\phi _{i + 1, j + 1, k}}} \right)/4 $ | (9) |
使用双时间方法[12]实现非定常流场求解，即在每一物理时间步引入一个虚时间迭代过程，气动方程在虚时间域内推进到定常解时即得到每一物理时间步的流场解，该迭代由下式表达
| $ {(\frac{{\partial \mathit{\boldsymbol{Q}}}}{{\partial \tau }})^{m + 1}} = - \frac{{3({\mathit{\boldsymbol{Q}}^{m + 1}} - {\mathit{\boldsymbol{Q}}^n}) - ({\mathit{\boldsymbol{Q}}^n} - {\mathit{\boldsymbol{Q}}^{n - 1}})}}{{2\Delta t}} + \mathit{\boldsymbol{R}}({\mathit{\boldsymbol{Q}}^{m + 1}}) $ | (10) |
式中τ和t分别为虚时间和物理时间，上标m和n分别表示虚时间步和物理时间步。虚时间域内的时间推进采用LU-SGS隐式算法。
对于运动边界的有限体积法流场计算，为避免控制体积变形引入守恒变量的误差，必须满足几何守恒律[13]
| $ {\left( {1/J} \right)_t} - \sum {\left\{ {{x_t},{y_t},{z_t}} \right\}} \cdot \mathit{\boldsymbol{n}}A = 0 $ | (11) |
式中1/J等于网格单元体积，∑表示对网格面求和。几何守恒率的物理意义表示单元变形引起的体积增量为运动网格面扫过的净体积。
2.5 气动网格变形气动阻尼由叶轮振动引起，而叶轮振动的计算使用的是有限元网格。为获得气动阻尼，就要求解气动网格运动变形下的流场，因而必须将结构表面变形数据传递至气动表面。本文研究关注叶轮模态振动引起的气动阻尼，为此需先将某阶模态的结构表面位移插值到气动表面，进而展开气动网格的变形。结构/气动数据传递和网格变形一般被视作不同性质的问题，然而借助径向基函数工具，两类问题可容易地按相同方法处理。
| $ s\left( \mathit{\boldsymbol{x}} \right)=\mathop {\mathop \sum \limits^N }\limits_{i = 1}\, {{\alpha }_{i}}\varphi \left( \left\| \mathit{\boldsymbol{x}}-{{\mathit{\boldsymbol{x}}}_{i}} \right\| \right)+p\left( \mathit{\boldsymbol{x}} \right) $ | (12) |
式中s (x)是在点x处被估计的函数，这里指变形量，φ是径向基函数，xi表示已知结构表面点。p (x)是多项式函数，可采用如下线性多项式
| $ p\left( \mathit{\boldsymbol{x}} \right) = {\gamma _0} + {\gamma _x}x + {\gamma _y}y + {\gamma _z}z $ | (13) |
对N个结构点列位移方程，并补充附加条件
| $ \begin{matrix} s\left( {{\mathit{\boldsymbol{x}}}_{j}} \right)=\mathop {\mathop \sum \limits^N }\limits_{i = 1}\, {{\alpha }_{i}}\varphi \left( \left\| {{{\mathit{\boldsymbol{x}}}_{j}}-{{\mathit{\boldsymbol{x}}}_{i}} } \right\| \right)+p\left( {{\mathit{\boldsymbol{x}}}_{j}} \right), \text{ }\!\!~\!\!\text{ }j=1, 2, \cdots , N \\ \mathop {\mathop \sum \limits^N }\limits_{i = 1}\, {{\alpha }_{i}}q\left( {{\mathit{\boldsymbol{x}}}_{i}} \right)=0 \\ \end{matrix} $ | (14) |
即可构造N+4个线代方程，求出插值系数αi，γ0，γx，γy和γz，式中q (x)是阶数小于等于p (x)的多项式。
然而径向基函数方法同样存在计算量和内存消耗较大的问题，为使求解矩阵呈稀疏和带状分布，以利于求解大型数据插值问题，Wendland[16]定义了如下的紧支撑径向基函数
| $ \phi \left( r \right) = \left\{ {\begin{array}{*{20}{l}} {{{\left( {1 - r} \right)}^k}g\left( r \right), \;\;\;\;\;\;\;\;\;\;0 \le r \le 1}\\ {0, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; r > 1} \end{array}} \right. $ | (15) |
式中k为正整数g (r)为关于r的多项式，
对于叶轮的气动网格变形，除叶轮表面外，还应考虑机匣面的限制。用紧支基函数将结构模态变形传递至气动网格后，采用下式进行网格变形的衰减控制
| $ {\rm{d}}\mathit{\boldsymbol{x'}} = {{\rm{e}}^{ - \sigma {d_{\rm{r}}}/{d_{\rm{s}}}}}{\rm{d}}\mathit{\boldsymbol{x}} $ | (16) |
式中dr，ds分别表示网格点到叶轮表面和机匣面的距离，σ为控制参数。
2.6 壁面距离的快速算法对于运动边界的非定常流动分析，在每一时间步做网格变形后，需重新计算各单元格心到壁面的最小距离，以供湍流计算使用。目前对壁面距离常用直接计算方法，即对于各流场网格点，计算其到每个壁面网格点的距离，通过逐点比较，取最小值为该流场网格点的壁面距离。对大规模网格，直接计算方法的计算量很大，会大大降低非定常流动的计算效率。本文采用一种二叉树搜索技术[17, 18]，实现了一种快速的壁面距离算法。
以二维空间为例，具体步骤为：对第一个壁面网格点建立一个树节点A，使该节点指向一块包含流场空间的长方区域R，并将此区域沿x坐标等分为二，一个子区域Rl为R左半部，另一子区域Rr为R右半部。对第二个壁面网格点，搜索包含该点的子区域Rl或Rr，建立子节点B指向该子区域，并将此区域再沿y坐标等分，得到两个子区域。对后续壁面网格点，遍历当前二叉树，搜索包含该点的最小子区域，建立子节点指向该子区域，再依次沿x和y坐标做区域分割，以此类推，直至建立包括所有壁面网格点的二叉树，其过程如图 3所示。
|
Fig. 3 Binomial tree model of geometric segmentation |
计算一流场网格点P1到壁面网格点的最小距离d1，任意流场网格点P2到壁面网格点的最小距离d2应满足d2 ≤ d1 + d1-2，其中d1-2为P1~P2的距离。令r = d1 + d1-2，则距点P2最近的壁面点必然位于以点P2为圆心，以r为半径的圆球内。为计算方便，可将搜索区域由圆球转换为其外切长方盒Q。
确定搜索范围后，遍历前面建立的二叉树，实施步骤为：若当前节点位于待搜索的长方盒Q内，则计算P2到该节点对应壁面网格点的距离；若左子节点非空，且方盒Q与当前节点的左子区域Rl相交，则遍历左子树；若右子节点非空，且方盒Q与当前节点的右子区域Rr相交，则遍历右子树。
3 算例验证 3.1 径向叶轮流场该算例为一径向叶轮的非定常流场，叶轮转速22790r/min，叶片数16，轮缘直径175mm，出口宽度10.5mm。气动网格由4块H型网格搭接而成，格点数为550420，如图 4所示。计算设置进口总压101.3kPa，总温293.15K。程序和软件计算均采用二阶格式和S-A湍流模型。计算设置了两个压力监控点，位于图 4所标示轴向位置的流道中心处。
|
Fig. 4 Aerodynamic grid of the radial impeller |
出口背压为135kPa时，计算得到的无叶扩压段出口平均周向速度的轴向分布如图 5所示，横轴表示轴向坐标。程序计算的周向速度分布与商业软件的结果符合很好，在两侧近壁面区域非常一致，显示程序的正确性。
|
Fig. 5 Comparison of average circumferential velocity at the diffuser outlet |
考虑对叶轮上游进口总压施加脉冲，下游扩压段出口背压设为130kPa，计算叶轮对压力脉冲的响应流场。压力脉冲呈现为一突降尖峰，变化幅度为初始总压的5%，脉冲宽度为1ms。定义无量纲时间t*为物理时间除以脉冲宽度，脉冲在t*=0.1时出现，t*=0.6时达到峰值，t*=1.1时脉冲结束。计算时间步长设为0.1ms。
图 6所示分别为叶轮的监控点压力和气动推力响应曲线。因未对轮盘背面间隙建模，气动推力只对叶盘的迎风面计算。由图可见程序计算的响应曲线与商业软件结果一致性较好，其中监控点1的压力曲线几乎重合。两者计算的响应曲线变化趋势完全相同，各峰值时刻完全对应。比较结果显示了程序模拟叶轮流场的正确性。
|
Fig. 6 Unsteady flow responses of the impeller to intake pressure pulse |
叶栅采用NACA 0012翼型，进口马赫数为0.593，雷诺数为3.21×106，进口气流角为0°，栅距/弦长等于1。叶片绕前缘做扭转振动，扭转振幅为2°，相邻叶片振动相位差为零，振动频率为64Hz。进、出口段均采用H型网格，叶片段采用C型网格，各块网格数分别为49×13×9，57×37×9和129×29×9。计算采用中心格式进行空间离散。
计算设定叶片表面网格做刚体随动，网格变形量按离开叶表距离逐渐衰减，直至进、出口边界和周期边界网格变形量衰减为零。计算得到的叶片一阶非定常压力系数如图 7所示，横轴表示相对弦长，显见本文计算结果与Huff[19]的结果符合较好，实部曲线几乎一致。比较结果验证了程序求解运动边界流场的正确性。
|
Fig. 7 1st unsteady pressure coefficients of the vibrating cascade |
考虑振动叶栅为单自由度气弹系统的情况，气动阻尼比定义为
| $ {\zeta _{{\rm{aero}}}} = \frac{{{W_{{\rm{damping}}}}}}{{2{\rm{ \mathit{ π} }}I\omega _0^2\theta _{{\rm{max}}}^2}} $ | (17) |
式中ω0为固有角频率，θmax为扭转振幅，I为转动惯量。Wdamping为一个周期内的气动阻尼功，由下式计算
| $ {W_{{\rm{damping}}}} = - \mathop \int_0^T Q\left( t \right){\rm{d}}t $ | (18) |
式中Q (t)为叶表气动功率。
气动边界条件同前设置，计算得到了不同频率下的气动阻尼功和气动阻尼比，如图 8所示。显见随着振动频率的增加，一个周期的气动阻尼功呈线性增长态势，但气动阻尼比不断下降。由图中数据，气动阻尼功与振幅的平方成正比，但不同振幅下的气动阻尼比几乎不变。文献[20]通过实验测量了一线性涡轮叶栅的气动阻尼，叶片绕其弾性轴振动，实验表明气动阻尼比与振幅无关。本例计算结果与此相符，验证了程序计算气动阻尼问题的正确性。
|
Fig. 8 Aero damping characteristics of the vibrating cascade |
考虑到离心压缩机实际工作状态复杂，硬件建设、实验设计和设备调试的成本较高，本文采用数值分析的手段，进行振动叶轮气动阻尼的计算研究工作。
研究对象为一单级离心压缩机的半开式径向叶轮，其出口直径280mm，出口宽度16.8mm，叶顶间隙0.8mm，叶片数16。为方便处理网格变形，本文计算对叶轮流场整体建模，气动网格由80块H型网格拼接而成，网格点数为8722384，最小网格角为8.6°，如图 9所示。显见数值模型考虑了轮盘与机匣间隙中泄漏流动的影响，常规分析在计算叶轮流场时只对轮盘迎风面的流道划分网格，对轮盘背面的泄漏流道不予建模，这对叶轮气动阻尼的计算有显著影响。使用上述开发的程序计算，对流项采用中心格式离散，湍流模型采用S-A模型。设定进口总压101.3kPa，总温293.15K，转速14250r/min，叶盘背面泄漏出口给定质量流量3g/s。非定常流场计算由定常流场初始化，虚时间域的推进步数设为30。
|
Fig. 9 Flow grid for the calculation of aerodynamic damping |
为获得叶轮在行波振动情况下的流场数据，流场分析前，首先应用有限元软件进行了叶轮的模态分析。结构网格如图 10所示，采用10节点四面体单元，节点数为173968。弹性模量取72GPa，泊松比0.33，密度2700kg/m3。图 11为计算得到的第23模态的质量矩阵归一化振型，振型位移为无量纲量。该模态为叶轮第3节径的第3模态，以轴向变形为主，是典型的由盘振动支配的模态。盘振动是离心叶轮的重要振动形式，这与轴流转子往往只考虑叶片振动不同。由振型分布可以看出该模态显著呈现出3条节径线，在行波振动条件下，相邻两扇区的振动相角差均为3π/8。对于离心叶轮，不同节径的振动模态具有不同的频率和周向波数，叶片和轮盘振动耦合较强，叶轮振动计算必须做整体分析，这与轴流叶片的情况有显著差异。
|
Fig. 10 Finite element mesh of the impeller |
|
Fig. 11 Modal shape of the impeller, Ω=3184Hz |
为实现上述盘振动模态在行波振动下的运动边界流场计算，需将模态变形由结构网格传递至气动网格。对于径向基函数法，要求解式（12）的插值系数，需要构造173972×173972的稠密系数矩阵，其内存需求PC机无法满足。采用紧支撑径向基函数法后，系数矩阵非零元素个数缩减为17051442。图 12为气动表面网格表达的第23模态振型，与图 11中的振型位移相比高度一致，插值精度得到了良好保证。
|
Fig. 12 Modal shape expressed by the aerodynamic grid |
应用前述的紧支径向基函数方法完成气动网格的变形，以上述模态为例，设叶轮的最大变形量为0.2mm，叶轮叶片出口区域的变形网格如图 13所示。红色实体为未变形叶轮，红色线条为叶顶间隙网格未变形时的边缘曲线，绿色透明实体为变形后叶轮，绿色线条为叶顶间隙网格变形后的边缘曲线。在由机匣到叶顶的间隙内，网格由零位移到较大变形，发生一定的扭曲，但变形后网格的最小网格角为8.5°，与原网格的8.6°相差很小，这表明网格变形算法的保角性较好。应用式（12）对每个气动网格点计算位移量时，都要计算该点到各结构表面点xi的欧氏距离
|
Fig. 13 Deformed aerodynamic grid |
对于运动边界的流场模拟，在每一物理时间步上要逐次进行网格变形、壁面距离计算和流动方程求解。完成网格数8722384规模的壁面距离计算，采用直接算法将会十分耗时，在双核四线程的单机上实测需17min 8s，而如果采用前述二叉树算法，则只需36s。每一物理时间步的流动方程求解实测需12min 45s，可见壁面距离的快速计算节省了50%以上的时间成本，这对本文实现叶轮气动阻尼的计算有重要帮助。
因叶片前缘为圆头，划分的叶顶间隙网格在叶片前缘具有图 2所示的拓扑结构，若使用偏移控制体积计算网格面元上的速度导数，块间几何数据交换将出现歧义，进而直接影响流场计算的收敛性。采用2.3节所述的算法较好地避免了该问题，使流场计算表现出良好的收敛性。
计算针对上述叶轮模态、背压135kPa的设计工况展开。设置一个振动周期52个物理时间步，叶轮最大变形量0.2mm，数值模拟了叶轮在前、后行波振动条件下的气动响应和阻尼特性。实际中离心叶轮可能受各种激振源作用而引起行波振动，在动坐标系中观察，前行波顺盘传播，后行波逆盘传播。图 14为叶轮第23模态在后行波振动下各瞬时的轴向变形情况，该变形是振型位移的轴向分量，为无量纲量。显见行波振动传播一周，叶轮附近的流场边界周期变化3次，行波振动将导致叶轮流场的非定常响应。
|
Fig. 14 Time marching of backward travelling wave vibration |
叶轮起振后，流场由初始定常状态经历一段瞬态响应并逐渐趋于准定常变化。计算在不同叶片流道选取了若干压力监控点，各点的压力响应均经历了相似的变化过程，最终呈现周期变化的特征。图 15所示为其中一监控点的压力响应曲线，非定常计算起始后经过一段时间的变化，压力振荡的幅度逐渐趋于稳定值，表明叶轮流场的响应开始进入周期变化的准定常状态。由图可见，前、后行波振动引起的压力响应不同，监控点在后行波振动下的压力变化幅值较小。
|
Fig. 15 Pressure curves of monitor point |
前、后行波振动导致了不同的流场响应，图 16（a）显示了叶轮表面的气动推力在计算过程中的收敛曲线，前行波和后行波振动引起的气动推力变化过程不同，并最终导致了不同的气动推力，后行波振动引起的气动推力更大。将叶表各面元的气动力矢量点乘振动速度，并对叶轮表面积分即可得到作用于叶轮的气动功率。图 16（b）所示为叶轮振动引起的气动功率在计算过程中的收敛曲线，显见后行波振动引起的气动功率较大，反映后行波振动具有较大的气动阻尼，气弹稳定性更好。
|
Fig. 16 Convergence curves of aerodynamic responses |
一个周期内叶轮的气动阻尼功同式（18）的定义。振型按质量矩阵正交归一化后，气动阻尼比由下式计算
| $ {\zeta _{{\rm{aero}}}} = \frac{{{W_{{\rm{damping}}}}}}{{2{\rm{ \mathit{ π} }}{\Omega ^2}A_{{\rm{max}}}^2}} $ | (19) |
式中Ω表示模态圆频率，Amax表示相对于模态的广义振幅。在本文研究的振动工况下，气动力对叶轮整体做负功，气动阻尼为正。前行波振动时一个振动周期的气动阻尼功为0.0832J，气动阻尼比为4.08×10-3；后行波振动时一个振动周期的气动阻尼功为0.133J，气动阻尼比为6.53×10-3。
图 17（a）为前行波振动下某时刻叶轮表面的气动功率密度分布，图中显示的叶轮表面为迎风面，轮盘背面的气动功率密度符号相反、分布相似。因该模态为轮盘振动支配，轮盘背面泄漏流道的气体必然对振动叶轮做功，进而影响叶轮的气动阻尼特性。由图可见，最大气动功率密度的位置均位于轮盘边缘处，若不考虑气动功率的正负，其分布云图将与图 10所示的叶轮模态振型比较近似，说明叶轮表面的气动功率密度受当地振动速度的影响较大。
|
Fig. 17 Aerodynamic power density on the impeller surface under the travelling wave vibrations |
图 17（b）显示了后行波振动下同一时刻叶轮表面的气动功率密度分布，显见该云图同样大致体现振型分布的特点，但气动功率密度的正负极值较前行波振动的情况更高，说明振动引起的气动阻尼更大。因该模态行波振动下轮盘的振动主要为轴向，这使叶片出口部分的压力面和吸力面气动功率密度近于零，气动力做功主要体现为气动力抵抗叶盘轴向变形的能力。图 17（a）和（b）所示均为叶轮在振动周期同一时刻的云图，该时刻前、后行波在叶轮周向处于同一位置，比较两者可见因振动传播方向相反，导致同一位置的气动功率密度符号亦相反。
行波振动条件下叶表的气动功率密度分布在任意时刻均体现出相应振型的特点，但对气动功率密度在一个振动周期上进行时间积分后，其云图则呈现循环对称特征。图 18为叶轮在后行波振动情况下一个振动周期内的平均气动功率密度分布，显见每一扇区气动功率密度的分布情况相同。平均气动功率反映一个周期内气动积累功的大小，体现气动阻尼的性质。由图可知，叶轮表面大部分区域的平均气动功率密度为负，轮盘背面区域的平均气动功率密度全为负，说明泄漏流动对叶轮振动起正阻尼作用。前行波振动下的平均气动功率密度绝对值略低，分布与该图相似，这里不再赘述。
|
Fig. 18 Average aerodynamic power density under the backward travelling wave vibration |
采用紧支撑径向基函数和二叉树技术，提升网格变形和壁面距离的计算效率，开发复杂运动边界的流场模拟程序，实现了振动叶轮气动响应数据的高效计算工作。对叶轮在行波振动下的气动阻尼计算表明：
（1）采用紧支径向基函数法进行结构-气动变形数据插值大大减小了内存需求，并具有足够的精度。应用二叉树搜索技术后，紧支基函数法的时间成本亦大大降低。实测结果表明二叉树算法使气动网格变形和湍流计算所需壁面距离的计算效率提高了至少一个数量级。
（2）考察的盘振动模态在行波振动情况下计算的气动阻尼比为正，气动力对叶轮整体做负功，阻碍叶轮的振动。前行波和后行波振动导致不同的流场响应和气动阻尼，后行波振动引起的气动推力和气动阻尼较大，气弹稳定性更好。
（3）行波振动条件下任意时刻的叶轮表面气动功率密度主要受当地振动速度的影响，其分布呈现类似振型的特点。叶轮一个振动周期的平均气动功率密度则呈现循环对称的特征，叶表大部分区域的平均气动功率密度为负。
除行波振动的传播方向外，离心叶轮的气动阻尼还可能受到诸多气流和结构参数的影响，如气动工况和结构约束等，其力学机制有待进一步的深入研究。此外，提高气动阻尼对减小叶轮受迫振动、进而改善疲劳寿命的影响，有待将来实验工作的开展和试验数据的支撑。
| [1] |
Crawley E F. Aerodynamic Damping Measurements in a Transonic Compressor[J]. ASME Journal of Engineering for Power, 1983, 105(3): 575-584. DOI:10.1115/1.3227456
(  0) |
| [2] |
Kielb J J, Abhari R S. Experimental Study of Aerodynamic and Structural Damping in a Full-Scale Rotating Turbine[J]. ASME Journal of Engineering for Gas Turbines and Power, 2003, 125(1): 102-112. DOI:10.1115/1.1496776
( 0) |
| [3] |
Kammerer A, Abhari R S. Experimental Study on Impeller Blade Vibration During Resonance, Part 2: Blade Damping[J]. ASME Journal of Engineering for Gas Turbines and Power, 2009, 131(2).
( 0) |
| [4] |
Kammerer A. Experimental Research into Resonant Vibration of Centrifugal Compressor Blades[D]. Zurich: Swiss Federal Institute of Technology, 2009. http://dx.doi.org/10.3929/ethz-a-006026548
( 0) |
| [5] |
Parthasarathy V. Computation of Aerodynamic Damping for Flutter Analysis of a Transonic Fan[R]. ASME GT 2011-46597. http://www.researchgate.net/publication/225007516_Computation_of_Aerodynamic_Damping_for_Flutter_Analysis_of_a_Transonic_Fan
( 0) |
| [6] |
Srivastava R, Bakhle M A, Keith T G. Numerical Simulation of Aerodynamic Damping for Flutter Analysis of Turbomachinery Bladed Rows[J]. Journal of Propulsion and Power, 2003, 19(2): 260-267. DOI:10.2514/2.6107
( 0) |
| [7] |
Seinturier E, Lombard J P, Dumas M, et al. Forced Response Prediction: Methodology for the Design of HP Compressors Bladed Disks[R]. ASME GT 2004-53372. http://www.researchgate.net/publication/290493077_Forced_response_prediction_methodology_for_the_design_of_HP_compressors_bladed_disks
( 0) |
| [8] |
Robin Elder, Ian Woods, Sunil Patil, et al. Investigation of Efficient CFD Methods for the Prediction of Blade Damping[R]. ASME GT 2013-95005. http://www.researchgate.net/publication/267504310_Investigation_of_Efficient_CFD_Methods_for_the_Prediction_of_Blade_Damping
( 0) |
| [9] |
Im H S, Zha G C. Flutter Prediction of a Transonic Fan with Travelling Wave Using Fully Coupled Fluid-Structure Interaction[R]. ASME GT 2013-94341. http://www.researchgate.net/publication/267650393_Flutter_Prediction_of_a_Transonic_Fan_With_Travelling_Wave_Using_Fully_Coupled_FluidStructure_Interaction
( 0) |
| [10] |
郭雪莲, 李琳. 亚声速叶片与跨声速叶片的气动阻尼比较[J]. 推进技术, 2010, 31(3): 296-308. (GUO Xue-lian, LI Lin. Comparison on Aerodynamic Damping between Subsonic Blade and Transonic Blade[J]. Journal of Propulsion Technology, 2010, 31(3): 296-308.)
( 0) |
| [11] |
杨晓东, 侯安平, 李漫露, 等. 频率失谐对跨声速压气机气弹稳定性的影响[J]. 推进技术, 2015, 36(4): 587-594. (YANG Xiao-dong, HOU An-ping, LI Man-lu, et al. Effects of Frequency Mistuning on Aeroelastic Stability of Transonic Compressor[J]. Journal of Propulsion Technology, 2015, 36(4): 587-594.)
( 0) |
| [12] |
Jameson A. Time Dependent Calculations Using Multigrid with Applications to Unsteady Flows past Airfoils and Wings[R]. AIAA 91-1596. http://arc.aiaa.org/doi/abs/10.2514/6.1991-1596
( 0) |
| [13] |
Sitaraman J, Baeder J D. Field Velocity Approach and Geometric Conservation Law for Unsteady Flow Simulations[J]. AIAA Journal, 2006, 44(9): 2084-2094. DOI:10.2514/1.5836
( 0) |
| [14] |
Buhmann M. Radial Basis Functions[M]. Cambridge: Cambridge University Press, 2003.
( 0) |
| [15] |
Allen C B, Rendall T C S. Unified Approach to CFD-CSD Interpolation and Mesh Motion Using Radial Basis Functions[R]. AIAA 2007-3804. https://arc.aiaa.org/doi/abs/10.2514/6.2007-3804
( 0) |
| [16] |
Wendland H. Piecewise Polynomial, Positive Definite and Compactly Supported Radial Basis Functions of Minimal Degree[J]. Advances in Computational Mathematics, 1995, 4(1): 389-396. DOI:10.1007/BF02123482
( 0) |
| [17] |
Bonet J, Peraire J. An Alternating Digital Tree (ADT) Algorithm for 3D Geometric Searching and Intersection Problems[J]. International Journal for Numerical Methods in Engineering, 1991, 31(1): 1-17. DOI:10.1002/(ISSN)1097-0207
( 0) |
| [18] |
Boger D A. Efficient Method for Calculating Wall Proximity[J]. AIAA Journal, 2001, 39(12): 2404-2406. DOI:10.2514/2.1251
( 0) |
| [19] |
Huff D L. Numerical Simulations of Unsteady, Viscous, Transonic Flow over Isolated and Cascaded Airfoils Using a Deforming Grid[C]. Hawaii: AIAA 19th Fluid Dynamics, Plasma Dynamics, and Lasers Conference, 1987.
( 0) |
| [20] |
Seeley C E, Wakelam C, Zhang X, et al. Investigations of Flutter and Aerodynamic Damping of a Turbine Blade: Experimental Characterization[J]. Journal of Turbomachiney, 2017, 139(8).
( 0) |