STUDY ON ENERGY BAND STRUCTURES AND ITS DYNAMICS CHARACTERISTICS IN CYCLIC-PERIODIC SYMMETRIC STRUCTURES
-
摘要: 针对叶轮这一类循环周期对称结构在服役中常出现的局部损坏及振动局部化问题, 在链式周期结构的能带理论基础上, 通过建立环向周期结构原胞中弹性波传播问题的计算方法, 得到了整体循环对称结构的能带结构及其独特动力学特性, 并揭示了振动局部化的形成机制. 首先, 采用Euler-Bernoulli梁模型和Kirchhoff板模型对离心叶轮较薄的区域进行力学建模, 同时根据对应的力学假设和边界条件对弯曲波方程进行求解; 进一步, 推导出了弯曲波沿周向波导传播的传递矩阵并计算了有限空间Lyapunov指数, 结果表明: 经典Euler-Bernoulli梁模型无法准确描述弯曲波沿周向波导的传播规律, 而Kirchhoff板模型可以在较宽的范围内给出5000 Hz附近的带隙; 最后, 采用有限元方法对离心叶轮结构的能带结构进行了数值计算, 并对计算结果进行了分析研究. 结果表明: 该离心叶轮中存在能带结构, 并具有3个“驻波型”带隙和1个“局域共振”带隙. 其中, “局域共振”带隙的存在是离心叶轮中出现振动局部化现象的原因之一, 即当外界作用在离心叶轮上的激振力频率落入“局域共振”带隙范围内时, 弯曲波在周向的传递受阻, 能量聚集于部分叶片上, 表现为“周期性循环对称破缺”. 研究结果揭示了循环对称结构动力学的独特特性, 能够为进一步完善离心叶轮的减振及气动/结构协同设计准则提供理论基础.Abstract: In view of the problems of local damage and vibration localization in cyclic-periodic symmetric structures such as impellers in service, the energy band structure and dynamic characteristics of the whole cyclic periodic symmetric structure are obtained by analyzing the dynamic characteristics of elastic wave propagation in the ring cell based on the band theory of chain periodic structures, meanwhile, formation mechanism of vibration localization is revealed. First, Euler-Bernoulli beam model and Kirchhoff plate model are used to model the thin region of the impeller, and the bending wave equation is solved with the corresponding mechanical assumptions and boundary conditions. Further, the transfer matrix of the bending wave propagating around the waveguide is derived in detail and the finite space Lyapunov indexes are given. The results show that the Euler-Bernoulli beam model cannot accurately describe the propagation law of curved wave along the circumferential waveguide, and the Kirchhoff plate model could give the band gap around 5000 Hz in a wide range. Finally, the finite element method is used to compute the energy band structure of the centrifugal impeller, and the results are analyzed and studied in detail. The results show that there is one energy band structure in the centrifugal impeller with 3 "standing wave" band gaps and 1 local resonance band gap. Among them, the existence of "local resonance" band gap is one of the reasons for vibration localization in the centrifugal impeller, that is, as the frequency of external excitation acting on the centrifugal impeller falls within the "local resonance" band gap, the circumferential transfer of bending wave is blocked, and the energy is collected on part of the blades, resulting in vibration localization, demonstrating “cyclic-periodic symmetry-breaking”. In summary, the results reveal the dynamic characteristics of cyclic symmetric structures, which can provide theoretical basis for further improving the design criteria and vibration suppression of centrifugal impellers in the future.
-
引言
近年来随着能源与动力领域透平机械向高参数(高压比、高转速和高效率等)趋势的发展, 这些变化在提高叶轮气动性能的同时也带来了许多新的结构方面的问题, 即: 在服役过程中频繁出现叶片的疲劳断裂损坏, 而现场振动症状、断裂位置及断口显示其破坏类型为振动局部化导致的结构疲劳断裂. 振动局部化现象是指弹性波在结构中传播时, 振动能量会聚集于局部区域而非均匀地分布在整个结构内, 使得局部区域的振动参与度明显高于其他部位. 对于该类问题, 以往的研究常从疲劳损坏、热蠕变、材料属性等方面对事故叶轮进行分析和改进[1-2], 但效果甚微.
事实上, 叶轮是一类典型的周期性循环对称结构, 其具有不同于非周期结构的许多独特的动力学性质, 如: 存在频率通带和禁带现象. 因此, 需要从理论上对该类循环周期结构振动局部化现象的形成机制进行研究, 该振动局部化表现为“周期性循环对称破缺”, 其中面临的首要问题是分析的理论方法.
在过去的研究中发现, 与链式周期波导结构类似, 叶轮这类循环周期对称结构也具有通带和禁带[3]. 而当激振力频率处于叶轮的禁带范围内时, 其振动响应呈现局部化、非对称性的特点. 这与固体物理领域声子晶体能带理论及晶格动力学理论中所描述的禁带响应特性, 尤其是“局域共振”带隙下的响应特性存在相似之处, 因此自然地联想到振动局部化与禁带响应特性之间应当存在关联. 另一方面, 由于叶轮在周向具有周期性, 同样属于周期结构, 只是周期延拓方向以极坐标系中的周向基矢量来描述. 且从数学角度, 该问题本质是链式周期结构中的定理和研究方法在坐标变换下的推广. 因此理论上该能带结构拓展并应用于叶轮是可行的. 下面首先对声子晶体能带理论的相关研究进行简要介绍.
声子晶体、超材料和超表面等结构在能带结构与振动传输特性方面的探索日益受到关注[4-13]. 声子晶体是一种具有周期性结构的人工复合材料[14-18], 主要可以分为两种类型: Bragg散射型[14]和局域共振型[19-20], 其中局域共振型通常适用于需要较高设计自由度的应用[21-24]. 2015年, 董立强[25]研究了广义声子晶体的弹性波带隙特性, 从理论、仿真和实验方面发现其能有效抑制特定频带内柱面波传播, 表明局部化因子分析方法对广义声子晶体的带隙研究是完全可行的. 由于叶轮的周期性延拓在周向, 属于θ类广义周期结构, 因此, 该类周期性研究方法与圆柱壳的周向研究方法类似; 2019年, 罗金雨等[26]利用有限元法计算了一种半径为0.1 m的圆柱壳声子晶体结构的能带结构, 发现其能带结构与板类结构类似, 具有良好的沿轴线和圆柱环两个方向的带隙或全禁带隙, 其振动传递损失也验证了该结论; 2020年, 黄洪赛等[27]提出一种通过在圆柱壳圆周方向布置弹簧振子实现的局域共振型圆柱壳类声子晶体结构, 由于圆柱壳和弹簧振子振动的耦合, 该结构能够形成两条低频带隙; 2023年, 姚凌云等[28]将周期分级排列的组合方式应用于圆柱壳类弹性波超材料的带隙拓宽中, 并利用有限元法进行能带结构和振动传输特性计算, 发现该方法可实现弯曲波带隙的拓宽.
基于以上研究可以发现: 能带理论在周向波导的能带结构计算中已经存在部分应用, 但主要集中于壳结构上, 而在离心叶轮类带叶片(子结构)的循环周期对称结构上尚未应用. 因此, 本文首先分别采用Euler-Bernoulli梁假设及Kirchhoff板假设下的传递矩阵法对弯曲波在离心叶轮类循环对称结构中的周向传播规律进行了研究; 之后采用有限元方法数值计算了离心叶轮结构的能带结构, 并对其中的带隙进行了详细分析. 以期从能量及弹性波的输运和迁移角度对振动局部化损坏现象进行解释, 同时一定程度上对叶轮类循环周期对称结构的减振控制规范的完善提供理论基础.
1. 循环对称结构的动力学计算模型
1.1 周向原胞的划分
为分析问题方便, 以典型的、具有普适性的循环对称半开式离心叶轮为例开展相关研究, 其具有周向循环对称、弱耦合和薄壁结构的特点, 由叶片及轮盘构成, 如图1所示. 根据链式周期结构能带结构的研究理论, 在计算循环周期结构的色散曲线前, 首先需确定单胞及其相关参数.
研究对象为含18扭曲叶片且直径为250 mm的离心叶轮, 叶轮边缘的叶-盘结构与链式结构中基梁-立柱结构类似, 故单胞的选择应包含叶片及轮盘. 理论上单胞内叶片选取数量可以包含18的任意整因数, 各单胞组内叶片数量应相等, 但考虑到由于圆形结构不存在链式结构所具有的明确起点与终点, 即弯曲波激励点无论选在轮盘边缘上任意空间位置均能够产生双向行波, 所以周期结构单胞应小于半圆. 为方便分析, 周期数应为整数, 此时单胞对应圆心角可选值如下列集合所示: $ \left\{ {{90}^ \circ },{{60}^ \circ },{{45}^ \circ }, {{30}^ \circ } ,\cdots , {{20}^ \circ } \right\} $.
另外, 该类结构的调制作用主要源于叶片与轮盘之间的动力学特性差异, 故单胞基体所对应的圆心角不可过小或过大. 根据有限元方法的模拟和分析, 在圆心角处于上述集合的第4至最后一项时, 不存在期望的区域振动, 叶片主导的振动模态被叶-盘耦合模态取代, 所以圆心角的最佳取值为$ {60^ \circ } $或$ {45^ \circ } $.
柱坐标xoy平面内的周向基矢量
$$ {{\boldsymbol{e}}_\varphi } = - {{\boldsymbol{e}}_x}\sin \theta + {{\boldsymbol{e}}_y}\cos \theta $$ (1) 该矢量为变矢量, 数值计算时角度参数可表示为单胞圆心角乘以参数化变量, 参数化变量一般取为k = (0:0.1:1). 由于轮盘上的叶片数为18, 因此在任意$ {45^ \circ } $倍数的单胞下, 均导致无法阵列得到原始叶轮, 其边界无法完全复制, 因此圆心角选择为$ {60^ \circ } $.
1.2 结构的力学模型构建
为建立基于有限元方法的数值分析方法, 以基梁-立柱模型的弯曲振动为例, 开展相应的力学建模和解析分析方法研究. 设置其周期延拓方向为x轴, 面内垂直轴为y轴, 垂直于xoy平面的面外方向为z轴. 理论上, 基梁存在垂直于传播方向(x方向)的对称振动平面, 即振幅在yoz平面内, 此时该弯曲波是非偏振的. 但实际梁振动模型中, y方向为主要弯曲振动模式, 因此二维模型下计算域在xoy平面内, 在y方向具有波传播并使用波方程描述, 而忽略其他方向的振幅, 且对于所研究的叶轮模型, 其边界部分的周向特征尺度为5π/144 mm, 厚度特征尺度为5 mm. 因此对于该类周期梁, 可按遵循以下假设的欧拉梁模型来进行解析计算:
(1) x方向尺寸远大于厚度方向和梁宽度的尺寸, 且挠度在z方向变形均匀;
(2) 变形前, 垂直于中性面的截面, 变形后仍垂直于中性面;
(3) 由上述条件可知, 截面不考虑剪切应力, 结合应力平衡条件, y方向各层间的横向应力为0.
基于以上假设的力学模型为Euler-Bernoulli梁, 弯曲波控制方程为
$$ YI\frac{{{\partial ^4}w}}{{\partial {x^4}}} + \rho A\frac{{{\partial ^2}w}}{{\partial {t^2}}} = 0 $$ (2) 其中, $ Y $为杨氏模量, $ I $为梁的截面惯性矩, $ \rho $为密度, $ A $为截面积, $ w $为横向位移即挠度.
考虑更一般的情况, 例如板结构的弯曲振动, 首先确定其主要弯曲振动模式. 与梁的振动模式对应, 板结构也具有扭转振动、径向振动和横向振动. 根据模态分析, 可知其主要振动模式为横向振动, 此时待解参数为挠度. 以Kirchhoff板基本假设进行推导, 对于板的界定或定义一般以厚宽比$ \delta $来确定
$$ \left\{\begin{aligned} & \delta < 1/6,\quad {\text{Kirchhoff plate}} \\ & 1/6 < \delta < 1/2,\quad {\text{Mindlin plate}} \end{aligned}\right. $$ 对于所研究叶轮的轮盘结构, 其盘厚度与参与弯曲振动部分的径向宽度比值小于1/6, 故在理论推导中应采用Kirchhoff板假设. 图2为轮盘结构图(不含叶片), 可以看到轮盘厚度在半径的2/3 ~ 1/2处激增, 因此整个轮盘可近似看作一个圆形薄板和一个厚空心圆柱体. 根据模态分析结果, 轮盘主导的模态与含中心孔的等厚度薄圆板相似. 因此, 下节将把轮盘等效视为薄圆板, 给出内环面约束条件, 并且由于扭曲叶片结构复杂, 采用文献[29]中的假设, 将叶片等效为Euler梁模型. 同时, 若将轮缘及此处的叶片振动模式投影至R = 0.125 m的柱面并展开, 其与周期梁模型下的弯曲振动类似, 故下节中将同时从该假设进行推导, 并对比Kirchhoff板模型与Euler梁模型的计算结果.
2. 结构的波动方程及传递矩阵
2.1 梁-立柱模型
对于梁的横向自由振动
$$ {\nabla ^4}w + \frac{{\rho {{S}}}}{{YI}}\frac{{{\partial ^2}w}}{{\partial {t^2}}} = 0 $$ (3) 纵波控制方程
$$ Y\frac{{{\partial ^2}u}}{{\partial {{{s}}^2}}} - \rho \frac{{{\partial ^2}u}}{{\partial {t^2}}} = 0 $$ (4) 其中, S为梁的横截面面积; 定义s方向为轮盘的环向, 弯曲波控制方程中的4阶算子同样定义在环向. 基于Euler梁的力学假设, 可设其通解为: $ \dfrac{{\partial w}}{{\partial r}} = 0 $时, 挠度解及纵向位移解为
$$ w(s) = \left( {{A_i}{{\mathrm{e}}^{{\mathrm{i}}{k_1}{{s}}}} + {B_i}{{\mathrm{e}}^{ - {\mathrm{i}}{k_1}{{s}}}} + {C_i}{{\mathrm{e}}^{{k_1}{{s}}}} + {D_i}{{\mathrm{e}}^{ - {k_1}{{s}}}}} \right){{\mathrm{e}}^{{\mathrm{i}}\omega t}} $$ (5) $$ u(s) = \left( {{E_i}{{\mathrm{e}}^{{\mathrm{i}}{k_2}{{s}}}} + {F_i}{{\mathrm{e}}^{ - {\mathrm{i}}{k_2}{{s}}}}} \right){{\mathrm{e}}^{{\mathrm{i}}\omega t}} $$ (6) 弯曲波解与纵波解相差一对非波解, 即振荡项. 式中$ {k_1^4} = \rho {\omega ^2}{{S}}/\left(YI \right) $和$ {k_2^2} = \rho {\omega ^2}/Y $分别为弯曲波波数及纵波波数.
假设弯曲波在径向分布均匀, 可将沿周向的薄圆环展开为链式结构. 叶片高度折算为平均等直梁, 单胞的特征长度取为$ {{s}} = 5\text{π} r/54 $.
下面推导单胞的传递矩阵. 取s方向为顺时针方向, 根据圆盘的对称性, 仅分析其1/2即可. 展开后取波传播方向为正方向, 并定义单胞内波从左端面传递至右端面.
定义各截面的挠度、纵向位移、剪切力、纵向力、弯矩、坡度共同构成状态向量
$$ {{\boldsymbol{v}}_i} = {\left( {{w_i},{u_i},{F_i},{T_i},{M_i},\alpha } \right)^\text{T}} $$ (7) 系数向量
$$ {{\boldsymbol{G}}_i} = {\left( {{A_i},{B_i},{C_i},{D_i},{E_i},{F_i}} \right)^\text{T}} $$ (8) 在某一截面上, 各变量相位相同, 可得状态向量和系数向量的关系
$$ {{\boldsymbol{v}}_i} = {{\boldsymbol{N}}_1}{{\boldsymbol{G}}_i} $$ (9) 其中, i表示单胞内左基梁、右基梁和等直立柱, 分别取值为(i = 1, 2, 3). 在波传播方向上, 左(右)基梁的两个截面上存在下列关系式
$$ {\boldsymbol{v}}_j^R = {{\boldsymbol{N}}_2}{\boldsymbol{G}}_j^L $$ (10) 其中, j = 1, 2; R和L分别代表单胞内某基梁的右、左端面; $ {{\boldsymbol{N}}_2} $为两截面间的传递矩阵. 受力情况与文献[30]规定一致, 如图3所示.
联立内部力平衡条件、边界条件和力矩平衡条件可得内部各截面间的系数向量关系
$$ \left.\begin{split} & {\boldsymbol{G}}_2^L = {{\boldsymbol{N}}_3}{\boldsymbol{G}}_1^R \\ & {{\boldsymbol{G}}_p} = {{\boldsymbol{N}}_4}{\boldsymbol{G}}_1^R\end{split}\right\} $$ (11) 进一步, 联立方程式(9) ~ 式(11), 可得
$$ {\boldsymbol{v}}_2^R = {{\boldsymbol{N}}_2}{{\boldsymbol{N}}_3}{\boldsymbol{N}}_1^{ - 1}{{\boldsymbol{N}}_2}{\boldsymbol{N}}_1^{ - 1}{\boldsymbol{v}}_1^L $$ (12) $$ {\boldsymbol{M}} = {{\boldsymbol{N}}_2}{{\boldsymbol{N}}_3}{\boldsymbol{N}}_1^{ - 1}{{\boldsymbol{N}}_2}{\boldsymbol{N}}_1^{ - 1} $$ (13) $$ {{\boldsymbol{N}}_1} = \left( {\begin{array}{*{20}{c}} 1&1&1&1&0&0 \\ 0&0&0&0&1&1 \\ 0&0&0&0&{{K_1}}&{ - {K_1}} \\ {{\mathrm{i}}{K_2}}&{ - {\mathrm{i}}{K_2}}&{ - {K_2}}&{{K_2}}&0&0 \\ { - {K_3}}&{ - {K_3}}&{{K_3}}&{{K_3}}&0&0 \\ {{\mathrm{i}}{k_1}}&{ - {\mathrm{i}}{k_1}}&{{k_1}}&{ - {k_1}}&0&0 \end{array}} \right) $$ (14) 其中, $ {K_1} = Y{{S}}{\mathrm{i}}{k_2} $, $ {K_2} = YIk_1^3 $, $ {K_3} = YIk_1^2 $.
由于$ {{\boldsymbol{N}}_2} $形式复杂, 采用分块矩阵来表示, 即原$ 6 \times 6 $矩阵写为两个$ 6 \times 3 $矩阵
$$ {{\boldsymbol{N}}_2} = \left( {\begin{array}{*{20}{c}} {{{\boldsymbol{P}}_1}}&{{{\boldsymbol{P}}_2}} \end{array}} \right) $$ (15) $$ {{\boldsymbol{P}}_1} = \left( {\begin{array}{*{20}{c}} {{{\mathrm{e}}^{{\mathrm{i}}{k_1}D/2}}}&{{{\mathrm{e}}^{ - {\mathrm{i}}{k_1}D/2}}}&{{{\mathrm{e}}^{{k_1}D/2}}} \\ 0&0&0 \\ 0&0&0 \\ {{\mathrm{i}}{K_2}{{\mathrm{e}}^{{\mathrm{i}}{k_1}D/2}}}&{ - {\mathrm{i}}{K_2}{{\mathrm{e}}^{ - {\mathrm{i}}{k_1}D/2}}}&{ - {K_2}{{\mathrm{e}}^{{k_1}D/2}}} \\ { - {K_3}{{\mathrm{e}}^{{\mathrm{i}}{k_1}D/2}}}&{ - {K_3}{{\mathrm{e}}^{ - {\mathrm{i}}{k_1}D/2}}}&{{K_3}{{\mathrm{e}}^{{k_1}D/2}}} \\ {{\mathrm{i}}{k_1}{{\mathrm{e}}^{{\mathrm{i}}{k_1}D/2}}}&{ - {\mathrm{i}}{k_1}{{\mathrm{e}}^{ - {\mathrm{i}}{k_1}D/2}}}&{{k_1}{{\mathrm{e}}^{{k_1}D/2}}} \end{array}} \right) $$ (16) $$ {{\boldsymbol{P}}_2} = \left( {\begin{array}{*{20}{c}} {{{\mathrm{e}}^{ - {k_1}D/2}}}&0&0 \\ 0&{{{\mathrm{e}}^{{\mathrm{i}}{k_2}D/2}}}&{{{\mathrm{e}}^{ - {\mathrm{i}}{k_2}D/2}}} \\ 0&{{K_1}{{\mathrm{e}}^{{\mathrm{i}}{k_2}D/2}}}&{ - {K_1}{{\mathrm{e}}^{ - {\mathrm{i}}{k_2}D/2}}} \\ {{K_2}{{\mathrm{e}}^{ - {k_1}D/2}}}&0&0 \\ {{K_3}{{\mathrm{e}}^{ - {k_1}D/2}}}&0&0 \\ { - {k_1}{{\mathrm{e}}^{ - {k_1}D/2}}}&0&0 \end{array}} \right) $$ (17) 自由振动条件、应力连续性条件、自由边界条件可表示为
$$ \left.\begin{split} & \sum\limits_j F = 0 \\ & \sum M = 0 \\ & \sigma _{jk}^n = \sigma _{jk}^{n + 1} \\ & u_{jk}^n = u_{jk}^{n + 1} \\ & {\sigma _{{\mathrm{free}}}} = {\tau _{{\mathrm{free}}}} = {\kappa _{{\mathrm{free}}}}{\text{ = 0}} \end{split}\right\} $$ (18) 其中, $ j和k $分别代表截面法向及力的作用力方向, $ n $和$ n + 1 $代表相邻截面, 自由面为叶片上端面.
在周期结构中, Bloch定理可以写为
$$ {\boldsymbol{v}}_2^R = {{\mathrm{e}}^{{\mathrm{i}}ks}}{\boldsymbol{v}}_1^L $$ (19) 联立式(12)和式(19), 可得齐次方程及有非零解条件
$$\qquad\qquad\quad ({\boldsymbol{M}} - {{\mathrm{e}}^{{\mathrm{i}}ks}}{\boldsymbol{E}}){\boldsymbol{v}}_1^L = {\boldsymbol{0}} $$ (20) $$ \qquad\qquad\quad \left| {{\boldsymbol{M}} - {{\mathrm{e}}^{{\mathrm{i}}ks}}{\boldsymbol{E}}} \right| = 0 $$ (21) 其中, k为Bloch波矢, 求解行列式方程即可得到波矢与频率间的关系.
取前3阶结果, 如图4所示. 图中纵坐标为频率, 横坐标为无量纲波数. 每阶频率都取10个离散点计算, 横坐标小于0.5时, 第一阶为弯曲波模式, 第三阶为纵波模式. 很明显, 图中前3阶无带隙.
![]()
图 4 波矢与频率间的关系(前3阶结果)Figure 4. Relationship between wave vector and frequency (first third order)在研究局部化问题时, 局部化因子同样是一项重要参数, 而根据动力系统中Lyapunov指数定义, 可得局部化因子的一般表达式
$$\qquad\qquad \gamma = \mathop {\lim }\limits_{n \to \infty } \frac{1}{n}\ln \sum\limits_{i = 1}^n {\left\| {\frac{{{{\boldsymbol{T}}_i}{{\boldsymbol{v}}_{i - 1}}}}{{\left\| {{{\boldsymbol{v}}_{i - 1}}} \right\|}}} \right\|} $$ (22) $$\qquad\qquad {{\boldsymbol{M}}_i}{\boldsymbol{v}}_i^L = {\boldsymbol{v}}_i^R = {\boldsymbol{v}}_{i + 1}^L $$ (23) 根据文献[31]给出的表达式, 可快速计算局部化因子
$$ {\gamma _m} = \mathop {\lim }\limits_{n \to \infty } \frac{1}{n}\ln \sum\limits_{i = 1}^n {\left\| {{\lambda _m}} \right\|} $$ (24) 其中, m为传递矩阵阶数即状态向量元素数. 利用式(22)计算时, 需选取m个正交标准化初始状态向量, 利用式(23)进行迭代, 迭代至待求单胞时, 对m个右状态向量进行标准正交化. 虽然表达式为极限形式, 但对有限周期可进行数值平均, 从而求得局部化因子.
图5为局部化因子的计算结果, 由于频率决定了弯曲波波数和纵波波数, 故当频率取0时, 局部化因子趋于不可取的极大值. 低频存在区域极大值, 这与色散曲线中低频不存在带隙结果不一致. 其原因是: 由于本文所讨论的周期结构的结构数低于10, 使得由平均因子所描述的波传播特性难以精确描述每个周期内的传播特性. 同时, 由于计算时, 为了考虑叶片高度沿其型线不断变化的特性, 对叶片高度进行了平均, 使得在所研究的宽度上, 叶片等效弯曲刚度降低, 弯曲振动模式在低频更易被激发. 由于叶片平均高度为梁厚度的3倍, 当叶片在低频被激发振动时, 将对梁产生剪切力, 抑制梁在低频的压缩波传播, 从而产生局部化因子低频极大值现象. 在高频位置, 存在弹性波传播受阻的频带, 这与实际情况相符合.
综上, 若采用Euler-Bernoulli梁的力学模型推导离心叶轮的周向传递矩阵, 其局部化因子计算结果振荡. 这意味着结构不存在稳定的通带, 与实际情况不符, 因此2.2节中将选取更复杂的力学模型来研究该问题.
2.2 Kirchhoff板模型
若考虑条件$ w'_r \ne 0 $, 则需在空间项中增加一个维度. 为简化计算, 将采用不考虑剪切效应并满足直法线假设的Kirchhoff板模型, 由于本文所研究对象为叶轮, 故将方程变换至柱坐标系
$$ \qquad\qquad {\nabla ^2}{\nabla ^2}w + \frac{{\rho h}}{D}\frac{{{\partial ^2}w}}{{\partial {t^2}}} = 0 $$ (25) $$ \qquad\qquad {\nabla ^2} = \frac{{{\partial ^2}}}{{\partial {r^2}}} + \frac{1}{r}\frac{\partial }{{\partial r}} + \frac{1}{{{r^2}}}\frac{{{\partial ^2}}}{{\partial {\theta ^2}}} $$ (26) $$ \qquad\qquad D = \frac{{E{h^3}}}{{12(1 - {v^2})}} $$ (27) 式中, $ D $为圆板的弯曲刚度, $ v $为泊松比.
圆盘具有内边界固定和外自由边界的条件, 且自由振动关于$ \theta = 0 $对称, 与梁类似, 取空间相位与时间相位相乘的分离变量解
$$\qquad\qquad w = G(\omega t)W(r,\theta ) $$ (28) $$\qquad\qquad W(r,\theta ) = F({k_b}r)\cos n\theta $$ (29) $$ \qquad\qquad\begin{split} & F({k_b}r) = A{{\mathrm{J}}_n}({k_b}r) + B{{\mathrm{N}}_n}({k_b}r) + \\ &\qquad C{{\mathrm{I}}_n}({k_b}r) + D{{\mathrm{K}}_n}({k_b}r) \end{split} $$ (30) 式中, $ {{\mathrm{J}}_n}({k_b}r) $和$ {{\mathrm{N}}_n}({k_b}r) $为实宗量的n阶第一、第二类Bessel函数; $ {{\mathrm{I}}_n}({k_b}r) $和$ {{\mathrm{K}}_n}({k_b}r) $为虚宗量的n阶第一、第二类Bessel函数. 而弯曲波波数表示为: $ {k_b^4} = \rho {\omega ^2}h/D $.
仅考虑横向振动时, 定义状态向量为
$$ {{\boldsymbol{v}}_i} = {\left( {{w_i},{T_i},{M_i},{\alpha _i}} \right)^\text{T}} $$ (31) 薄板上的叶片与前文一致, 等效为梁结构, 此时板将产生扭转波与弯曲波的耦合模态. 若将连接条件设置为简支, 由力平衡条件则叶片的自由边界将产生弯矩与剪力. 同时, 叶片的挠度即弯曲振动也不由板内周向力激励, 而是由弯曲振动产生的偏心力激励. 其波方程为
$$ YI\frac{{{\partial ^4}{w_p}}}{{\partial {z^4}}} + \rho A\frac{{{\partial ^2}{w_p}}}{{\partial {t^2}}} = 0 $$ (32) $$ Y\frac{{{\partial ^2}{u_p}}}{{\partial {z^2}}} - \rho \frac{{{\partial ^2}{u_p}}}{{\partial {t^2}}} = 0 $$ (33) $$ {w_{pi}} = {{\mathrm{e}}^{{\text{i}}\omega t}}(E_1^i{{\mathrm{e}}^{{\text{i}}kz}} + E_2^i{{\mathrm{e}}^{ - {\text{i}}kz}} + E_3^i{{\mathrm{e}}^{kz}} + E_4^i{{\mathrm{e}}^{ - kz}}) $$ (34) $$ {u_{pi}} = {{\mathrm{e}}^{{\text{i}}\omega t}}(F_1^i{{\mathrm{e}}^{{\text{i}}{k_2}z}} + F_2^i{{\mathrm{e}}^{ - {\text{i}}{k_2}z}}) $$ (35) $$ {k^4} = \rho {\omega ^2}{{S}}/\left(YI \right) $$ (36) $$ {k_2^2} = \rho {\omega ^2}/Y $$ (37) $$ z = 0:\left\{ \begin{aligned} & {w_p} = 0 \\ & \kappa = 0 \end{aligned} \right. $$ (38) 对单胞中心节点, 存在平衡条件
$$ \left.\begin{split} & {w_1} = {w_2} = - {u_p} \\ & EA\frac{{\partial {u_p}}}{{\partial z}} - D\left(d + \frac{1}{{{r_0}}}c - \frac{1}{{{r_0}^2}}b + 2\frac{{{n^2}}}{{{r_0}^3}}\right)\left| {\begin{array}{*{20}{c}} {w = {w_1}} \\ {r = {r_0}} \end{array}} \right. = \\ &\qquad - D\left(d + \frac{1}{{{r_0}}}c - \frac{1}{{{r_0}^2}}b + 2\frac{{{n^2}}}{{{r_0}^3}}\right)\left| {\begin{array}{*{20}{c}} {w = {w_2}} \\ {r = {r_0}} \end{array}} \right. \\ & EI\frac{{{\partial ^2}{w_p}}}{{\partial {z^2}}} - D\left(c + \frac{b}{{{r_0}}} - \frac{{{n^2}}}{{{r_0}^2}}\right)\left| {\begin{array}{*{20}{c}} {w = {w_1}} \\ {r = {r_0}} \end{array}} \right. = \\ &\qquad - D\left(c + \frac{b}{{{r_0}}} - \frac{{{n^2}}}{{{r_0}^2}}\right)\left| {\begin{array}{*{20}{c}} {w = {w_2}} \\ {r = {r_0}} \end{array}} \right. \\ & \frac{{\partial {w_p}}}{{\partial z}} = b\left| {\begin{array}{*{20}{c}} {w = {w_1}} \\ {r = {r_0}} \end{array}} \right. = b\left| {\begin{array}{*{20}{c}} {w = {w_2}} \\ {r = {r_0}} \end{array}} \right. \end{split}\right\} $$ (39) 各状态参数表达式如下
$$ w = G(\omega t)F({k_b}{r_0})({K_1}\cos n\theta + {K_2}\sin n\theta ) $$ (40) $$ \alpha (\theta ) = {{\partial w} {\partial r}}\big| {r = {r_0}} $$ (41) $$ M(\theta ) = - D{\nabla ^2}w\big| {r = {r_0}} $$ (42) $$ T(\theta ) = - D\frac{\partial }{{\partial r}}({\nabla ^2}w)\big| {r = {r_0}} $$ (43) 其中, $ a = F({k_b}r)\big| {r = }{r_0}, $ $ b = \dfrac{{\partial F({k_b}r)}}{{\partial r}}\Big| {r = } {r_0}, $ $ c = \dfrac{{{\partial ^2}F({k_b}r)}}{{\partial {r^2}}}\Big| {r = } {r_0} ,$ $ d = \dfrac{{{\partial ^3}F({k_b}r)}}{{\partial {r^3}}}\Big| {r = } {r_0}. $
基于上述理论进行计算时, 由于弯曲波主要沿周向传播, 因此可认为其具有一维周期性. 采用与文献[32]类似的处理, 在某一半径的邻域内选择平均半径$ {r_0} $作为计算的特征参数, 并认为在$ {r_0} $附近的环周期在半径方向挠度分布相同, 如图6所示. 这样的假设使得可以通过分析某一半径领域附近的环周期得到整个扇形的弹性波周期传播特性. 对应的拉普拉斯算子的径向导数为$ {r_0} $处的导函数值, 此时为常数.
基于上述条件和公式, 可推导出每一周期左右截面间的传递矩阵. 首先, 将与转角有关的解部分写为更一般的形式(仅考虑空间解系, 忽略时间因子)
$$ w = F({k_b}{r_0})({K_1}\cos n\theta + {K_2}\sin n\theta ) $$ (44) 式中, $ {K_1} $和$ {K_2} $为解系数, 满足关系$ {K_1}^2 + {K_2}^2 = 1 $, $ {r_0} $为所研究的环向对应的某一确定半径, 为常数. 由于解形式为2阶, 故状态向量取两个二维向量.
下面给出状态向量与对应位置系数向量的关系
$$\qquad\quad \left.\begin{aligned} & {{\boldsymbol{v}}_{i,1}} = {\left( {{w_i},{\alpha _i}} \right)^\text{T}} \\ & {{\boldsymbol{v}}_{i,2}} = {({T_i},{M_i})^\text{T}} \end{aligned}\right\} $$ (45) $$\qquad\quad {{\boldsymbol{v}}_{i,j}} = {{\boldsymbol{N}}_{1,j}}{\left( {\cos n\theta ,\sin n\theta } \right)^\text{T}}{{\boldsymbol{K}}_{i,j}} $$ (46) 其中, $ i $为截面序号, $ j $为状态向量序号.
$ {{\boldsymbol{N}}_{1,1}} $可写为
$$ {{\boldsymbol{N}}_{1,1}} = \left( {\begin{array}{*{20}{c}} a&a \\ b&b \end{array}} \right) $$ (47) 与梁-立柱模型相似, 将单胞分为三部分, 第一部分中考虑单胞左截面系数矩阵与中截面状态向量的关系
$$ {\boldsymbol{v}}_{i,j}^R = {{\boldsymbol{N}}_{2,j}} \cdot \left( {\cos n\theta ,\sin n\theta } \right) \cdot {{\boldsymbol{K}}^L_{i,j}} $$ (48) 式中, $ {{\boldsymbol{N}}_{2,j}} $为右边界状态向量与左边界系数向量间的传递矩阵.
设空间单胞相位因子为$ {\theta {'}} $, 此时考虑同横向位移与系数向量的关系
$$\qquad w_i^R = \left( {\begin{array}{*{20}{c}} {\cos {\theta {'}}}&{\sin {\theta {'}}} \\ { - \sin {\theta {'}}}&{\cos {\theta {'}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{K_1}} \\ {{K_2}} \end{array}} \right) $$ (49) $$ \qquad {\theta {'}} = \frac{{{k_b}{r_0}{\theta _0}}}{2} $$ (50) 在最关心的4节径横向振动中, 此时取n = 4, $ \omega = \left( {4500, 5000, 5500} \right){\text{Hz}} $. 由式(49)计算的单胞圆心角为$ {\theta _0} = \left( {{{63}^ \circ }{{,60}^ \circ }{{,57}^ \circ }} \right) $. 可将周期单胞确定为$ {60^ \circ } $, 这样在有限元计算中运用Bloch定理时, 可以保证计算精度.
同时, 观察叶轮的单胞发现, 主要起调制作用的为一个扭曲叶片, 故在运用Bloch定理计算时, 按实际等直梁数计算; 而在局部化因子推导时仅考虑一个立柱. 另外, 还需注意在式(49)中, 极坐标参量$ {\theta _0} $并非单周期的波空间相位, 而是表示旋转的相位因子, 所以
$$ n{\theta {'}} \ne 0.5{\theta _0} $$ (51) 同理可以推导出叶片上的相位矩阵, 此处不再展开推导. 定义沿周向的Bloch波矢, 并在每个单胞边界运用Bloch定理
$$ {\boldsymbol{v}}_i^R = {\boldsymbol{v}}_{i + 1}^L = {{\mathrm{e}}^{{\text{i}}k{r_0}{\theta _0}}}{\boldsymbol{v}}_i^L $$ (52) 根据单胞边界连续性条件, 给定n即可数值计算对应阶的色散曲线.
关于局部化因子的计算, 可按以下步骤进行. 由式(48)且假设初始截面为$ \theta = 0 $, 可得$ {{\boldsymbol{N}}_{2,j}} $
$$ {{\boldsymbol{N}}_{2,{\text{1}}}} = \left( {\begin{array}{*{20}{c}} {a\left( {\cos {\theta {'}} - \sin {\theta {'}}} \right)}&{a(\sin {\theta {'}} + \cos {\theta {'}})} \\ {b( - \cos {\theta {'}} - \sin {\theta {'}})}&{b(\cos {\theta {'}} - \sin {\theta {'}})} \end{array}} \right) $$ (53) $$ {{\boldsymbol{N}}_1} = \left( {\begin{array}{*{20}{c}} a&a \\ { - D\left(d + \dfrac{1}{{{r_0}}}c - \dfrac{1}{{{r_0}^2}}b + 2\dfrac{{{n^2}}}{{{r_0}^3}}\right)}&{ - D\left(d + \dfrac{1}{{{r_0}}}c - \dfrac{1}{{{r_0}^2}}b + 2\dfrac{{{n^2}}}{{{r_0}^3}}\right)} \\ { - D\left(c + \dfrac{b}{{{r_0}}} - \dfrac{{{n^2}}}{{{r_0}^2}}\right)}&{ - D\left(c + \dfrac{b}{{{r_0}}} - \dfrac{{{n^2}}}{{{r_0}^2}}\right)} \\ b&b \end{array}} \right) $$ $ {{\boldsymbol{N}}_{1,1}} $为第1, 4行; $ {{\boldsymbol{N}}_{1,2}} $为第2, 3行
$$ {{\boldsymbol{N}}_{2,2}} = \left( {\begin{array}{*{20}{c}} { - D\left(d + \dfrac{1}{{{r_0}}}c - \dfrac{1}{{{r_0}^2}}b + 2\dfrac{{{n^2}}}{{{r_0}^3}}\right)}&{ - D\left(d + \dfrac{1}{{{r_0}}}c - \dfrac{1}{{{r_0}^2}}b + 2\dfrac{{{n^2}}}{{{r_0}^3}}\right)} \\ { - D\left(c + \dfrac{b}{{{r_0}}} - \dfrac{{{n^2}}}{{{r_0}^2}}\right)}&{ - D\left(c + \dfrac{b}{{{r_0}}} - \dfrac{{{n^2}}}{{{r_0}^2}}\right)} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {\cos {\theta {'}} - \sin {\theta {'}}}&0 \\ 0&{\sin {\theta {'}} + \cos {\theta {'}}} \end{array}} \right) $$ 与前文相似, 根据状态向量间的关系及连续条件可求得总传递矩阵$ {{\boldsymbol{M}}_j} $, 由于传递矩阵为2阶, 故选取一标准正交初始向量$ {\boldsymbol{u}}_{1,j}^0 $, 经过m次迭代得到状态向量, 结合式(24)选择较小的结果可得
$$ {\boldsymbol{v}}_1^{m + 1} = {M_{m + 1}}{\boldsymbol{u}}_1^m\quad (m = 0,1,2,\cdots, N) $$ (54) 由此可得如图7所示的局部化因子图.
图7为n = 4时的局部化因子图, 可以看出: 采用板模型后, 由于局部化因子按n来计算, 因此n = 4时低频的波峰被消除. 在4节径振动部分的局部化因子较大, 此时若在时间历程上依次输入2个对应波长的激励波(左右行波共4个输入波), 可以发现振动响应在周向的传播均匀并形成驻波, 叶片参与度较大(即叶片被激发对应的弯曲模态, 以弱耦合系统形式在随叶盘的振动, 而运动的同时自身也进行弯曲振动).
该节主要讨论了以传递矩阵描述的解析解, 可以看到在叶轮的振动问题中, 虽然采取板方程描述的弹性波传输特性更为精准, 但仍未能准确描述离心叶轮类循环对称结构中普遍存在的节径型振动模式应具有的“驻波”特性. 因此在下一节, 将针对该问题利用有限元方法进行求解计算并分析其数值结果.
3. 基于有限元的求解方法及数值结果
3.1 有限元方法求解弹性波动方程
对于一般动力系统, 其动力学方程
$$\qquad {\boldsymbol{Mx''}} + {\boldsymbol{Cx'}} + {\boldsymbol{Kx}} = {\boldsymbol{Q}} $$ (55) $$\qquad {\boldsymbol{x}} = \left( {\begin{array}{*{20}{c}} {{x_1}} \\ {{x_2}} \\ {\begin{array}{*{20}{c}} {. ..} \\ {{x_n}} \end{array}} \end{array}} \right),\quad {\boldsymbol{Q}} = \left( {\begin{array}{*{20}{c}} {{q_1}} \\ {{q_2}} \\ {\begin{array}{*{20}{c}} {. ..} \\ {{q_n}} \end{array}} \end{array}} \right) $$ (56) 式中$ {\boldsymbol{M}} $, $ {\boldsymbol{C}} $, $ {\boldsymbol{K}} $和$ {\boldsymbol{Q}} $分别为系统质量矩阵、阻尼矩阵、刚度矩阵和外力矩阵. 在自由振动问题中, 方程右边为0. 若考虑简谐形式通解, 则有特征方程
$$ \left( {{\boldsymbol{K}} - {\omega ^2}{\boldsymbol{M}}} \right){\boldsymbol{x}} = {\boldsymbol{0}} $$ (57) 其中位移向量为广义向量. 单元内某点处的位移可由节点位移表示为
$$ {{\boldsymbol{x}}_e} = {\boldsymbol{\tilde Nx}} $$ (58) 多自由度方程通过解耦成为多个单自由度节点方程, 刚度矩阵、质量矩阵、载荷向量正交变换为
$$\qquad\qquad\quad {{\boldsymbol{M}}_n} = {{\boldsymbol{A}}^\text{T}}{\boldsymbol{MA}} $$ (59) $$\qquad\qquad\quad {{\boldsymbol{K}}_n} = {{\boldsymbol{A}}^\text{T}}{\boldsymbol{KA}} $$ (60) $$\qquad\qquad\quad {{\boldsymbol{Q}}_n} = {{\boldsymbol{A}}^\text{T}}{\boldsymbol{Q}} $$ (61) 变换向量$ {\boldsymbol{A}} $为简谐形式振幅向量, 满足正交条件; $ {{\boldsymbol{M}}_n} $和$ {{\boldsymbol{K}}_n} $分别为正定或半正定矩阵. 通过解式(56), 可得到n阶固有频率和模态振型.
对周期结构运用Bloch定理, 可以得到随各方向波矢量变化的各模态矩阵, 进一步求解可得色散曲线即能带结构.
叶轮为典型的循环周期结构, 沿界面传播的弯曲波、扭转波的波矢量为周向和径向, 此时需与笛卡尔坐标系下的波矢量进行基变换. 由式(1)知
$$ {{\boldsymbol{k}}_{\boldsymbol{\varphi }}} = - \sin {\theta _0}{{\boldsymbol{k}}_x} + \cos {\theta _0}{{\boldsymbol{k}}_y} $$ (62) 设置参数化扫描k = (0:0.05:1), $ {\theta _0} $为单胞圆心角, 即可得到色散曲线.
3.2 能带结构的数值结果
图8给出了以$ {60^ \circ } $为圆心角的单胞能带结构. 从图8可以看到, 轮盘所属低阶曲线频散性极弱, 群速度变化很小; 叶片部分则出现群速度变化率较高的色散曲线. 为了区别振动模式, 按文献[33]给出的方法来确定模式极化率
![]()
图 8 以$ {60^ \circ } $为圆心角的单胞能带结构Figure 8. Single cell band structure with a center angle of $ {60^ \circ } $$$ \beta = \frac{{\displaystyle\iiint {{u_3}u_3^*{\mathrm{d}}V}}}{{\displaystyle\iiint {\sum\limits_{i = 1}^3 {{u_i}u_i^*} {\mathrm{d}}V}}} $$ (63) 上式可用来识别平板类声子晶体的面外模式及面内模式. 最新研究表明, 在具有弯曲边界的薄板结构中, 面外与面内模式耦合[34], 且振动模式常对应偏振范围内的耦合振动模式, 这使得极化率除了文献所定义的识别范围外, 还存在$ 0.01 \leqslant \beta \leqslant 0.94 $的情况, 这也是弱耦合系统的典型特性. 另外, 在叶片的弯曲振动时, 其面外位移场累加变量为xz或yz耦合, 故某些偏振模式可以识别为弯曲模式, 此时极化率的弯曲振动模式识别范围应适当改变. 图9为极化率示意图, 颜色越接近红色意味着极化率越大, 反之极化率越小.
图9中存在很明显的模式转换现象. 不包括曲线斜率为0的局域化振动区域, 其他部分均存在面内模式、耦合模式、面外模式的转换. 文献[33]给出的识别范围为: $ \beta \leqslant 0.01 $振动模式为面内振动; $ \beta \geqslant 0.94 $振动模式为面外振动. 结合图9和过去对离心叶轮模态振型的研究结果, 给出本研究涉及循环周期结构的振动模式识别范围: $ \beta < 0.124 $时为严格面内模式, 叶片基本不参与弯曲振动, 轮盘面内扭转为主振动模式; $ 0.124 \leqslant \beta < 0.361 $时为轮盘和叶片的耦合弯曲振动模式, 叶片的弯曲振动参与度大于轮盘的弯曲振动参与度; $ 0.361 \leqslant \beta < 0.599 $为严格叶片弯曲振动模式, 叶盘不参与弯曲振动. 其极化率比轮盘的面外模式极化率低的原因为叶片的弯曲振动同时包含xy和z方向的贡献; $ 0.599 \leqslant \beta < 0.718 $时为轮盘和叶片的耦合弯曲振动模式, 轮盘的弯曲振动参与度大于叶片的弯曲振动参与度; $ 0.718 \leqslant \beta < 0.955 $时为轮盘严格面外模式, 此时叶片几乎不参与弯曲振动; 而当$ \beta > 0.955 $时整个单胞处于剧烈弯曲振动模式, 轮盘的横向位移远超轮缘厚度.
确定极化率后, 可以通过分析带隙两侧的振动模式来确定带隙的类型及形成机理.
第一个带隙位于2762.4 ~ 3111.6 Hz, 下边界为轮盘主导的扭转振动模式, 上边界为轮盘1节径横向振动模式. 此处带隙应为“Bragg带隙”, 产生的原因在于圆盘的任意一条直径为完全对称轴, 同时满足轴对称和中心对称. 在圆盘圆心角大于180°时, 就会产生上述带隙, 其来自于叠加驻波解. 事实上若仅计算锐角单胞, 0斜率振动曲线以下的各阶振动曲线会形成一条斜率不断变化的曲线来替代该类带隙(不一定全部消失). 虽然我们计算时选取的另一个单胞同样非周向连通, 但它已经可以表现出左右行波在周向叠加为驻波模态的性质. 由此我们可以大胆假设当波在循环对称结构中传播时, 并非传统意义上人们所认识的沿周向无“物理边界”, 而是以完全贯穿结构的直径或半径作为对称轴时, 如果此轴同时满足两种对称性, 那么我们可计算得到其具有的驻波带隙. 换句话说, 非连通周向循环结构的物理反射壁面在180°附近.
第二个带隙位于3115.8 ~ 3773.9 Hz,下边界为轮盘0波数横向振动模式, 上边界为1节径轮盘横向振动模式. 带隙类型与上述带隙相同, 为驻波模态型的Bragg带隙.
第三带隙处于3863.9 ~ 4743.6 Hz, 下边界为2节径轮盘横向振动模式, 上边界为4节径轮盘横向振动模式. 此处带隙为驻波模态型的Bragg带隙.
第四带隙处于4771.9 ~ 5386.2 Hz, 该带隙内包含群速度为0的色散曲线, 整个结构处于局域共振模式. 由上文分析可知: 系统处于叶片主导的弯曲振动模式, 此时带隙属于局域共振带隙. 处于该频率范围的激励源作用于单胞时, 能量将汇集到叶片上.
下面以机械能流示意图和振动位移图(如图10和图11所示)来验证能带结构对弯曲波的调制能力.
为保证在上述响应图所示时刻, 周向弯曲波具有相同的时间相位, 均取第9.6个周期时刻, 此时激振力幅值相同. 激励探针设置在轮盘边缘激励位置以左的轮缘上某点.
结果表明: 当激励载荷的频率处于带隙外时, 除激励点外, 机械能流密度在整个结构内都处于均匀且极小的分布状态. 该时刻的位移图可以直观显示叶轮内部的位移响应分布: 位移响应在轮盘上的分布均匀, 叶片虽有一定的不对称性, 但会发现整体呈现周向对称的特征; 在禁带频率内的激励载荷作用于叶轮时, 机械能流密度几乎只集中于叶片部分, 结构进入叶片弯曲振动模式(子结构局域共振模式). 此时叶片布置形式带来的影响更加显著, 可以看到完全相同的左右行波所激发的叶片响应并不相同: 左行波作用下的叶片机械能流密度相对较低, 其分布相对均匀, 基本不存在驻波叶片; 右行波作用下的叶片出现了驻波叶片且整体机械能流密度更大. 位移图表明: 此时局部叶片(激励点左侧第二个叶片)位移响应极大. 对应叶片处的加速度、弹性应变能流密度在整个结构中处于最大值, 损坏风险激增. 该最大加速度存在于叶片顶部, 大小与激励载荷幅值成正比. 同时, 参考动能密度会发现其在叶轮上的分布与机械能流分布一致. 可能发生动能损坏的危险位置主要集中于叶轮一侧的各孤立叶片处.
上述分析表明, 激励载荷的频率处于带隙内时, 该激励对部分叶片的寿命可产生负面影响, 严重时将导致断裂. 比较激励载荷的频率处于带隙内和带隙外时的位移响应能够发现: 振动幅值之比大于2.5. 因此, 在叶轮设计中, 应当充分考虑参数变化时叶轮所具有的能带结构的影响, 而不能仅从谐响应结果和静力学强度考虑几何构型. 原因在于: 谐响应结果中的幅值与激励施加方式相关. 为得到完整的共振危险频率点需要考虑各种可能的力的施加方式, 例如, 图12为叶轮的谐响应结果.
图12直观显示了横向激励和x方向激励的位移响应峰, 两方向的激励下最大响应峰均在3100 Hz左右. 而仅x向激励显示了5400 Hz的响应峰, 其幅值相较横向激励下的第一响应峰相差甚远. 从上述结果可以发现: 在常见的激励方式下, 所识别的响应峰与实际的危险频率存在一定的出入, 甚至有可能出现漏判的情况. 图中所示在3100 Hz时的响应峰对应1节径轴对称振动, 整个结构的平均响应达到最大.
为了考察谐响应法对叶片响应峰的识别效果, 将响应点选取在某叶片, 激励类型选择为横向激励. 图13为谐响应结果.
图13表明, 即使将响应点选在叶片上, 对危险频率识别的准确性仍十分有限. 结合谐响应结果可得结论: 在含子结构的弱耦合系统中, 谐响应法识别危险模态的准确性将受到影响, 可能会忽略危险频率, 对系统的几何结构设计产生不当的指导. 综上所述, 对于含子结构的弱耦合系统, 能带结构法对危险频率识别的准确性更高.
4. 结论
本文结合链式周期结构的能带结构理论, 提出了循环周期对称结构能带结构的分析方法, 并揭示了其振动局部化的机制, 得出以下结论.
(1) 基于梁-立柱力学模型的弯曲波方程和纵波方程, 利用力学平衡、力矩平衡、截面应力连续和切应力互等条件, 建立了波的输入截面及输出截面间的传递矩阵. 同时利用相空间Lyapunov指数概念和周期结构中的Bloch定理, 计算了局部化因子和能带结构. 结果表明, 梁-立柱模型得到的结果与实际情况出入较大, 在类似的动力学分析中, 应舍弃该力学模型.
(2) 基于薄板力学模型的弯曲波方程, 将其通解写为分离变量法得到的周向和径向的乘积形式. 通过对解的拓展, 得到在某一半径邻域内的解形式, 与梁模型相同, 推导该单胞结构的传递矩阵并计算局部化因子. 结果表明, 该方法可以预测范围在5000 Hz左右的弯曲波带隙.
(3) 结合有限元方法, 提出了计算单胞能带结构的数值方法, 得到了色散曲线. 通过极化率可以识别对应曲线上的振动模式, 据此确定了4771.9 ~ 5472 Hz的带隙为弯曲波带隙. 响应结果验证了该叶轮存在此弯曲波带隙, 同时发现叶片的布置方式将导致激励点两侧叶片上的能量分布不均匀. 右行波作用的一侧整体动能密度大于左侧, 相同服役条件下, 该侧更易发生叶片顶部的断裂.
对比本文提出的能带结构法与传统谐响应法, 结果表明: 传统谐响应法在弱耦合系统的危险频率识别准确性上有所欠缺. 同时, 谐响应法要得到相对正确的结果需要设置不同位置的响应点及不同的激励方式; 本文提出的利用有限元计算能带结构法可以快速确定危险频率的范围, 从色散曲线斜率可以判断是否存在局域共振模式, 再利用极化率可以准确识别振动模式从而确定带隙类型. 此结果可以应用于叶轮设计, 准确确定危险频率范围, 从而提前规避部分叶片断裂风险.
-
![]()
图 4 波矢与频率间的关系(前3阶结果)
Figure 4. Relationship between wave vector and frequency (first third order)
![]()
图 8 以$ {60^ \circ } $为圆心角的单胞能带结构
Figure 8. Single cell band structure with a center angle of $ {60^ \circ } $
-
[1] Zhang JZ. Dynamics and Fault Diagnosis of Nonlinear Rotors and Impellers. Cham: Springer Nature Switzerland AG, 2022
[2] 高锋, 马杰. 200MW机组叶片断裂分析与控制//全国火电100-200 MW级机组技术协作会年会论文集. 桂林: 中国电力企业联合会, 2009: 46-51 (Gao Feng, Ma Jie. Analysis and control of blade fracture in a 200MW unit//Proceedings of the 2009 Annual Conference of the National Technical Cooperation Conference for 100-200 MW Thermal Power Units. Guilin: China Electricity Council, 2009: 46-51 (in Chinese) Gao Feng, Ma Jie. Analysis and control of blade fracture in a 200MW unit//Proceedings of the 2009 Annual Conference of the National Technical Cooperation Conference for 100-200 MW Thermal Power Units. Guilin: China Electricity Council, 2009: 46-51 (in Chinese)
[3] 赵问银, 张家忠, 周成武. 大型离心叶轮振动模态局部化特性研究. 应用力学学报, 2012, 29(6): 699-704, 775 (Zhao Wenyin, Zhang Jiazhong, Zhou Chengwu. Study on vibration mode localization characteristics of large centrifugal impeller. Chinese Journal of Applied Mechanics, 2012, 29(6): 699-704, 775 (in Chinese) Zhao Wenyin, Zhang Jiazhong, Zhou Chengwu. Study on vibration mode localization characteristics of large centrifugal impeller. Chinese Journal of Applied Mechanics, 2012, 29(6): 699-704, 775 (in Chinese)
[4] Pal RK, Ruzzene M. Edge waves in plates with resonators: An elastic analogue of the quantum valley hall effect. New Journal of Physics, 2016, 19(2): 025001
[5] Andrea C, Craster R, Daniel C, et al. Elastic wave controlbeyond band-gaps: Shaping the flow of waves in plates and half-spaces with subwavelength resonant rods. Frontiers in Mechanical Engineering, 2017, 3: 10 doi: 10.3389/fmech.2017.00010
[6] Anufriev R, Yanagisawa R, Nomura M. Aluminium nanopillars reduce thermal conductivity of silicon nanobeams. Nanoscale, 2017, 9(39): 15083 doi: 10.1039/C7NR05114J
[7] Chaunsali R, Chen CW, Yang J. Subwavelength and directional control of flexural waves in zone-folding induced topological plates. Physical Review B, 2018, 97(5): 054307 doi: 10.1103/PhysRevB.97.054307
[8] Chen HF, Fu Y, Ling L, et al. Design of locally resonant acoustic metamaterials with specified band gaps using multi-material topology optimization. Materials, 2024, 17(14): 3591 doi: 10.3390/ma17143591
[9] Wang C, Zhao HG, Wang Y, et al. Topology optimization of chiral metamaterials with application to underwater sound insulation. Applied Mathematics and Mechanics, 2024, 45(7): 1119-1138 doi: 10.1007/s10483-024-3162-8
[10] Xin JY, Li YJ, Li DX, et al. Comprehensive analysis of band gap modulation of hexagonal fan blade and optimized ligament structure in the low-frequency range. Micro and Nanostructures, 2024, 193: 207918 doi: 10.1016/j.micrna.2024.207918
[11] Li ZH, Yang SJ, Liu Q, et al. A novel sandwich structure for integrated sound insulation and absorption. International Journal of Mechanical Sciences, 2024, 279: 109526 doi: 10.1016/j.ijmecsci.2024.109526
[12] Alicia G, Roger R, James W, et al. An adjustable acoustic metamaterial cell using a magnetic membrane for tunable resonance. Scientific Reports, 2024, 14(1): 15044 doi: 10.1038/s41598-024-65819-2
[13] Cao W, Shi J, Xiong R, et al. Thermal transport in 2D nanophononic metamaterials embedded with cylindrical arrays. Physics Letters A, 2023, 481: 128997 doi: 10.1016/j.physleta.2023.128997
[14] Sigalas MM, Economou EN. Elastic and acoustic wave band structure. Journal of Sound and Vibration, 1992, 158(2): 377-382 doi: 10.1016/0022-460X(92)90059-7
[15] Kushwaha MS, Halevi P, Dobrzynski L, et al. Acoustic band structure of periodic elastic composites. Physical Review Letters, 1993, 71(13): 2022-2025 doi: 10.1103/PhysRevLett.71.2022
[16] Jin Y, Djafari-Rouhani B, Torrent D. Gradient index phononic crystals and metamaterials. Nanophotonic, 2018, 8(5): 683-701
[17] 温熙森. 声子晶体. 北京: 国防工业出版社, 2009 (Wen Xisen. Phononic Crystals. Beijing: National Defense Industry Press, 2009 (in Chinese) Wen Xisen. Phononic Crystals. Beijing: National Defense Industry Press, 2009 (in Chinese)
[18] 高南沙, 沈礼, 侯宏. 声子晶体减振降噪特性分析研究及应用. 西安: 西北工业大学出版社, 2017 (Gao Nansha, Shen Li, Hou Hong. Analysis and Application of the Vibration and Noise Reduction Properties of Phononic Crystals. Xi’an: Northwestern Polytechnical University Press, 2017 (in Chinese) Gao Nansha, Shen Li, Hou Hong. Analysis and Application of the Vibration and Noise Reduction Properties of Phononic Crystals. Xi’an: Northwestern Polytechnical University Press, 2017 (in Chinese)
[19] Liu ZY, Zhang XX, Mao YW, et al. Locally resonant sonic materials. Science, 2000, 289(5485): 1734-1736 doi: 10.1126/science.289.5485.1734
[20] Wang G, Liu YZ, Wen JH, et al. Formation mechanism of the low-frequency locally resonant band gap in the two-dimensional ternary phononic crystals. Chinese Physics, 2006, 15(2): 407-411 doi: 10.1088/1009-1963/15/2/029
[21] Yu DL, Liu YZ, Qiu J, et al. Experimental and theoretical research on the vibrational gaps in two-dimensional three-component composite thin plates. Chinese Physics Letters, 2005, 22(8): 1958 doi: 10.1088/0256-307X/22/8/038
[22] Hsu JC, Wu TT. Lamb waves in binary locally resonant phononic plates with two-dimensional lattices. Applied Physics Letters, 2007, 90(20): 201904 doi: 10.1063/1.2739369
[23] Assouar MB, Oudich M. Enlargement of a locally resonant sonic band gap by using double-sides stubbed phononic plates. Applied Physics Letters, 2012, 100(12): 123506
[24] Serrano Ó, Zaera R, Fernández-Sáez J. Band structure analysis of a thin plate with periodic arrangements of slender beams. Journal of Sound and Vibration, 2018, 420: 330-345 doi: 10.1016/j.jsv.2017.11.016
[25] 董立强. 一维广义声子晶体圆板振动特性研究. [硕士论文]. 哈尔滨: 哈尔滨工程大学, 2015 (Dong Liqiang. Research on vibrational characteristics of circular plate of generalized phononic crystals. [Master Thesis]. Harbin: Harbin Engineering University, 2015 (in Chinese) Dong Liqiang. Research on vibrational characteristics of circular plate of generalized phononic crystals. [Master Thesis]. Harbin: Harbin Engineering University, 2015 (in Chinese)
[26] 罗金雨, 姚凌云, 江国期等. 一种圆柱壳类声子晶体振动带隙及振动特性研究. 振动与冲击, 2019, 38(8): 133-138 (Luo Jinyu, Yao Lingyun, Jiang Guoqi, et al. A study on the vibration band gap and vibration characteristics of a cylindrical shell phononic crystal. Journal of Vibration and Shock, 2019, 38(8): 133-138 (in Chinese) Luo Jinyu, Yao Lingyun, Jiang Guoqi, et al. A study on the vibration band gap and vibration characteristics of a cylindrical shell phononic crystal. Journal of Vibration and Shock, 2019, 38(8): 133-138 (in Chinese)
[27] 黄洪赛, 冉冀林, 陈凯伦等. 局域共振型圆柱壳类声子晶体带隙特性研究. 人工晶体学报, 2020, 49(6): 1078-1082, 1106 (Huang Hongsai, Ran Jilin, Chen Kailun, et al. Study on band gap of locally resonant cylindrical shell phononic crystals. Journal of Synthetic Crystals, 2020, 49(6): 1078-1082, 1106 (in Chinese) Huang Hongsai, Ran Jilin, Chen Kailun, et al. Study on band gap of locally resonant cylindrical shell phononic crystals. Journal of Synthetic Crystals, 2020, 49(6): 1078-1082, 1106 (in Chinese)
[28] 姚凌云, 姚敦辉. 圆柱壳弹性波超材料分级排列的带隙拓宽方法研究. 动力学与控制学报, 2023, 21(7): 38-42 (Yao Lingyun, Yao Dunhui. On the cell-dependent vibrations and wave propagation in uniperiodic cylindrical shells. Journal of Dynamics and Control, 2023, 21(7): 38-42 (in Chinese) Yao Lingyun, Yao Dunhui. On the cell-dependent vibrations and wave propagation in uniperiodic cylindrical shells. Journal of Dynamics and Control, 2023, 21(7): 38-42 (in Chinese)
[29] 毕红霞, 王艾伦, 曹旭辉. 基于应变模态的自带冠叶盘结构振动局部化问题. 华东理工大学学报, 2012, 38(4): 524-528 (Bi Hongxia, Wang Ailun, Cao Xuhui. Vibration localization of bladed disk with shrouds based on strain model. Journal of East China University of Science and Technology, 2012, 38(4): 524-528 (in Chinese) Bi Hongxia, Wang Ailun, Cao Xuhui. Vibration localization of bladed disk with shrouds based on strain model. Journal of East China University of Science and Technology, 2012, 38(4): 524-528 (in Chinese)
[30] 王芳隆, 沈一舟, 徐艳龙等. 弯曲波彩虹捕获效应及其在能量俘获中的应用. 力学学报, 2022, 54(10): 2695-2707 (Wang Fanglong, Shen Yizhou, Xu Yanlong, et al. Rainbow trapping of flexural waves and its application in energy harvesting. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(10): 2695-2707 (in Chinese) Wang Fanglong, Shen Yizhou, Xu Yanlong, et al. Rainbow trapping of flexural waves and its application in energy harvesting. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(10): 2695-2707 (in Chinese)
[31] 吴锋. 快速计算随机声子晶体的局部化因子//第十一届全国随机振动理论与应用学术会议论文集. 北京, 2015: 1-8 (Wu Feng. Computing efficiently the localization factors for the random phononic crystals//Proceedings of the 11th National Conference on the Theory and Application of Random Vibration. Beijing, 2015: 1-8 (in Chinese) Wu Feng. Computing efficiently the localization factors for the random phononic crystals//Proceedings of the 11th National Conference on the Theory and Application of Random Vibration. Beijing, 2015: 1-8 (in Chinese)
[32] Yao DH, Xiong MK, Luo JY, et al. Flexural wave mitigation in metamaterial cylindrical curved shells with periodic graded arrays of multi-resonator. Mechanical Systems and Signal Processing, 2022, 168: 108721 doi: 10.1016/j.ymssp.2021.108721
[33] Ma TX, Fan QS, Li ZY, et al. Flexural wave energy harvesting by multi-mode elastic metamaterial cavities. Extreme Mechanics Letters, 2020, 41: 101073 doi: 10.1016/j.eml.2020.101073
[34] Mannattil M, Santangelo CD. Geometric localization of waves on thin elastic structures. American Physical Society, 2024, 109(3): 035001
下载:
计量
- 文章访问数: 126
- HTML全文浏览量: 30
- PDF下载量: 35
