Modeling and Dynamics Analysis of High-Temperature Magnetic Bearing-Rotor System

传统线性颤振计算分析理论能较好地预测桥梁断面等钝体结构的颤振临界风速，其理论基础是Scanlan等^[1]所提出的基于小振幅假定的线性自激力模型. 对于目前风洞试验中所发现的非线性颤振（或“软颤振”）现象^[2-4]，由于自激气动力的非线性特性，颤振位移不会以指数型增长发散，而表现出类似于涡激振动的极限环稳态运动^[5]，故以线性颤振导数为基础的传统理论分析方法并不适用，也无法准确计算线性颤振发散临界点后的稳态振幅.

非线性颤振现象是由非线性气动力主导的自限幅运动，类似于美国塔科马大桥风毁事故之前的大振幅运动，即存在较强的非线性效应. Ying等^[6]利用节段模型风洞试验详细研究了流线型箱梁非线性颤振，并提出具体的单自由度非线性颤振计算方法. 朱乐东等^[7-8]研究几种典型断面的软颤振特性，并针对双边肋主梁的非线性颤振现象作出了理论解释. 王骑^[9]基于泰勒展开的非线性自激力模型，研究大跨度桥梁在大振幅运动下的颤振稳定性问题. 伍波等^[10]研究扁平箱梁颤振后状态的振幅依存性问题，指出同一风速下具有多个稳定振幅点的可能性，并随后详细研究了钢桁架梁在不同初始条件下发生的滞回响应动力学特性及演化机理^[11]. 林思源^[2]利用强迫振动测压试验和自由振动风洞试验研究5∶1矩形断面非线性气动力高次谐波特性、气动力叠加性和振幅依存性. Gao等^[12]利用自由振动天平测力试验详细研究了双边肋主梁在非线性颤振状态下的气动力特性，并对该现象作出机理解释. 以往学者的研究增进了钝体断面非线性气动力和非线性颤振问题的理解，较为全面地阐释了钝体断面非线性问题的动力学机理.

综合以往研究可以发现：针对非线性颤振的研究主要基于现有气动力模型的拓展，而对气动力模型的参数识别主要依赖于测力或测压试验. 现有技术条件下，仍然无法摆脱测力或测压试验进行气动力辨识. 尽管Wu等^[13]实现了基于自由振动试验的参数识别，但仍依赖于颤振计算闭合解理论；由于方法的局限性，无法一次性识别出所有待识别参数，且未对非线性气动力进行时域描述. 综合以上气动力辨识问题，所使用的处理方法是确定系统输入到输出的映射函数，此方法的前提条件是已知风致振动系统的模型结构（即自激力模型），且模型是待估计参数的线性函数. 若将上述方法推广到非线性参数识别问题或模型结构未知的非线性系统，不仅处理起来较为繁琐，且容易面临无法收敛到全局最优解的问题. 由此，上述围绕均匀流下的非线性颤振研究并不能拓展到紊流条件下的大振幅运动计算分析，此问题目前仍未有较好的解决办法. 近些年，基于神经网络的深度学习技术逐渐应用到各个领域，并取得了显著成效. 神经网络方法的优点主要在于通过对样本数据的训练即可完成系统输入到输出的映射关系描述，而无需过多关心系统的模型结构. Chu等^[14-15]早在1990年将神经网络理论应用到系统辨识领域中，并取得了较好的效果. Masri等^[16]动力学研究者又进一步将神经网络方法应用到动力学系统辨识当中，但其仅限于系统非线性识别，而无法对系统外荷载进行辨识.

基于以上讨论，本文以神经网络方法和运动方程数值求解方法为基础，提出一种基于Encoder-decoder的自编码器模型. 该模型着眼于解决现有气动力辨识方法的不足，旨在摆脱现有的测力或测压技术，仅依靠自由振动测振试验进行非线性气动力辨识，为钝体断面非线性颤振计算提供参考. 选取5∶1矩形断面为研究对象，首先利用无风条件下的模型自由衰减运动，获取弹簧悬挂系统结构非线性阻尼的多项式表达，并获取该断面在两组动力参数下不同折算风速的唯一稳态振幅. 利用同步测压测振技术，获取各风速下运动时程数据和压力时程数据. 基于风洞试验数据和所提出的自编码器模型，实现了5∶1矩形断面在非线性颤振状态下的非线性自激力辨识和非线性颤振运动时程预测.

