A method for modeling flexible multi-body dynamics for all-movable rudder system
-
摘要: 全动舵系统作为航天飞行器控制飞行姿态、调整飞行方向的部件，其动态特性对飞行器的正常工作起重要作用。为了开展带有电机伺服系统和舵轴间隙的全动舵系统动力学特性分析，提出基于柔性多体动力学方法的全动舵系统建模方法：采用Craig-Bampton方法建立典型舵面刚柔耦合降阶模型，采用多体动力学方法建立电动舵机连接机构与全动舵面连接，采用偶极子网格法建立基于模态的广义气动力模型。仿真结果表明，自建的模型预测颤振速度为1270 m/s，与商用软件预测对比的偏差小于2%，验证了该建模方法的正确、可行。研究表明，伺服系统的存在会令典型舵面响应存在较大跳跃现象，而舵轴间隙的存在则极大降低了舵面产生极限环振荡的临界速度。Abstract: For the all-movable rudder system as a component to control the flight attitude and adjust the flight direction, its dynamic characteristics play an important role for the normal operation of a spacecraft. In order to analyze the dynamic characteristics of the all-movable rudder with motor servo system and rudder shaft gap, a method for modeling all-movable rudder system based on flexible multi-body dynamics was proposed. The rigid-elastic coupling reduced-order model for typical rudder surface was established via Craig-Bampton method. The connection mechanism between the motor and the all-movable rudder surface was created by multi-body dynamics method. The modal-based generalized aerodynamic model was produced using doublet lattice method. The simulation results show that the predicted flutter speed by the proposed model is 1270 m/s, with a deviation of less than 2% compared with a commercial software, verifying the correctness and feasibility of the modeling method. The research shows that the presence of servo system would cause a large jump in the response of the typical rudder surface, while the presence of the rudder shaft gap would greatly reduce the critical speed of the rudder surface to generate limit cycle oscillations.
-
0. 引言
全动舵系统是由舵机、传动机构与气动力作用下的弹性舵面共同形成的耦合系统，存在着流、固、电、磁等多物理场的耦合作用，可能在一定的条件下引发颤振失稳[1-3]；而一旦发生颤振失稳，将会对航天飞行器造成非常严重的不良后果：结构的振动响应会被不断放大，甚至会造成结构破坏、飞行器姿态失控，乃至飞行失利。
目前国内外的相关研究主要集中于气动力对于舵面颤振特性的影响方面，常采用计算流体动力学/计算结构动力学（CFD/CSD）方法开展全动舵面气动弹性响应计算。针对舵轴间隙对颤振特性的影响[4-7]，学者们多基于虚拟质量法把含间隙的非线性系统划分为3个线性子系统，采用虚拟质量模态作为统一的坐标去表示各子系统。黄程德等[8]发展了基于CFD 和虚拟质量法耦合的间隙非线性气动弹性分析方法，发现全动舵面极限环幅值随着间隙角的增大而增大，且跨声速极限环的临界速度比线性颤振的跨声速凹坑下降20%。但基于虚拟质量法建模仅能考察舵轴间隙的影响，难以拓展到舵机等其他耦合系统。针对舵机对颤振特性影响的研究主要采用假设模态法，即将舵机系统简化为具有复刚度的弹簧系统[9-10]。Shin等[11]对带有间隙环节和电动舵机的舵面进行了气动弹性分析，结果表明：考虑舵机动刚度与不考虑时计算得到的颤振速度和颤振频率均存在差别，而且非线性的影响会改变线性颤振特性，使得舵面出现极限环振荡现象。丁伟涛等[12]基于典型舵面二阶局部刚化模态和弹性模态建立了考虑舵机复刚度的舵系统动力学模型，研究热环境、电动舵机设计参数以及指令信号幅值对颤振速度的影响，并提出在舵机电流环加入超前滞后环节的颤振抑制措施。但将舵机简化为复刚度弹簧时，其机电耦合特性、传动特性等需进行线性化处理，会降低仿真模型与真实模型的匹配程度。
为了能够较为简单且有效地耦合舵机模型、舵轴间隙等因素的影响，本文采用刚柔耦合多体动力学方法建立典型舵面动力学模型，并采用地面颤振试验的集中式气动力模拟方法开展气动弹性耦合仿真，分析舵机模型、舵轴间隙等因素对动力学响应的影响。
1. 全动舵系统建模方法
1.1 刚柔耦合多体动力学模型
多体动力学的研究对象是由多个具有运动学约束、存在大范围相对运动的刚性或柔性部件组成的复杂系统，主要研究内容是这类系统的动力学建模、计算和控制。对于构件在运动过程中同时存在弹性变形的问题，基于Craig-Bampton方法[13]，采用固连在结构体上的体坐标系运动来描述结构体空间大范围运动，相对于体坐标系来描述柔性变形引起的运动，从而建立能够考虑弹性变形与非线性大范围运动相互耦合影响的柔性多体动力学模型。刚柔耦合多体动力学模型为基于模态截断法获得的降阶有限元模型，假设有限元模型为
$$ {\boldsymbol{M}}{{\boldsymbol{\ddot u}}_{\text d}} + {\boldsymbol{K}}{{\boldsymbol{u}}_{\text d}} = {\boldsymbol{F}} \text{，} $$ (1) 其中：M为有限元模型质量矩阵；K为有限元模型刚度矩阵；üd和ud分别为加速度矢量和位移矢量；F为载荷矢量。选取模型受外载荷点和存在位移约束点为边界点，将式中M、K的各行列重新整理，得到刚柔耦合动力学模型
$$ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{M}}_{{\text{bb}}}}}&{{{\boldsymbol{M}}_{{\text{bi}}}}} \\ {{{\boldsymbol{M}}_{{\text{ib}}}}}&{{{\boldsymbol{M}}_{{\text{ii}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\ddot u}}}_{\text{b}}}} \\ {{{{\boldsymbol{\ddot u}}}_{\text{i}}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{K}}_{{\text{bb}}}}}&{{{\boldsymbol{K}}_{{\text{bi}}}}} \\ {{{\boldsymbol{K}}_{{\text{ib}}}}}&{{{\boldsymbol{K}}_{{\text{ii}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{u}}_{\text{b}}}} \\ {{{\boldsymbol{u}}_{\text{i}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{F}}_{\text{b}}}} \\ {\boldsymbol{0}} \end{array}} \right] \text{，} $$ (2) 其中，下标b代表边界点，i代表内部节点。通过Craig-Bampton方法[13]，保持边界点6自由度刚体位移，其余节点位移采用模态化处理，则刚柔耦合动力学模型的位移矢量转化为
$$ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{u}}_{\text{b}}}} \\ {{{\boldsymbol{u}}_{\text{i}}}} \end{array}} \right] = {\boldsymbol{H}}\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{u}}_{\text{b}}}} \\ {\boldsymbol{\eta }} \end{array}} \right] \text{，} $$ (3) 其中：η为边界点固定时的模态自由度位移；H为Craig-Bampton转化矩阵，
$$ {\boldsymbol{H}} = \left[ {\begin{array}{*{20}{c}} {\boldsymbol{\varXi }}&{\boldsymbol{\varPhi }} \end{array}} \right] \text{，} $$ (4) 其中：Ξ为边界点约束模态，由有限元模型刚体模态给定；Φ为固定边界点模态阵型，由有限元固定边界点后的弹性模态给定。
假设动力学分析频率最高为fmax，在H矩阵中对Φ阵型展开模态截断，仅保留模态频率低于fmax的模态，则刚柔耦合动力学模型的位移矢量转化为
$$ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{u}}_{\text{b}}}} \\ {{{\boldsymbol{u}}_{\text{i}}}} \end{array}} \right] \simeq {{\boldsymbol{H}}^ * }\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{u}}_{\text{b}}}} \\ {{{\boldsymbol{\eta }}^ * }} \end{array}} \right] \text{，} $$ (5) 其中：
${\boldsymbol{\eta}}^ * $ 为降阶模态自由度位移；H*为降阶Craig-Bampton转化矩阵。最终降阶子结构刚柔耦合动力学模型为$$ \left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\hat M}}}_{{\text{bb}}}}}&{{{{\boldsymbol{\hat M}}}_{{\text{bm}}}}} \\ {{{{\boldsymbol{\hat M}}}_{{\text{mb}}}}}&{{{{\boldsymbol{\hat M}}}_{{\text{mm}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\ddot u}}}_{\text{b}}}} \\ {{{{\boldsymbol{\ddot \eta }}}^*}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\hat K}}}_{{\text{bb}}}}}&{{{{\boldsymbol{\hat K}}}_{{\text{bm}}}}} \\ {{{{\boldsymbol{\hat K}}}_{{\text{mb}}}}}&{{{{\boldsymbol{\hat K}}}_{{\text{mm}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{u}}_{\text{b}}}} \\ {{{\boldsymbol{\eta }}^{\boldsymbol{*}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{F}}_{\text{b}}}} \\ {\boldsymbol{0}} \end{array}} \right] \text{，} $$ (6) 其中：下标m代表降阶模态自由度，
$ {{\boldsymbol{\hat M}}_{{\text{**}}}} $ 为降阶质量矩阵$ {\boldsymbol{\hat M}} $ 的对应分量，$ {{\boldsymbol{\hat K}}_{{\text{**}}}} $ 为降阶刚度矩阵$ {\boldsymbol{\hat K}} $ 的对应分量。根据$$\left\{ \begin{gathered} \left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\hat M}}}_{{\text{bb}}}}}&{{{{\boldsymbol{\hat M}}}_{{\text{bm}}}}} \\ {{{{\boldsymbol{\hat M}}}_{{\text{mb}}}}}&{{{{\boldsymbol{\hat M}}}_{{\text{mm}}}}} \end{array}} \right] = {{\boldsymbol{H}}^ * }^{\text{T}}{\boldsymbol{M}}{{\boldsymbol{H}}^ * } \\ \left[ {\begin{array}{*{20}{c}} {{{{\boldsymbol{\hat K}}}_{{\text{bb}}}}}&{{{{\boldsymbol{\hat K}}}_{{\text{bm}}}}} \\ {{{{\boldsymbol{\hat K}}}_{{\text{mb}}}}}&{{{{\boldsymbol{\hat K}}}_{{\text{mm}}}}} \end{array}} \right] = {{\boldsymbol{H}}^ * }^{\text{T}}{\boldsymbol{K}}{{\boldsymbol{H}}^ * } \\ \end{gathered} \text{，}\right. $$ (7) 可获得降阶质量和降阶刚度的具体数值。
令
$$ {\boldsymbol{x}} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{\eta }}^*}} \\ {{{{\boldsymbol{\dot \eta }}}^*}} \end{array}} \right] ,\;\;{\boldsymbol{u}} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{u}}_{\text{b}}}} \\ {{{{\boldsymbol{\dot u}}}_{\text{b}}}} \\ {{{{\boldsymbol{\ddot u}}}_{\text{b}}}} \end{array}} \right],\;\;{\boldsymbol{y}} ={\text -} {{\boldsymbol{F}}_{\text{b}}} \text{，} $$ (8) 可将子结构动力学方程转化为状态空间方程的形式
$$ \left\{ \begin{gathered} {\boldsymbol{\dot x}} = {\boldsymbol{Ax}} + {\boldsymbol{Bu}} \\ {\boldsymbol{y}} = {\boldsymbol{Cx}} + {\boldsymbol{Du}} \\ \end{gathered} \right.{\kern 1pt} \text{，} $$ (9) 其中：
$$ \begin{gathered} {\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} 0&{\boldsymbol{I}} \\ { {\text{-}} {{{\boldsymbol{\hat M}}}_{{\text{mm}}}}^{ {\text{-}} 1}{{{\boldsymbol{\hat K}}}_{{\text{mm}}}}}&0 \end{array}} \right] \text{；} \\ {\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} 0&0&0 \\ { {\text{-}} {{{\boldsymbol{\hat M}}}_{{\text{mm}}}}^{ {\text{-}} 1}{{{\boldsymbol{\hat K}}}_{{\text{mb}}}}}&0&{ {\text{-}} {{{\boldsymbol{\hat M}}}_{{\text{mm}}}}^{ {\text{-}} 1}{{{\boldsymbol{\hat M}}}_{{\text{mb}}}}} \end{array}} \right] \text{；} \\ {\boldsymbol{C}} = \left[ {\begin{array}{*{20}{c}} { {\text{-}} ( {{{{\boldsymbol{\hat K}}}_{{\text{bm}}}} - {{{\boldsymbol{\hat M}}}_{{\text{bm}}}}{{{\boldsymbol{\hat M}}}_{{\text{mm}}}}^{ {\text{-}} 1}{{{\boldsymbol{\hat K}}}_{{\text{mm}}}}} )}&0 \end{array}} \right] \text{；} \\ {\boldsymbol{D}} = \left[ {\begin{array}{*{20}{c}} { {\text{-}} ( {{{{\boldsymbol{\hat K}}}_{{\text{bb}}}} - {{{\boldsymbol{\hat M}}}_{{\text{bm}}}}{{{\boldsymbol{\hat M}}}_{{\text{mm}}}}^{ {\text{-}} 1}{{{\boldsymbol{\hat K}}}_{{\text{mb}}}}} )} & 0 & { {\text{-}} ( {{{{\boldsymbol{\hat M}}}_{{\text{bb}}}} - {{{\boldsymbol{\hat M}}}_{{\text{bm}}}}{{{\boldsymbol{\hat M}}}_{{\text{mm}}}}^{ {\text{-}} 1}{{{\boldsymbol{\hat M}}}_{{\text{mb}}}}} )} \end{array}} \right] \text{。} \\ \end{gathered} $$ 1.2 舵轴间隙碰撞模型
建立多个子结构模型状态方程，以顶杆子结构X和舵轴拨片子结构Y为例，它们相互耦合时的位移及受力关系为
$$ \left\{ \begin{gathered} {u_X} - {u_Y} = \sigma \\ {f_{Xu}} = g\left( \sigma \right) \\ {f_{Yu}} = {\text{-}} g\left( \sigma \right) \\ {{\dot u}_X} - {{\dot u}_Y} = \dot \sigma \\ {f_{X\ddot u}} = g\left( {\dot \sigma } \right) \\ {f_{Y\ddot u}} = {\text{-}} g\left( {\dot \sigma } \right) \\ \end{gathered} \right. \text{，} $$ (10) 其中：下标u、ü代表由位移、加速度引起的内力；g为非线性映射关系；σ和
$ \dot \sigma $ 为子结构模型的相对位移和相对速度。弹性接触问题采用罚函数模拟[14]。假设顶杆与舵轴的初始间隙为δ，顶杆与舵轴拨片发生碰撞的判断式为
$$ \left\lfloor \sigma \right\rfloor - \delta {\text{＜}} 1 \times {10^{ {\text{-}} 5}} \text{，} $$ (11) 当顶杆与舵轴拨片的相对位移满足该判断条件，则存在碰撞弹性力Fs与碰撞阻尼力Fd：
$$ {F_{\text{s}}} = \left\{ {\begin{split}& {1 \times {{10}^6} \times \left| \sigma \right|,{\text{ }}\left| \sigma \right| {\text{＜}} {\sigma _{\max }}} \\ &{1 \times {{10}^6} \times {\sigma _{\max }},{\text{ }}\left| \sigma \right| {\text{≥}} {\sigma _{\max }}} \end{split}} \right. \text{；} $$ (12) $$ {F_{\text{d}}} = \left\{{\begin{split}& {1 \times {{10}^3} \times \left| {\dot \sigma } \right|,{\text{ }}\left| {\dot \sigma } \right| {\text{＜}} {{\dot \sigma }_{\max }}} \\ &{1 \times {{10}^3} \times {{\dot \sigma }_{\max }},{\text{ }}\left| {\dot \sigma } \right| {\text{≥}} {{\dot \sigma }_{\max }}} \end{split}} \right. \text{。} $$ (13) 其中最大碰撞弹性力约束σmax和最大碰撞阻尼力约束
$ {\dot \sigma _{\max }} $ 为人工定义参数，碰撞弹性力作用于顶杆与舵轴双模型上，且方向与碰撞点相反；而碰撞阻尼力仅作用于发生碰撞时相对速度较大的模型中。根据顶杆−舵轴拨片间隙碰撞模型，可由顶杆−舵轴拨片的相对位移σ和相对速度$ \dot \sigma $ ，获得碰撞弹性力Fs与碰撞阻尼力Fd。1.3 伺服舵机机电耦合模型
采用传递函数形式建立舵机机电模型，舵机输入电压信号V为伺服系统给定信号，负载转速ωB（减速齿轮的主动齿轮转速）与负载功率W的关系式为
$$ W = \frac{{{k_{\text{t}}}}}{R}\left( {V - {k_{\text{v}}}{\omega _{\text{B}}}} \right){\omega _{\text{B}}} \text{，} $$ (14) 其中，扭矩常数kt、反电动势常数kv和电阻R均为给定参数。
舵机系统传动模型采用传递函数形式建模，其中包含减速齿轮、滚珠丝杠和顶杆装置。包含摩擦影响的减速齿轮模型输出为
$$ {\omega _{\text{F}}} = \frac{{{r_{\text{B}}}{\omega _{\text{B}}}}}{{{r_{\text{F}}}}} \text{，} $$ (15) 其中：ωB、ωF分别为主动、从动齿轮转速；rB、rF分别为主动、从动齿轮半径。从动齿轮的输出功率WF为
$$ {W_{\text{F}}} = W - {W_{{\text{loss}}}} \text{，} $$ (16) 其中Wloss为摩擦损耗功率，
$$ {W_{{\text{loss}}}} = {r_{\text{F}}}\left( {1 - \eta } \right)\tanh \left( {4{\omega _{\text{F}}}/{\omega _{\text{B}}}} \right) \text{，} $$ (17) 其中η为摩擦损耗系数，由人工指定。
滚珠丝杠模型输出为
$$ v = \frac{{{\omega _{\text{F}}}L}}{{2{\text{π}}}} \text{，} $$ (18) 其中：v为丝杠的平动速度；L为丝杠长度。假设滚珠丝杠为理想模型，则丝杠的输出功率与从动齿轮输出功率WF相同。
同样假设顶杆装置为理想模型且与丝杠连接处不存在间隙摩擦，则顶杆装置模型输出速度为v，输出功率与丝杠输出功率WF相同，此时顶杆的输出力为
$$ F = \frac{{{W_{\text{F}}}}}{v} 。 $$ (19) 舵机一般采用电流环、速度环、位置环三闭环控制，如图1所示，电流环居内，速度环居中，位置环为最外环。电流环输入为舵机内电流，速度环输入为舵机输出速度，位置环输入为舵轴偏转角度，根据PID控制算法调整相应控制参数，可令舵机伺服特性满足需求。
1.4 集中气动力模型
基于偶极子网格法和最小状态法有理函数近似拟合获得[15]广义气动力模型的一般表达式
$$ {\boldsymbol{\bar Q}}\left( s \right) = {{\boldsymbol{A}}_0} + \frac{L}{V}{{\boldsymbol{A}}_1}s + \frac{{{L^2}}}{{{V^2}}}{{\boldsymbol{A}}_2}{s^2} + {\boldsymbol{D}}{\left( {{\boldsymbol{I}}s - \frac{V}{L}{\boldsymbol{R}}} \right)^{ {\text{-}} 1}}{\boldsymbol{E}}s \text{，} $$ (20) 其中：s为拉普拉斯域变量；A0代表气动刚度；A1代表气动阻尼；A2代表非定常气动力表观质量；D、R和E为气动力响应系数矩阵；I为单位矩阵。根据1.1节多体动力学模型中边界点选取获得边界点前n阶模态阵型为Ψ，其维度为n×N，其中N为舵面所受集中气动力个数。则舵面响应位置上所受集中气动力为
$$ {{\boldsymbol{F}}_{\text{a}}} = {{\boldsymbol{\varPsi }}^{{\text{-T}}}}{\boldsymbol{\bar Q}} 。 $$ (21) 2. 全动舵系统模型
本章验证采用刚柔耦合降阶模型建立典型舵面动力学模型的正确性，首先开展刚柔耦合降阶模型的静载荷激励下响应与商用有限元软件计算结果对比。
典型舵面的几何尺寸如图2所示，据此建立的几何模型如图3所示，其有限元模型如图4所示，有限元单元个数为1375，节点个数为1460。
根据有限元模型，固定节点后求解典型舵面模态振型，保留其前6阶模态，构造典型舵面的刚柔耦合降阶模型。固定该模型舵轴处的自由度，给定舵面翼尖后缘静载荷激励，获得典型舵面翼尖后缘处节点的时域响应，并与商用有限元软件计算结果进行对比。为了提升计算效率并降低数值求解发散的概率，经不同算法对比后，本文基于刚柔耦合降阶模型的动力学数值算法均采用MatLab中的变步长ode23t算法。对比结果如图5所示，基于Craig-Bampton方法建立的典型舵面刚柔耦合降阶模型与基于商用软件建立的有限元模型之动力学响应间的偏差较小，证明了采用基于Craig-Bampton方法建立典型舵面降阶模型的正确性。
![]()
图 5 无阻尼典型舵面静载荷下时域响应的计算结果对比Figure 5. Comparison of calculation results for the time domain response of a typical rudder surface under static load with no damping典型舵面边界点取6个，以翼根前缘为坐标原点，6个节点的坐标分别为(0, 0, 0)(550, 0, 0)(200, 250, 0)(550, 250, 0)(280, -80, 0)(280, -40, 20)，其中：节点1～4为激振器激励节点，连接方式采用多体动力学连接模型；节点5为舵轴拨片节点，与顶杆连接，开展无间隙仿真时连接方式采用线性多体动力学连接模型，开展间隙仿真时采用顶杆−舵轴拨片间隙碰撞模型；节点6为舵轴根部节点，与支撑结构相连，连接方式采用多体动力学连接模型。
3. 数值算例
3.1 全动舵面颤振结果分析
建立典型舵面计及阻尼的刚柔耦合降阶模型，其阻尼由如下比例给定，
$$ {{\boldsymbol{C}}_i} = \alpha {{\boldsymbol{M}}_i} + \beta {{\boldsymbol{K}}_i}\text{，} $$ (22) 其中，α取值为0，β取值为0.001。计及阻尼后，典型舵面基于偶极子网格法计算获得的颤振速度为1256 m/s，基于本文建立的颤振模型在不同来流速度下的时域响应如图6所示，仿真获得颤振速度为1270 m/s，两者间的偏差仅为1.11%，证明基于刚柔耦合多体动力学建立的颤振模型精度符合要求。
![]()
图 6 典型舵面在不同来流速度下的第二阶广义位移时域响应Figure 6. Time domain response of the second-order generalized displacement of a typical rudder surface under different free stream velocities3.2 计及舵轴间隙的全动舵面颤振结果分析
舵轴转动自由度释放自由，舵轴拨片与拨叉采用碰撞模型，拨叉保持固定，即舵机模型不参与耦合，碰撞距离δ设定为舵轴转动角度（0.38°），采用传感器检测舵轴角度θ和舵面4节点位移
$\tilde x$ ，则舵面广义位移由各节点的修正位移x与模态逆矩阵相乘获得，各节点的修正位移为$$ x = \tilde x + l\sin \theta \text{，} $$ (23) 其中l为节点到舵轴沿舵面的位移。
基于计及舵轴间隙的全动舵面颤振模型，典型舵面在不同来流速度下的时域响应如图7和图8所示，可以发现：存在间隙时，在来流速度远低于颤振速度的工况下，全动舵面会产生结构极限环振荡现象（图7(a)、图8(a)）；在来流速度高于颤振速度的工况下，最初在短时间内出现既不收敛也不发散的结构极限环振荡现象，随后振动幅值逐步增大直至振荡发散（图7(b)、图8(b)）。
![]()
图 7 δ=0.38°时典型舵面在不同来流速度下的第二阶广义位移时域响应Figure 7. Time domain response of the second-order generalized displacement of a typical rudder surface under different free stream velocities for δ=0.38°![]()
图 8 δ=0.38°时典型舵面在不同来流速度下的舵轴转角时域响应Figure 8. Time domain response of rotating angle of a typical rudder surface under different free stream velocities for δ=0.38°采用上述相同参数，碰撞距离δ设定改为0.95°时，基于计及舵轴间隙颤振模型，典型舵面在不同来流速度下的时域响应如图9和图10所示，可以发现：此时来流速度高于颤振速度后出现既不收敛也不发散的结构振荡，且并未随仿真时间的延长而出现振荡发散现象，推测产生该现象的原因是仿真时间较短。
![]()
图 9 δ=0.95°时典型舵面在不同来流速度下的第二阶广义位移时域响应Figure 9. Time domain response of the second-order generalized displacement of a typical rudder surface under different free stream velocities for δ=0.95°![]()
图 10 δ=0.95°时典型舵面在不同来流速度下的舵轴转角时域响应Figure 10. Time domain response of rotating angle of a typical rudder surface under different free stream velocities for δ=0.95°综上可知，存在舵轴间隙时，全动舵面会在来流速度较低时产生极限环振荡现象，但对于颤振现象的发生具有减缓作用，且随着间隙的增大，全动舵面产生颤振发散现象的时刻会推后。
3.3 全动舵系统颤振结果分析
随后开展舵机系统对于全动舵面颤振特性的研究，此时舵轴转动自由度释放，舵轴拨片与拨叉采用多体动力学连接模型，拨叉运动由舵机给定，舵机模型参与振动耦合。同样采用传感器检测舵轴角度θ和舵面4节点位移
$\tilde x$ ，用于求解广义位移（参见式(23)）。基于全动舵系统颤振模型的不同来流速度下第二阶广义位移时域响应如图11和图12所示，分析发现：由于舵轴拨片未固定，采用伺服舵机驱动导致某些时刻出现舵轴转角偏离平衡位置较大距离的现象，并引发弹性振动幅值急剧增加；但随着舵轴转角逐步至0°，弹性振动的收敛和发散现象由来流速度决定。
![]()
图 11 全动舵系统在不同来流速度下的第二阶广义位移时域响应Figure 11. Time domain response of the second-order generalized displacement of an all-movable rudder under different free stream velocities![]()
图 12 全动舵系统在不同来流速度下的舵轴转角时域响应Figure 12. Time domain response of rotating angle of an all-movable rudder under different free stream velocities仿真表明：在舵轴偏离平衡位置时，舵机伺服系统给定控制指令迫使舵轴偏转角度逐步靠近0°，该控制效果由舵机伺服系统的控制参数给定，因此伺服控制系统中参数调整不合理时可能加剧颤振现象的发生，甚至降低颤振速度。
4. 结束语
本文基于柔性多体动力学方法建立全动舵系统颤振模型，采用Craig-Bampton方法建立典型舵面刚柔耦合降阶模型，采用激振器式集中气动力模拟基于偶极子网格法获得广义气动力。数值仿真结果分析表明，该系统模型能够快速与间隙、伺服系统等多种模型相互耦合，适用于研究各种因素对全动舵面动力学响应的影响。本文主要研究舵轴间隙、舵机对全动舵面颤振特性的影响，得到主要结论如下：
1）采用Craig-Bampton方法建立全动舵面动力学模型与激振器式集中气动力模型耦合可获得舵面颤振临界速度，且与商用软件预测对比的偏差小于2%；
2）全动舵面的舵轴存在间隙时，发生极限环振荡的临界速度远远低于模型不存在间隙时的情况，且来流速度大于颤振速度后，短时间内出现既不收敛也不发散的结构极限环振荡现象，且极限环振荡持续时间随着间隙增大而延长；
3）全动舵系统的舵机影响将导致某些时刻出现舵轴转角偏离平衡位置较大距离，由此引发弹性振动幅值急剧增加，意味着伺服控制系统中参数调整不合理时将可能加剧结构颤振，甚至降低颤振速度。
-
![]()
图 5 无阻尼典型舵面静载荷下时域响应的计算结果对比
Figure 5. Comparison of calculation results for the time domain response of a typical rudder surface under static load with no damping
![]()
图 6 典型舵面在不同来流速度下的第二阶广义位移时域响应
Figure 6. Time domain response of the second-order generalized displacement of a typical rudder surface under different free stream velocities
![]()
图 7 δ=0.38°时典型舵面在不同来流速度下的第二阶广义位移时域响应
Figure 7. Time domain response of the second-order generalized displacement of a typical rudder surface under different free stream velocities for δ=0.38°
![]()
图 8 δ=0.38°时典型舵面在不同来流速度下的舵轴转角时域响应
Figure 8. Time domain response of rotating angle of a typical rudder surface under different free stream velocities for δ=0.38°
![]()
图 9 δ=0.95°时典型舵面在不同来流速度下的第二阶广义位移时域响应
Figure 9. Time domain response of the second-order generalized displacement of a typical rudder surface under different free stream velocities for δ=0.95°
![]()
图 10 δ=0.95°时典型舵面在不同来流速度下的舵轴转角时域响应
Figure 10. Time domain response of rotating angle of a typical rudder surface under different free stream velocities for δ=0.95°
![]()
图 11 全动舵系统在不同来流速度下的第二阶广义位移时域响应
Figure 11. Time domain response of the second-order generalized displacement of an all-movable rudder under different free stream velocities
-
[1] 李家旭, 田玮, 谷迎松. 间隙非线性气动弹性系统颤振及控制问题研究进展[J]. 强度与环境, 2021, 48(4): 18-25 LI J X, TIAN W, GU Y S. Advance in the study on nonlinear flutter and control of aeroelastic system with freeplay nonlinearity[J]. Structure & Environment Engineering, 2021, 48(4): 18-25
[2] 杨享文, 武洁, 叶坤, 等. 高超声速全动舵面的热气动弹性研究[J]. 力学学报, 2014, 46(4): 626-630 DOI: 10.6052/0459-1879-13-415 YANG X W, WU J, YE K, et al. Study on aerothermoelasticity of a hypersonic all-movable control surface[J]. Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(4): 626-630 DOI: 10.6052/0459-1879-13-415
[3] 陈香言, 李明, 刘元清, 等. 高速风洞全动舵面颤振模型设计方法与试验研究[J]. 强度与环境, 2021, 48(5): 22-29 DOI: 10.19447/j.cnki.11-1773/v.2021.05.004 CHEN X Y, LI M, LIU Y Q, et al. Research on design method and experiment about equivalent flutter model for all-moving control surface of supersonic aircraft[J]. Structure & Environment Engineering, 2021, 48(5): 22-29 DOI: 10.19447/j.cnki.11-1773/v.2021.05.004
[4] TANG D, DOWELL E H. Aeroelastic response induced by free play, part 1: theory[J]. AIAA Journal, 2011, 49(11): 2532-2542 DOI: 10.2514/1.J051055
[5] TANG D, DOWELL E H. Aeroelastic response induced by free play, part 2: theoretical experimental correlation analysis[J]. AIAA Journal, 2011, 49(11): 2543-2554 DOI: 10.2514/1.J051056
[6] 李治涛, 韩景龙, 员海玮. 带有双线性非线性全动舵面气动弹性的数值与试验研究[J]. 振动与冲击, 2020, 39(19): 234-242 LI Z T, HAN J L, YUAN H W. Simulation and tests for aeroelasticity of a fully moving rudder surface with bilinear nonlinearity[J]. Journal of Vibration and Shock, 2020, 39(19): 234-242
[7] 苏华昌, 孙颖, 王晓晖, 等. 间隙非线性对舵面小幅值频率特性影响研究[J]. 强度与环境, 2020, 47(6): 31-378 SU H C, SUN Y, WANG X H, et al. Research on the influence of backlash nonlinearity on the frequency characteristics of small amplitude of rudder surface[J]. Structure & Environment Engineering, 2020, 47(6): 31-378
[8] 黄程德, 郑冠男, 杨国伟, 等. 基于CFD/CSD耦合含间隙三维全动舵面气动弹性研究[J]. 应用力学学报, 2018, 35(1): 1-7 HUANG C D, ZHENG G N, YANG G W, et al. Aeroelastic study of a three dimensional all-movable wing with free play using CFD/CSD coupling[J]. Chinese Journal of Applied Mechanics, 2018, 35(1): 1-7
[9] 杨超, 吴志刚. 导弹气动伺服弹性稳定性分析[J]. 飞行力学, 2000, 18(4): 1-5 DOI: 10.3969/j.issn.1002-0853.2000.04.001 YANG C, WU Z G. Aeroservoelastic stability of missile[J]. Flight Dynamics, 2000, 18(4): 1-5 DOI: 10.3969/j.issn.1002-0853.2000.04.001
[10] ZHANG X T, WU Z G, YANG C. New flutter-suppression method for a missile fin with an actuator[J]. Journal of Aircraft, 2013, 50(3): 989-994 DOI: 10.2514/1.C031961
[11] SHIN W H, LEE I, SHIN Y S, et al. Nonlinear aeroelastic analysis for a control fin with an actuator[J]. Journal of Aircraft, 2007, 44(2): 597-605 DOI: 10.2514/1.24721
[12] 丁伟涛, 吴志刚, 黄玉平, 等. 高超声速飞行器空气舵系统耦合特性分析及颤振抑制研究[J]. 导弹与航天运载技术, 2019(1): 84-88 DING W T, WU Z G, HUANG Y P, et al. Interaction analysis and flutter suppression research of hypersonic vehicles actuator-fin system[J]. Missiles and Space Vehicles, 2019(1): 84-88
[13] CRAIG R R, BAMPTON M C. Coupling of substructures for dynamic analyses[J]. AIAA Journal, 1968, 6(7): 1313-1319 DOI: 10.2514/3.4741
[14] OHTAKE K, ODEN J T, KIKUCHI N. Analysis of certain unilateral problems in von Karman plate theory by a penalty method-part 1. a variational principle with penalty[J]. Computer Methods in Applied Mechanics and Engineering, 1980, 24(2): 187-213 DOI: 10.1016/0045-7825(80)90045-6
[15] CHEN P C, LIU D D. Unsteady supersonic computations of arbitrary wing-body configurations including external stores[J]. Journal of Aircraft, 1990, 27(2): 108-116 DOI: 10.2514/3.45905
-
期刊类型引用(1)
1. 陈文权，马志赛，符杰，李亚南，张忠，丁千. 折叠变体飞行器的降阶建模与颤振分析. 强度与环境. 2024(05): 26-34 . 
百度学术
其他类型引用(1)
下载:
计量
- 文章访问数: 171
- HTML全文浏览量: 21
- PDF下载量: 33
- 被引次数: 2