1.   节段模型风洞试验

1.1   试验概况

本次试验在西南交通大学XNJD-2风洞开展，试验流场为均匀流场，采用传统的弹簧悬挂系统进行节段模型颤振试验（图1），并在模型中安装DSM4000电子压力扫描阀模块（图2）. 试验支架、弹簧、激光位移计等仪器均置于风洞外侧，TFI眼镜蛇风速测量探头安装在模型前端0.3 m处，以保证不干扰模型周围风场. 利用信号同步触发装置，实现脉动风压与位移的同步测量，各套采集系统的采样频率均设置为256 Hz.

图  1  悬挂在风洞中的模型

Figure  1.  Model in wind tunnel

下载: 全尺寸图片 幻灯片

图  2  测压模型内部扫描阀

Figure  2.  Scan valves in model for pressure test

下载: 全尺寸图片 幻灯片

由于5∶1矩形断面非线性颤振现象明显，且测压孔布置更为便捷，因此，以5∶1矩形断面为研究对象. 模型长度L=1.1 m，宽度B=0.3 m，高度H=0.06 m. 沿模型长度方向分别布置5排（A～E）表面脉动压力测点，每排包含一个完整测试条带，共计50个测压点，测点布置详情如图3所示. 图中，小红点为测压点. 模型两端安装有端板，以保证模型周围流场的二维流动.

图  3  测点布置

Figure  3.  Layout of measuring points

下载: 全尺寸图片 幻灯片

无风条件下，通过大振幅激励的方式，获取大振幅运动下动力系统的频率和阻尼参数. 表1列出了2种动力参数组合（D1和D2）. 表中：$ m $为质量，$ I $为质量惯性矩，$ f_{{\mathrm{h}}} $和$ f_{{\text{α}}} $分别为竖向和扭转频率. 通过试验发现：不同振幅下系统频率变化较小，扭转振幅为0.5°～6.0°，扭转频率变化为0.8%，竖向振幅为1～20 mm，竖向频率变化为0.7%，因此，可不考虑系统刚度非线性对系统频率结果的影响.

表  1  节段模型试验动力参数

Table  1.  Dynamic parameters of sectional model test

动力参数 类型	$ m $/kg	$ I $/（kg•m2•m−1）	$ f_{{\mathrm{h}}} $/Hz	$ f_{\text{α}}$/Hz
D1	18.56	0.4928	1.768	3.038
D2	15.36	0.2152	1.875	4.435

下载: 导出CSV 

| 显示表格

1.2   动力响应测试

试验测试了2种动力参数下，两自由度耦合系统（2DOF）在不同风速下的稳态振幅响应（攻角为0°）. 参照文献[11]，非线性颤振稳态响应特征依赖于系统初始运动状态，不同初始激励条件会导致稳态幅值的变化，故本次试验采取2种不同振幅水平的初始激励对系统进行激振，以获取不同初始条件下系统运动响应发展时程. 经检验：无论给定大振幅还是小振幅初始激励（$ A_{{{\text{α}}}, 0} $），两者均衰减到基本相同的限环运动响应稳态极幅值（$ A_{{{{\text{α}}, {\mathrm{st}}}}} $），故该系统不存在由于初始运动状态的不同而导致滞回响应现象. 图4给出了2种不同动力参数系统的典型运动响应发展时间历程. 据此，获取了不同折算风速$ V=U /\left(f_{{\text{α}} } B\right) $条件下2DOF系统的唯一稳态振幅，如图5所示. 其中：U为来流风速.

图  4  不同初始激励下的运动响应时程

Figure  4.  Motion response time-history under different initial excitations

下载: 全尺寸图片 幻灯片

图  5  稳态扭转响应振幅随折算风速的变化

Figure  5.  Variation of steady-state torsional response amplitude with reduced wind speed

下载: 全尺寸图片 幻灯片

由图4可知，2DOF系统竖向运动参与较少，耦合响应基本由扭转分支主导，不同折算风速下的竖向和扭转响应无量纲振幅比$ 180 h /(B {\text{π}} ) $约为0.031～0.033，其中，$ h $和$ \alpha $分别为竖向位移时程和扭转位移时程. 故可以判断：不同风速下的非线性颤振主要由扭转自激力主导，2DOF系统可简化描述为仅考虑扭转自由度的单自由度（SDOF）运动系统，因此，可忽略竖向运动对耦合系统的影响，动力方程可简化描述为

式中：$ \omega_{a} $为扭转圆频率，$ \xi_{\alpha}(\alpha) $为与扭转位移相关的非线性结构阻尼，$ M_{{\mathrm{se}}}(\alpha, \dot{\alpha}) $为自激力矩.

1.3   非线性结构阻尼

结构阻尼的非线性是影响非线性颤振稳态振幅大小的关键因素，精确地获取非线性结构阻尼是辨识非线性气动力和非线性颤振响应的前提，故在此对非线性结构阻尼的提取方法进行详细地说明. 结合上文可知，5∶1矩形断面非线性颤振形态由扭转分支主导，竖向运动参与较少，因此，仅对扭转模态阻尼的提取进行陈述. 结构非线性阻尼可表示成扭转运动振幅的函数：

式中：$ A_{\alpha} $为扭转运动瞬时振幅，$ \xi_{A_{n}}$(·)为等效非线性结构阻尼，$ B_{a i}$(i=0～3) 为待拟合系数.

实际上，结构非线性阻尼$ \xi_{\alpha}(\alpha) $可表示成与位移有关的多项式形式^[17]，如式（3）.

式中：$ b_{\alpha k} $(k=0～3)为待拟合系数.

利用谐波平衡法（式4）将式（3）进行等效线性化，可获得2种表达形式的转换关系如式（5）.

利用希尔伯特变换（Hilbert transform, HT）获得衰减位移的$ A_{\alpha} $和相位$ \vartheta(t) $，并根据式（6）拟合瞬时振幅，结果如图6所示. 根据频率和阻尼比的定义式（式7）和式（8）可获得时变频率$ f_{\alpha}(t) $和阻尼比$ \xi_{\alpha}(t) $. 经检验，利用HT获得的时变频率与表1所示的结构频率吻合. 与此同时，获取了2种动力参数下的扭转模态阻尼为$ \xi_{A_{\alpha}, {\mathrm{D}} 1} $和$ \xi_{A_{\alpha}, {\mathrm{D}} 2} $（图7）. 利用上述多项式系数转换关系，可得到2种动力参数下扭转模态非线性结构阻尼随$ \alpha $变化的表达式（式（9））.

图  6  自由振动衰减时程及其包络线

Figure  6.  Free vibration decaying time-history and its envelope

下载: 全尺寸图片 幻灯片

式中：c₁～c₄均为待拟合系数.

图  7  非线性结构阻尼随扭转振幅的变化

Figure  7.  Variation of nonlinear structural damping with torsional amplitude

下载: 全尺寸图片 幻灯片

2.   气动力辨识方法

2.1   自编码器模型

众多学者提出了适用于非线性条件下的自激力模型，如基于泰勒展开的多项式模型^[18]，颤振导数随振幅变化的等效线性模型等^[13]. 利用以上模型能有效地计算非线性振幅，难点均在于参数识别，高精度的参数识别是有效提取自激气动力的关键.

神经网络方法高精度描述动力系统内部未知隐藏状态，而无需进行显式表达. 其本质是由众多基本神经元组成的复杂网络结构，具体形式可根据问题复杂度灵活选取. 然而，目前诸多关于神经网络的研究均着眼于描述系统输入、输出映射^[19-21]，其网络参数的构成严重依赖训练样本数据. 若仅从输入、输出关系对气动力进行辨识，需要先利用测力试验或数值模拟获取气动力，再利用神经网络建立映射关系（输入与输出分别为位移和气动力），但是利用神经网络的拟合能力，并不能获取系统内部的隐藏状态（气动力）.

本文提出一种基于神经网络的自编码器模型，如图8所示. 该模型由2部分组成：编码器（Encoder）和解码器（Decoder）. Encoder为神经网络模型，用于输出结构第$ i + 1 $个时间步的系统内部状态（非线性自激力矩），可表示成$ \hat{M}_{{\mathrm{se}}, i + 1}=E\left(\alpha_{i}, \dot{\alpha}_{i}\right) $，E（·）代表Encoder，即气动力由第$ i $个时间步的位移$ \alpha_{i} $和速度$ \dot{\alpha}_{i} $表达. 限于篇幅，本文不再对神经网络理论进行陈述，相关知识可参考文献[22]；Decoder基于$ \hat{M}_{{\mathrm{se}}, i + 1} $计算第$ i + 1 $个时间步的位移$ \hat{\alpha}_{i + 1} $、速度$ \hat{\dot{\alpha}}_{i + 1} $和加速度$ \hat{\ddot{\alpha}}_{i + 1} $，可由常用的数值方法求解动力方程（式（1））获取，如Newmark-β法，故Decoder可表示成$ \hat{\alpha}_{i + 1}, \hat{\dot{\alpha}}_{i + 1}, \hat{\ddot{\alpha}}_{i + 1} = D\left(\hat{M}_{{\mathrm{se}}, i + 1}, {\alpha}_{i}, \dot{\alpha}_{i}, \ddot{\alpha}_{i}, \xi_{\alpha}(\alpha)\right) $，D(·)代表Decoder. 进一步地，将$ \alpha_{i + 1} $和$ \dot{\alpha}_{i + 1} $输入到Encoder和Decoder即可计算下一个时间步的位移、速度和加速度. 构建如式（10）所示的损失函数，并利用梯度下降法进行损失最小化，即可获取Decoder中神经元的权重$ {w} $和偏置$ b $. 值得注意的是：该自编码器模型中，仅Encoder包含待识别的未知参数（$ {w} $和$ b $），且仅与当前折算风速有关；Decoder中各项参数均为定值，仅与动力系统参数有关，无待识别参数.

图  8  自编码器模型

Figure  8.  Autoencoder model

下载: 全尺寸图片 幻灯片

式中：$ s $为预测步长，可灵活设定；$ p $为速度损失折减系数，以保证位移损失和速度损失量级相当.

由以上描述可知：自编码器模型的输入和输出均为与位移相关的状态，而无需输入隐藏在模型内部的气动力状态. 整套识别系统可将位移状态（测量状态）和气动力状态（隐藏状态）分开，由Encoder输出隐藏状态，Decoder输出测量状态. 因此，该模型可以同时实现气动力的辨识和位移的预测.

2.2   模型训练

首先利用D1所获取的某一折算风速下的响应时程进行模型训练. 考虑到均匀流下的非线性气动力构成相对简单，构建如图8所示的包含4层神经网络的Encoder模型，该模型包含2个输入神经元（$ \alpha_i, \dot{\alpha}_i $），2层分别包含有4个隐藏神经元的隐藏层以及有一个输出神经元的输出层. 激活函数选用$ \operatorname{Re} {\mathrm{LU}}=\max \{0, x\} $，其中，$ {x} $为神经元的输入参数，学习率$ \eta $大小（即参数更新的增量大小）采用随训练轮次逐步衰减的形式，$ \eta=\eta_{0} c^{{e}} $，其中$ \eta_{0}=0.001 $，为初始学习率，$ c=0.98 $，e为训练轮次.

该自编码器模型训练可由以下步骤实现：

步骤1　给定$ s, m, \alpha_{1}, \dot{\alpha}_{1}, \ddot{\alpha}_{1} $，给定训练样本数N，给定总训练轮次E₁；

步骤2　初始化Encoder网络参数；

步骤3　对于第一个训练样本，输入$ \alpha_{1}{\text{、}} \dot{\alpha}_{1} $至Encoder计算$ \hat{M}_{{\mathrm{se}}, 2} $；

步骤4　输入$ \hat{M}_{{\mathrm{se}}, 2} $至Decoder （Newmark-β法）计算下一个时间步$ \hat{\alpha}_{2}{\text{、}}\hat{\dot{\alpha}}_{2}{\text{、}} \hat{\ddot{\alpha}}_{2} $；

步骤5　利用步骤 4计算的位移与速度预测值重复步骤 3和步骤4 直至完成包含$ s $步长（即一个训练样本）的位移与速度时序预测；

步骤6　计算长度为$ s $的位移与速度时序损失，以损失函数梯度方向作为参数更新的增量方向，更新Encoder中各项待识别参数，并作为下一次训练的网络参数初值；

步骤7　以预测值作为下一轮训练的输入，重复步骤 3～ 6，直至完成训练样本个数为N的网络训练，此为一轮训练；

步骤8　反复利用训练样本进行Encoder训练，直至完成总训练轮次E₁或$ J<1 \times 10^{-8} $；

训练完毕的自编码器模型已对各项待训练参数进行全局最优搜索识别，系统隐藏气动力状态可由Encoder部分输出，其本质是利用神经网络中的各项参数对气动力进行描述. 对于Decoder部分，其参数仅依赖系统动力参数，而与气动力各项参数无关，因此，可拓展到任意动力参数组合下的非线性气动力辨识和非线性颤振计算当中.

3.   识别结果

基于上述步骤，本文方法分训练和验证进行. 训练过程：首先利用D1测量的非线性颤振时程数据训练模型，获取特定折算风速下的Encoder；验证过程：将训练完毕的Encoder与D2组成的Decoder进行组合，同时获取结构非线性气动力和响应时程，并与试验结果进行对比，由此可验证训练完成的Encoder中各项识别的参数是否能完备描述不同动力参数系统中非线性气动力与速度和位移的映射关系. 需要注意的是：本文的训练和验证方法不同于传统的神经网络模型训练，训练和验证不体现在数据集的划分（如训练集和测试集等），而体现在同一节段模型在更换动力参数后所训练的Encoder能否有效输出气动力时程并预测颤振响应.

本次模型训练基于Pytorch深度学习框架，优化器选用加入动量项的Adam^[23]，预测步长$ s=20 $，速度损失折减系数$ p=0.01 $. 设定Newmark-β法中参数$ \gamma=0.5, \beta=0.25 $。总训练轮次E₁=200，每一折算风速下的训练样本数为19980个时间步的位移数据和微分得到的速度数据. 训练过程中，利用式（9）实时更新Decoder中的非线性结构阻尼$ \xi_{\alpha}(\alpha) $.

3.1   训练结果

模型训练基于D1下的非线性颤振响应时程数据. 构建不同折算风速下的自编码器模型，并基于Encoder模型进行气动力输出，由此利用Decoder反算出D1动力参数模型在不同折算风速下的稳态响应时程. 给定2种不同的初始条件（$ A_{{\text{α}}, 0}>A_{{\text{α}}, {\mathrm{st}}} $和$ A_{{\text{α}}, 0}<A_{{\text{α}}, {\mathrm{st}}} $）以模拟图4所示的2种颤振发展状态. 需同时检验位移和气动力辨识结果以核验自编码器模型的有效性，风洞试验所测得的气动力可由测压系统所获取的表面压力积分获取.

由同步测得的压力时程数据，并沿测压条带积分，可得到自激力矩为

式中： $ p_{j} $（t）为单个测点压力值，$ l_{j} $为第$ j $个测点的积分长度，$ \theta_{j} $为第$ j $个测点的内法线与模型水平面的夹角，$ X_{j} $与$ Y_{j} $为第$ j $个测点与模型轴心的横向距离与垂直距离.

由于均匀流下自激力沿跨向相关性较强^[2]，此处忽略跨向相关性的影响，并取靠近模型中心（条带C）断面的数据作为研究对象.

图9给出了V=8.24下本文提出的自编码器模型的预测值和试验值对比：自编码器模型可以很好地预测位移时程，其幅值和相位均能较好地拟合. 图10给出了气动力的预测值和试验值对比结果：两者在幅值上基本保持一致，且气动力预测值更为平滑和稳定；结合频域结果可知，风洞试验测压结果包含一定量的噪声，两者基频成分在频率、幅值和相位角上基本一致，高次谐波成分幅值和相位角有一定的差异，且所预测的气动力噪声成分更少. 图10中：$f_0 $为基频. 整体来说，本文方法不仅可以很好地辨识非线性气动力，还具有一定的滤波和降噪功能^[24]，其原因在于Encoder模型仅提取了与非线性颤振相关的气动力成分，过滤了无关的噪声成分和对非线性颤振不做功的偶数次谐波成分^[12,24]，这一特点在后文中将进一步证实.

图  9  训练预测值和试验值对比 （$ V=8.24 $）

Figure  9.  Comparison of prediction results in training and testing results ($ V=8.24 $)

下载: 全尺寸图片 幻灯片

图  10  气动力训练预测值和试验值对比（$ V=8.24 $）

Figure  10.  Comparison of prediction results in aerodynamic force training and testing results ($ V=8.24 $)

下载: 全尺寸图片 幻灯片

不同折算风速下的唯一稳态振幅响应如图11所示. 相同折算风速下相对误差最大为4.83% （V=3.17），平均误差为1.15%. 利用训练完成的自编码器模型可以很好地反算出不同折算风速下的稳态振幅.

图  11  不同折算风速下非线性颤振响应 （D1）

Figure  11.  Nonlinear flutter responses at different reduced wind speeds （D1）

下载: 全尺寸图片 幻灯片

3.2   验证结果

尽管上述训练过程已对模型的有效性和精度进行了严格的反算证明，但仍需验证的问题是：如何保证在相同折算风速下，Encoder模型能有效输出不同Decoder所需的非线性气动力，即训练完毕的Encoder模型是否适用于其他动力参数系统，由此进一步验证Encoder模型是否能准确捕获非线性气动力与位移和速度的关系. 值得说明的是：训练完成的Encoder模型需保证在相同外荷载激励的条件下进行Decoder验证，即需要保证风环境（均匀流）和风攻角相同. 基于以上条件，选取仅改变动力参数的D2系统非线性颤振响应和非线性气动力进行验证.

利用上节训练完成的Encoder模型，仅改变Decoder中与动力系统相关的各项参数，并给定初始条件（$ A_{{\text{α}}, 0}>A_{{\text{α}}, {\mathrm{st}}} $和$ A_{{\text{α}}, 0} < A_{{\text{α}}, {\mathrm{st}}} $），即可获取相同折算风速下的非线性颤振响应时程. 图12给出了V=8.24下D2系统的位移时程预测结果，图13展示了基于Encoder模型的非线性气动力时域和频域辨识结果. 显然，Encoder模型能较精确地描述非线性自激力矩时程，其幅值和相位与风洞试验结果保持一致；由Decoder模型所输出的位移时程均能较好地模拟风洞试验结果，进一步验证了上节所训练模型的有效性. 与此同时，从气动力频域结果可知，Encoder模型能很好地捕捉非线性气动力奇数项高次谐波分量，基频成分的频率、幅值和相位与试验结果保持一致，其余倍频分量在幅值上有一定的差异.

图  12  非线性颤振运动时程验证预测值和试验值对比（$ V=8.24 $）

Figure  12.  Comparison of prediction results and testing results of motion time-history of nonlinear flutter ($ V=8.24 $)

下载: 全尺寸图片 幻灯片

图  13  气动力验证预测值和试验值对比（$ V=8.24 $）

Figure  13.  Comparison of prediction results and testing results of aerodynamic force ($ V=8.24 $)

下载: 全尺寸图片 幻灯片

图14给出了与D1系统相同折算风速下的稳态振幅响应，其中V=8.24时扭转振幅预测误差为0.26%，其他折算风速时无法与风洞结果进行对应，故只能根据图示结果进行定性比较. 该预测结果与风洞试验结果具有较好的一致性，验证了本文所提出模型的有效性.

图  14  不同折算风速下非线性颤振响应 （D2）

Figure  14.  Nonlinear flutter responses at different reduced wind speeds （D2）

下载: 全尺寸图片 幻灯片

4.   结　论

1） 基于自由振动风洞试验获取了5∶1矩形在不同风速下的非线性颤振响应，明确了不同初始激励下颤振振幅的唯一性；

2） 弹簧悬挂系统结构阻尼存在振幅依存性，非线性结构阻尼可由振幅或时变位移进行3次多项式表达；

3） 所提出的自编码器模型仅依赖位移响应时程即可较好地辨识非线性气动力和开展非线性颤振分析，不同动力参数下的预测结果进一步验证了所提出方法的有效性；

4） 所提出的方法通用性较强，是解决紊流条件下非线性颤振分析的潜在途径，后续研究将会对此进行论证.