CN110889169A - 基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法 - Google Patents

基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法 Download PDF

Info

Publication number
CN110889169A
CN110889169A CN201911155418.9A CN201911155418A CN110889169A CN 110889169 A CN110889169 A CN 110889169A CN 201911155418 A CN201911155418 A CN 201911155418A CN 110889169 A CN110889169 A CN 110889169A
Authority
CN
China
Prior art keywords
equation
control surface
formula
nonlinear
bending
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201911155418.9A
Other languages
English (en)
Other versions
CN110889169B (zh
Inventor
陈东阳
顾超杰
孙振业
曹九发
李迺璐
朱卫军
杨华
庄舒青
刘锟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Yangzhou University
Original Assignee
Yangzhou University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Yangzhou University filed Critical Yangzhou University
Priority to CN201911155418.9A priority Critical patent/CN110889169B/zh
Publication of CN110889169A publication Critical patent/CN110889169A/zh
Application granted granted Critical
Publication of CN110889169B publication Critical patent/CN110889169B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Feedback Control In General (AREA)

Abstract

本发明公开了一种基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法，包括以下过程：基于多体系统传递矩阵法推导弯扭耦合梁传递矩阵；建立系统的总传递方程，并求解圆频率和振型；使用Theodorsen非定常流理论，建立舵面系统的运动控制方程；考虑间隙非线性和摩擦非线性，建立基于MSTMM的体动力学方程，得到系统非线性颤振模型；求解系统非线性颤振模型，得到舵面系统振动时域响应。本发明解决了线性颤振计算方法不能准确预测非线性系统颤振响应的问题，实现舵面系统非线性颤振响应的快速求解。

Description

基于多体系统传递矩阵法的舵面系统非线性颤振模型建模 方法
技术领域
本发明属于气动弹性仿真技术领域，具体涉及一种基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法。
背景技术
工程实际中，如飞机尾舵，潜艇舵翼等舵面系统作为动力学结构，在高速流场中可能会发生颤振现象。颤振是指弹性结构在均匀流体中受到流动力、弹性力和惯性力的耦合作用而发生的大幅度振动。颤振的产生会使结构产生疲劳损耗，甚至可能在短时间内导致结构破坏。因此，研究舵面系统的颤振响应具有一定意义。
同时，实际系统连接件中往往是存在间隙的，轴与轴承之间也可能发生间隙非线性运动以及摩擦引起的迟滞非线性。例如，舵面系统的操纵系统部分磨损而产生间隙，这种间隙可能会引起结构持续发生微弱的、不衰减的振荡，从而引起噪声并降低的隐蔽性。传统线性颤振计算方法不能准确解决非线性系统的问题，因此需要针对非线性系统建立非线性颤振模型。
发明内容
本发明的目的在于克服现有技术中的不足，提出了一种基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法，解决了线性颤振计算方法不能准确预测非线性系统颤振响应的问题，实现舵面系统非线性颤振响应的快速求解。
为解决上述技术问题，本发明提供了一种基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法，其特征是，包括以下过程：
基于多体系统传递矩阵法推导弯扭耦合梁传递矩阵；
建立系统的总传递方程，并求解圆频率和振型；
使用Theodorsen非定常流理论，建立舵面系统的运动控制方程；
考虑间隙非线性和摩擦非线性，建立基于MSTMM的体动力学方程，得到系统非线性颤振模型；
求解系统非线性颤振模型，得到舵面系统振动时域响应。
进一步的，所述基于多体系统传递矩阵法推导弯扭耦合梁传递矩阵包括：
将舵面简化为n段不同参数的弯扭耦合梁，建立弯扭耦合梁的振动微分方程：。
Figure BDA0002284671650000021
式中，m为单位长度质量，Iα为单位长度转动惯量，EI为弯曲刚度，GJ为扭转刚度，x表示x轴，y为弯曲方向位移，θx为扭转角度，t为时间，b为舵面横截面弦长的一半，xα质心到弹性轴的距离，且当质心在Z正方向，xα为正；
将物理坐标转化为模态坐标，令
y(x,t)＝Y(x)sinωt,θx(x,t)＝Θx(x)sinωt (2)
式中，ω为圆频率，Y(x)、Θx(x)为y(x,t)、θx(x,t)通过模态坐标转换得到。
将(2)带入(1)中，(1)式变为：
Figure BDA0002284671650000022
消去式(3)中的Y(x)或者Θx(x)，得到
Figure BDA0002284671650000031
式中，
W＝Y或者Θ (5)
引入无量纲长度：
ξ＝x/L (6)
式(6)中L为划分后单元弯扭耦合梁长度；
式(4)可改写成无量纲形式：
(D6+aD4-bD2-abc)W＝0 (7)
其中，
Figure BDA0002284671650000032
六阶微分方程(7)的通解可以表示为：
W(ξ)＝C1coshαξ+C2sinhαξ+C3cosβξ+C4sinβξ+C5cosγξ+C6sinγξ (9)
式中C1,......,C6是常数，并且
Figure BDA0002284671650000033
Figure BDA0002284671650000034
Figure BDA0002284671650000035
q＝b+a2/3
Figure BDA0002284671650000036
式(9)中的W(ξ)代表弯曲位移Y和扭转角度Θx在不同常数下的解；因此，
Y(ξ)＝A1coshαξ+A2sinhαξ+A3cosβξ+A4sinβξ+A5cosγξ+A6sinγξ (11)
Θx(ξ)＝B1coshαξ+B2sinhαξ+B3cosβξ+B4sinβξ+B5cosγξ+B6sinγξ (12)
式中A1,......,A6和B1,......,B6是两组不同的常数；
将式(11)和(12)代入式(3)可以确定常数有如下规律：
Figure BDA0002284671650000041
式中，
Figure BDA0002284671650000042
从式(11)、(12)可以得到弯曲角度Θz(ξ)，弯矩Mz(ξ)，剪切力Qy(ξ)和扭矩Mx(ξ)的表达式：
Figure BDA0002284671650000043
Figure BDA0002284671650000044
Figure BDA0002284671650000045
Figure BDA0002284671650000046
因此可得
Figure BDA0002284671650000047
即Z(ξ)＝B(ξ)·a,a＝[A1,A2,A3,A4,A5,A6]T，其中Z(ξ)为状态矢量，即方程(19)等式左边项，B(ξ)代表方程(19)等式右边的第一个矩阵；令ξ＝0，且根据方程(19)和多体系统传递矩阵法可将方程(19)写成ZI＝B(0)·a形式；令ξ＝1得到ZO＝B(1)·a＝B(1)B(0)-1ZI＝UiZI；其中ZO代表模态坐标下元件输出端的状态矢量，ZI代表模态坐标下元件输入端的状态矢量；
所以第i段弯扭耦合梁的传递矩阵为：
Ui＝B(1)·B-1(0) (20)
其中，
Figure BDA0002284671650000051
Figure BDA0002284671650000052
进一步的，所述建立系统的总传递方程包括：
根据各段弯扭耦合梁传递矩阵，建立系统总传递方程为：
Zo＝UallZi (21)
式中Uall＝Un…U3U2U1，U1,U2,U3…Un分别代表各段弯扭耦合梁的传递矩阵，Zi,Zo分别代表系统模态坐标下输入端、输出端的状态矢量；输入端
Figure BDA0002284671650000053
输出端Zo＝[Y,Θz,0,0,Θx,0]T。
进一步的，所述建立舵面系统的运动控制方程包括：
建立舵面系统的运动控制方程：
Figure BDA0002284671650000054
其中，
Figure BDA0002284671650000055
n表示将舵面系统沿x轴分为n段，Ln、Tαn表示舵面上第n段的升力和俯仰力矩；M、C、K分别为质量矩阵、阻尼矩阵、刚度矩阵；v为舵面坐标系中的物理坐标，
Figure BDA0002284671650000056
分别为v对时间t的一阶导数和二阶导数。
进一步的，所述系统非线性颤振模型包括：
根据多体系统传递矩阵法，由于气动力不解耦，令v＝V·q＝[V1 V2][q1 q2]T，可以将其转化到模态坐标系中，其中V1为系统第1阶振型，V2为第2阶振型，q为舵面系统的广义坐标；同时，由于系统连接件中可能存在间隙造成间隙非线性运动以及摩擦引起的迟滞非线性，在式(24)中加入非线性项Fnonlinear，同时将式(24)转换到模态坐标下，建立的考虑间隙、摩擦非线性的体动力学方程表示为：
Figure BDA0002284671650000061
式中，
Figure BDA0002284671650000062
取前两阶模态
Figure BDA0002284671650000063
Figure BDA0002284671650000064
k1、k2为比例阻尼系数；其中ω1、ω2分别为结构的第一阶圆频率和第二阶圆频率；
式中，
Figure BDA0002284671650000065
其中h是舵面系统的弯曲位移，α舵面系统的扭转角度，fh是弹性恢复力，fα是弹性恢复力矩，hs是沉浮方向间隙位移，αs2αs1代表两个角度，αs2到αs1是俯仰方向的间隙，其中fα中包含了干摩擦引起的非线性力矩，kh表示液压弹簧刚度系数，kα表示扭转刚度系数。
进一步的，所述求解系统非线性颤振模型，得到舵面系统振动时域响应包括：
由于式(22)是Theodorsen理论的任意运动时域气动力模型，是一个积分-微分方程，其中存在积分项，因此直接数值积分很繁琐。引入下列新的状态变量ωv以简化计算：
Figure BDA0002284671650000071
Figure BDA0002284671650000072
采用式(27)将式(22)的积分项展开，式(22)中将只包括ωv项，没有积分项存在。然后利用含有参变量积分的导数公式(28)对式(27)进行转换，可以得到状态向量ωv应满足的微分方程
Figure BDA0002284671650000073
Figure BDA0002284671650000074
将式(25)中
Figure BDA0002284671650000075
中含q,
Figure BDA0002284671650000076
的项移至方程的左侧合并同类项，剩下的F(t)留在方程的右侧，得到新的质量矩阵Mnew，新的刚度矩阵Dnew，新的阻尼矩阵Knew，最后方程可以变换为如式(29)所示形式：
Figure BDA0002284671650000077
其中：G为方程化简过程中产生的系数矩阵，状态空间的气动弹性方程为：
γi(t)＝[ωv1(t)ωv2(t)ωv3(t)ωv4(t)]T,i＝1,2…n
ωv(t)＝[γ1(t)γ2(t)…γn(t)]T,
Figure BDA0002284671650000081
Figure BDA0002284671650000082
Figure BDA0002284671650000083
Figure BDA0002284671650000084
根据式(29)建立状态空间的气动弹性方程为：
Figure BDA0002284671650000085
采用龙格库塔法求解式(30)，可以得到广义坐标q1,q2，结合由MSTMM计算出的振型V1,V2，进而求出物理坐标随时间的响应。
与现有技术相比，本发明所达到的有益效果：
(1)本发明提出一种基于多体系统传递矩阵法的舵面系统非线性颤振模型快速建模思路，相比于其他方法过长的计算时间，更能满足工程上快速计算的需求；
(2)本发明采用多体系统传递矩阵法可建立气动、水弹性体动力学模型和仿真方法，给工程上类似的舵面系统气动弹性、水弹性分析提供参考。
(3)本发明在舵面系统非线性颤振模型的建模中，同时考虑了间隙非线性和摩擦非线性，使颤振模型更接近实际，解决了传统线性颤振计算方法不能准确预测非线性系统颤振响应的问题。
附图说明
图1是本发明方法流程示意图；
图2是舵面系统模型示意图；
图3是弯扭耦合梁模型示意图；
图4是基于MSTMM计算出的舵面系统圆频率计算结果图；
图5是舵面系统在弯曲方向上的功率密度谱；
图6是舵面系统在扭转方向上的功率密度谱；
图7是舵面系统弯曲随时间变化的响应图；
图8是舵面系统扭转随时间变化的响应图；
图9是舵面系统弯曲二维相图；
图10是舵面系统扭转二维相图。
附图标记：1、舵面；2、扭簧。
具体实施方式
下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。
现有技术中，多体系统传递矩阵法(MSTMM)由芮筱亭院士及其团队建立，用于多体动力学分析。该方法能实现一般多体系统动力学研究，其高计算效率方便运用于工程中。运用于本发明中对舵面系统振动特性的计算，可以实现颤振模型的快速建模与快速计算。
本发明基于MSTMM计算出舵面系统固有振动特性，提出舵面系统非线性颤振模型建模思路。同等计算条件下，采用CFD/CSD(计算流体力学/计算结构动力学双向流固耦合)计算时间过长，难以满足工程上快速计算的需求，本方法更方便运用于工程中对舵面系统颤振时域响应的快速计算，为实际舵面结构的颤振计算提供参考。本发明的一种基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法，参见图1所示。具体实施步骤如下：
步骤一：建立舵面弯扭耦合梁模型，基于多体系统传递矩阵法推导弯扭耦合梁传递矩阵；
现以舵面左端加一根扭转刚度为1000N·m/rad的扭簧的系统为例，建立非线性颤振模型并进行时域分析。系统模型如图2所示，系统主要由舵面1和扭簧2组成。其中空间直角坐标系，以舵面弹性轴为x轴，以与扭簧连接的舵面一端方向为z轴，x轴与z轴交点为坐标原点，y轴过坐标原点垂直于xz平面。
根据舵面的变形特性，将舵面简化为n段不同参数的弯扭耦合梁。如图3所示为舵面弯扭耦合梁模型，建立坐标系以弯扭耦合梁弹性轴为x轴，以梁端面弦线为z轴，x轴与z轴交点为坐标原点，y轴过坐标原点垂直于xz平面。弯扭耦合梁模型中长度为L，xα是质心到弹性轴的距离，质心在Z轴正方向，则xα为正。以弯曲和扭转振动为主，忽略其剪切效应，建立弯扭耦合梁的振动微分方程：
Figure BDA0002284671650000101
式中，m为单位长度质量，Iα为单位长度转动惯量，EI为弯曲刚度，GJ为扭转刚度，x表示x轴，y为弯曲方向位移，θx为扭转角度，t为时间，b为舵面横截面弦长的一半，xα质心到弹性轴的距离，且当质心在Z正方向，xα为正。
将物理坐标转化为模态坐标，令
y(x,t)＝Y(x)sinωt,θx(x,t)＝Θx(x)sinωt (2)
式中，ω为圆频率，Y(x)、Θx(x)为y(x,t)、θx(x,t)通过模态坐标转换得到。
将(2)带入(1)中，(1)式变为：
Figure BDA0002284671650000111
消去式(3)中的Y(x)或者Θx(x)，得到
Figure BDA0002284671650000112
式中，
W＝Y或者Θ (5)
引入无量纲长度：
ξ＝x/L (6)
式(6)中L为划分后单元弯扭耦合梁长度。
式(4)可改写成无量纲形式：
(D6+aD4-bD2-abc)W＝0 (7)
其中，
Figure BDA0002284671650000113
六阶微分方程(7)的通解可以表示为：
W(ξ)＝C1coshαξ+C2sinhαξ+C3cosβξ+C4sinβξ+C5cosγξ+C6sinγξ (9)
式中C1,......,C6是常数，并且
Figure BDA0002284671650000114
Figure BDA0002284671650000121
Figure BDA0002284671650000122
q＝b+a2/3
Figure BDA0002284671650000123
式(9)中的W(ξ)代表弯曲位移Y和扭转角度Θx在不同常数下的解。因此，
Y(ξ)＝A1coshαξ+A2sinhαξ+A3cosβξ+A4sinβξ+A5cosγξ+A6sinγξ (11)
Θx(ξ)＝B1coshαξ+B2sinhαξ+B3cosβξ+B4sinβξ+B5cosγξ+B6sinγξ (12)
式中A1,......,A6和B1,......,B6是两组不同的常数。
将式(11)和(12)代入式(3)可以确定常数有如下规律：
Figure BDA0002284671650000124
式中，
Figure BDA0002284671650000125
从式(11)、(12)可以得到弯曲角度Θz(ξ)，弯矩Mz(ξ)，剪切力Qy(ξ)和扭矩Mx(ξ)的表达式：
Figure BDA0002284671650000126
Figure BDA0002284671650000127
Figure BDA0002284671650000128
Figure BDA0002284671650000129
因此可得
Figure BDA0002284671650000131
即Z(ξ)＝B(ξ)·a,a＝[A1,A2,A3,A4,A5,A6]T，其中Z(ξ)为状态矢量，即方程(19)等式左边项，B(ξ)代表方程(19)等式右边的第一个矩阵。令ξ＝0，且根据方程(19)和多体系统传递矩阵法可将方程(19)写成ZI＝B(0)·a形式。令ξ＝1得到ZO＝B(1)·a＝B(1)B(0)-1ZI＝UiZI。其中ZO代表模态坐标下元件输出端的状态矢量，ZI代表模态坐标下元件输入端的状态矢量。
所以第i段弯扭耦合梁的传递矩阵为：
Ui＝B(1)·B-1(0) (20)
其中，
Figure BDA0002284671650000132
Figure BDA0002284671650000133
U的下标i为上文中对舵面系统的分段，分了n段就可以取到n，i表示具体的某一段，各基本参数都表示第i段的参数，所以不标出i了，此处仅为推导传递矩阵形式，各段传递矩阵形式都相同，但参数根据不同形状的弯扭耦合梁都不同。
步骤二：建立系统的总传递方程，并求解固有频率和振型；
状态矢量写为[Y,Θz,Mz,Qy,Θx,Mx]T，其中Y为沿坐标轴y位移对应的模态坐标列阵，Θz为该点处相对于平衡位置在z轴的角位移对应的模态坐标列阵，Mz为沿坐标轴z内力矩对应的模态坐标列阵，Qy为沿坐标轴y内力对应的模态坐标列阵，Θx为该点处相对于平衡位置在x轴的角位移对应的模态坐标列阵，Mx为沿坐标轴x内力矩对应的模态坐标列阵。
根据各段弯扭耦合梁传递矩阵，建立系统总传递方程为：
Zo＝UallZi (21)
式中Uall＝Un…U3U2U1，U1,U2,U3…Un分别代表各段弯扭耦合梁的传递矩阵，Zi,Zo分别代表系统模态坐标下输入端、输出端的状态矢量；根据系统两端边界点，固定端位移为0，自由端力和力矩为0，因此可以确定边界条件为：输入端Zi＝[0,0,Mz,Qy,0,Mx]T，输出端Zo＝[Y,Θz,0,0,Θx,0]T。
求解总传递方程，得到舵面系统的振动特性。其中计算得到舵面系统的圆频率值如图4所示。图中横坐标为圆频率，纵坐标表示|Δ|值的大小，当|Δ|值接近于0时即可以求出圆频率，每条竖线下对应的就是该阶模态的圆频率。计算得到第一阶圆频率ω1为3.27，第二阶圆频率ω2为11.88。计算得到的圆频率和振型V运用于下一步颤振模型建模中。
步骤三：使用Theodorsen非定常流理论，建立舵面系统的运动控制方程；
根据Theodorsen非定常任意气动力流体理论，在不可压缩流中，舵面系统单位展长的升力L(t)(向上为正)和对刚心的俯仰力矩Tα(t)(迎风抬头为正)分别为：
Figure BDA0002284671650000151
Figure BDA0002284671650000152
Figure BDA0002284671650000153
式中，ρ为不可压缩流体的密度，σ为表征时间的量，
Figure BDA0002284671650000154
为无量纲时间，U为来流速度，b为半弦长，a为刚心距离弦线中点的距离与半弦长的比值，并在中点后为正。ω为圆频率，y表示沉浮位移，
Figure BDA0002284671650000155
为沉浮速度，
Figure BDA0002284671650000156
为沉浮加速度；θ表示俯仰位移，
Figure BDA0002284671650000157
表示俯仰速度，
Figure BDA0002284671650000158
表示俯仰加速度。
式中，
Figure BDA0002284671650000159
建立舵面系统的运动控制方程：
Figure BDA00022846716500001510
其中，
Figure BDA00022846716500001511
n表示将舵面系统沿x轴分为n段，Ln、Tαn表示舵面上第n段的升力和俯仰力矩；M、C、K分别为质量矩阵、阻尼矩阵、刚度矩阵；v为舵面坐标系中的物理坐标，
Figure BDA00022846716500001512
分别为v对时间t的一阶导数和二阶导数。
步骤四：考虑间隙非线性和摩擦非线性，建立基于MSTMM的体动力学方程，即系统非线性颤振模型；
根据多体系统传递矩阵法，由于气动力不解耦，令v＝V·q＝[V1 V2][q1 q2]T，可以将其转化到模态坐标系中，其中V1为系统第1阶振型，V2为第2阶振型，q为舵面系统的广义坐标。同时，由于系统连接件中可能存在间隙造成间隙非线性运动以及摩擦引起的迟滞非线性，在式(24)中加入非线性项Fnonlinear，同时将式(24)转换到模态坐标下，建立的考虑间隙、摩擦非线性的体动力学方程表示为：
Figure BDA0002284671650000161
式中，
Figure BDA0002284671650000162
取前两阶模态
Figure BDA0002284671650000163
Figure BDA0002284671650000164
k1、k2为比例阻尼系数。其中ω1、ω2分别为结构的第一阶圆频率和第二阶圆频率。
式中，
Figure BDA0002284671650000165
其中h是舵面系统的弯曲位移，α舵面系统的扭转角度，fh是弹性恢复力，fα是弹性恢复力矩，hs是沉浮方向间隙位移，αs2αs1代表两个角度，αs2到αs1是俯仰方向的间隙，其中fα中包含了干摩擦引起的非线性力矩，kh表示液压弹簧刚度系数，kα表示扭转刚度系数。
步骤五：求解平面非线性颤振模型，得到平面振动时域响应；
由于式(22)是Theodorsen理论的任意运动时域气动力模型，是一个积分-微分方程，其中存在积分项，因此直接数值积分很繁琐。引入下列新的状态变量ωv以简化计算：
Figure BDA0002284671650000171
Figure BDA0002284671650000172
采用式(27)将式(22)的积分项展开，式(22)中将只包括ωv项，没有积分项存在。然后利用含有参变量积分的导数公式(28)对式(27)进行转换，可以得到状态向量ωv应满足的微分方程。
Figure BDA0002284671650000173
Figure BDA0002284671650000174
将式(25)中
Figure BDA0002284671650000175
中含q,
Figure BDA0002284671650000176
的项移至方程的左侧合并同类项，剩下的F(t)留在方程的右侧，得到新的质量矩阵Mnew，新的刚度矩阵Dnew，新的阻尼矩阵Knew，最后方程可以变换为如式(29)所示形式：
Figure BDA0002284671650000177
其中：G为方程化简过程中产生的系数矩阵，状态空间的气动弹性方程为：
γi(t)＝[ωv1(t)ωv2(t)ωv3(t)ωv4(t)]T,i＝1,2...n
ωv(t)＝[γ1(t)γ2(t)…γn(t)]T,
Figure BDA0002284671650000181
Figure BDA0002284671650000182
Figure BDA0002284671650000183
Figure BDA0002284671650000184
根据式(29)建立状态空间的气动弹性方程为：
Figure BDA0002284671650000185
采用龙格库塔法求解式(30)，可以得到广义坐标q1,q2，结合由MSTMM计算出的振型V1,V2，进而求出物理坐标随时间的响应(位移，扭角随时间变化的响应)。
实施例
针对来流速度为25m/s时，计算扭转间隙为0.001rad的舵面系统动力学响应。得到舵面系统弯曲和扭转方向上功率谱密度如图5、6所示，从图中可看出舵面振动频率主要在0.78处。舵面系统弯曲和扭转方向上的动力学响应如图7、8所示，结合两张图可以得到舵面系统位移随时间的变化。图9、10是舵面系统颤振二维相图。从图中可以看出，舵面发生了极限循环振荡，说明在该条件下舵面已发生颤振。
同样，本方法对其他舵面系统建立非线性颤振模型，只需修改简化多体动力学模型和其中具体参数即可。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变型，这些改进和变型也应视为本发明的保护范围。

Claims (6)

1.一种基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法，其特征是，包括以下过程：
基于多体系统传递矩阵法推导弯扭耦合梁传递矩阵；
建立系统的总传递方程，并求解圆频率和振型；
使用Theodorsen非定常流理论，建立舵面系统的运动控制方程；
考虑间隙非线性和摩擦非线性，建立基于MSTMM的体动力学方程，得到系统非线性颤振模型；
求解系统非线性颤振模型，得到舵面系统振动时域响应。
2.根据权利要求1所述的一种基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法，其特征是，所述基于多体系统传递矩阵法推导弯扭耦合梁传递矩阵包括：
将舵面简化为n段不同参数的弯扭耦合梁，建立弯扭耦合梁的振动微分方程：
Figure FDA0002284671640000011
式中，m为单位长度质量，Iα为单位长度转动惯量，EI为弯曲刚度，GJ为扭转刚度，x表示x轴，y为弯曲方向位移，θx为扭转角度，t为时间，b为舵面横截面弦长的一半，xα质心到弹性轴的距离，且当质心在Z正方向，xα为正；
将物理坐标转化为模态坐标，令
y(x,t)＝Y(x)sinωt,θx(x,t)＝Θx(x)sinωt (2)
式中，ω为圆频率，Y(x)、Θx(x)为y(x,t)、θx(x,t)通过模态坐标转换得到。
将(2)带入(1)中，(1)式变为：
Figure FDA0002284671640000021
消去式(3)中的Y(x)或者Θx(x)，得到
Figure FDA0002284671640000022
式中，
W＝Y或者Θ (5)
引入无量纲长度：
ξ＝x/L (6)
式(6)中L为划分后单元弯扭耦合梁长度；
式(4)可改写成无量纲形式：
(D6+aD4-bD2-abc)W＝0 (7)
其中，
Figure FDA0002284671640000023
六阶微分方程(7)的通解可以表示为：
W(ξ)＝C1coshαξ+C2sinhαξ+C3cosβξ+C4sinβξ+C5cosγξ+C6sinγξ (9)
式中C1,......,C6是常数，并且
Figure FDA0002284671640000024
Figure FDA0002284671640000025
Figure FDA0002284671640000026
q＝b+a2/3
Figure FDA0002284671640000031
式(9)中的W(ξ)代表弯曲位移Y和扭转角度Θx在不同常数下的解；因此，
Y(ξ)＝A1coshαξ+A2sinhαξ+A3cosβξ+A4sinβξ+A5cosγξ+A6sinγξ (11)
Θx(ξ)＝B1coshαξ+B2sinhαξ+B3cosβξ+B4sinβξ+B5cosγξ+B6sinγξ (12)
式中A1,......,A6和B1,......,B6是两组不同的常数；
将式(11)和(12)代入式(3)可以确定常数有如下规律：
Figure FDA0002284671640000032
式中，
Figure FDA0002284671640000033
从式(11)、(12)可以得到弯曲角度Θz(ξ)，弯矩Mz(ξ)，剪切力Qy(ξ)和扭矩Mx(ξ)的表达式：
Figure FDA0002284671640000034
Figure FDA0002284671640000035
Figure FDA0002284671640000036
Figure FDA0002284671640000037
因此可得
Figure FDA0002284671640000038
即Z(ξ)＝B(ξ)·a,a＝[A1,A2,A3,A4,A5,A6]T，其中Z(ξ)为状态矢量，即方程(19)等式左边项，B(ξ)代表方程(19)等式右边的第一个矩阵；令ξ＝0，且根据方程(19)和多体系统传递矩阵法可将方程(19)写成ZI＝B(0)·a形式；令ξ＝1得到ZO＝B(1)·a＝B(1)B(0)-1ZI＝UiZI；其中ZO代表模态坐标下元件输出端的状态矢量，ZI代表模态坐标下元件输入端的状态矢量；
所以第i段弯扭耦合梁的传递矩阵为：
Ui＝B(1)·B-1(0) (20)
其中，
Figure FDA0002284671640000041
Figure FDA0002284671640000042
3.根据权利要求1所述的一种基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法，其特征是，所述建立系统的总传递方程包括：
根据各段弯扭耦合梁传递矩阵，建立系统总传递方程为：
Zo＝UallZi (21)
式中Uall＝Un…U3U2U1，U1,U2,U3…Un分别代表各段弯扭耦合梁的传递矩阵，Zi,Zo分别代表系统模态坐标下输入端、输出端的状态矢量；输入端
Figure FDA0002284671640000043
输出端Zo＝[Y,Θz,0,0,Θx,0]T。
4.根据权利要求1所述的一种基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法，其特征是，所述建立舵面系统的运动控制方程包括：
建立舵面系统的运动控制方程：
Figure FDA0002284671640000051
其中，
Figure FDA0002284671640000052
n表示将舵面系统沿x轴分为n段，Ln、Tαn表示舵面上第n段的升力和俯仰力矩；M、C、K分别为质量矩阵、阻尼矩阵、刚度矩阵；v为舵面坐标系中的物理坐标，
Figure FDA0002284671640000053
分别为v对时间t的一阶导数和二阶导数。
5.根据权利要求4所述的一种基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法，其特征是，所述系统非线性颤振模型包括：
根据多体系统传递矩阵法，由于气动力不解耦，令v＝V·q＝[V1 V2][q1 q2]T，可以将其转化到模态坐标系中，其中V1为系统第1阶振型，V2为第2阶振型，q为舵面系统的广义坐标；同时，由于系统连接件中可能存在间隙造成间隙非线性运动以及摩擦引起的迟滞非线性，在式(24)中加入非线性项Fnonlinear，同时将式(24)转换到模态坐标下，建立的考虑间隙、摩擦非线性的体动力学方程表示为：
Figure FDA0002284671640000054
式中，
Figure FDA0002284671640000055
取前两阶模态
Figure FDA0002284671640000056
Figure FDA0002284671640000057
k1、k2为比例阻尼系数；其中ω1、ω2分别为结构的第一阶圆频率和第二阶圆频率；
式中，
Figure FDA0002284671640000061
其中h是舵面系统的弯曲位移，α舵面系统的扭转角度，fh是弹性恢复力，fα是弹性恢复力矩，hs是沉浮方向间隙位移，αs2αs1代表两个角度，αs2到αs1是俯仰方向的间隙，其中fα中包含了干摩擦引起的非线性力矩，kh表示液压弹簧刚度系数，kα表示扭转刚度系数。
6.根据权利要求1所述的一种基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法，其特征是，所述求解系统非线性颤振模型，得到舵面系统振动时域响应包括：
Theodorsen理论的任意运动时域气动力模型，是一个积分-微分方程，其中存在积分项，引入下列新的状态变量ωv以简化计算：
Figure FDA0002284671640000062
Figure FDA0002284671640000063
采用式(27)将式(22)的积分项展开，式(22)中将只包括ωv项，没有积分项存在；然后利用含有参变量积分的导数公式(28)对式(27)进行转换，可以得到状态向量ωv应满足的微分方程
Figure FDA0002284671640000064
Figure FDA0002284671640000065
将式(25)中
Figure FDA0002284671640000071
中含q,
Figure FDA0002284671640000072
的项移至方程的左侧合并同类项，剩下的F(t)留在方程的右侧，得到新的质量矩阵Mnew，新的刚度矩阵Dnew，新的阻尼矩阵Knew，最后方程可以变换为如式(29)所示形式：
Figure FDA0002284671640000073
其中：G为方程化简过程中产生的系数矩阵，状态空间的气动弹性方程为：
γi(t)＝[ωv1(t) ωv2(t) ωv3(t) ωv4(t)]T,i＝1,2...n
ωv(t)＝[γ1(t) γ2(t)…γn(t)]T,
Figure FDA0002284671640000074
Figure FDA0002284671640000075
Figure FDA0002284671640000076
Figure FDA0002284671640000077
根据式(29)建立状态空间的气动弹性方程为：
Figure FDA0002284671640000078
采用龙格库塔法求解式(30)，可以得到广义坐标q1,q2，结合由MSTMM计算出的振型V1,V2，进而求出物理坐标随时间的响应。
CN201911155418.9A 2019-11-22 2019-11-22 基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法 Active CN110889169B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911155418.9A CN110889169B (zh) 2019-11-22 2019-11-22 基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911155418.9A CN110889169B (zh) 2019-11-22 2019-11-22 基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法

Publications (2)

Publication Number Publication Date
CN110889169A true CN110889169A (zh) 2020-03-17
CN110889169B CN110889169B (zh) 2020-10-16

Family

ID=69748436

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911155418.9A Active CN110889169B (zh) 2019-11-22 2019-11-22 基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法

Country Status (1)

Country Link
CN (1) CN110889169B (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111597715A (zh) * 2020-05-18 2020-08-28 中国核动力研究设计院 一种变截面管路自由振动特性的分析方法
CN112818580A (zh) * 2021-02-07 2021-05-18 上海机电工程研究所 基于扩充模态矩阵的间隙结构动力学模型降阶方法及系统
CN113050596A (zh) * 2021-03-12 2021-06-29 北京强度环境研究所 一种随机激励下空气舵模态参数准确获取方法
CN113408072A (zh) * 2021-06-25 2021-09-17 扬州大学 风力机柔塔系统固有振动特性快速建模及仿真方法
CN113591178A (zh) * 2021-07-07 2021-11-02 南京航空航天大学 变高度非对称梁纵-弯耦合振动的动力学建模方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040248182A1 (en) * 2003-06-09 2004-12-09 Protein Mechanics, Inc. Efficient methods for multibody simulations
CN104793497A (zh) * 2015-04-20 2015-07-22 天津理工大学 基于多体系统离散时间传递矩阵法的机器人动力学建模方法
CN104965991A (zh) * 2015-07-08 2015-10-07 中国人民解放军军械工程学院 基于传递函数的机翼颤振速度确定方法
CN107103103A (zh) * 2016-02-22 2017-08-29 上海机电工程研究所 二维间隙结构非线性气动弹性模型建模方法
CN108052772A (zh) * 2017-12-30 2018-05-18 北京航空航天大学 一种基于结构降阶模型的几何非线性静气动弹性分析方法
CN109656149A (zh) * 2018-12-10 2019-04-19 上海卫星装备研究所 星箭耦合多体系统动力学计算试验方法及系统

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040248182A1 (en) * 2003-06-09 2004-12-09 Protein Mechanics, Inc. Efficient methods for multibody simulations
CN104793497A (zh) * 2015-04-20 2015-07-22 天津理工大学 基于多体系统离散时间传递矩阵法的机器人动力学建模方法
CN104965991A (zh) * 2015-07-08 2015-10-07 中国人民解放军军械工程学院 基于传递函数的机翼颤振速度确定方法
CN107103103A (zh) * 2016-02-22 2017-08-29 上海机电工程研究所 二维间隙结构非线性气动弹性模型建模方法
CN108052772A (zh) * 2017-12-30 2018-05-18 北京航空航天大学 一种基于结构降阶模型的几何非线性静气动弹性分析方法
CN109656149A (zh) * 2018-12-10 2019-04-19 上海卫星装备研究所 星箭耦合多体系统动力学计算试验方法及系统

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
CHEN DONGYANG 等: "Dynamic modeling of sail mounted hydroplanes system- Part I: Modal characteristics from a transfer matrix method", 《OCEAN ENGINEERING》 *
CHEN DONGYANG 等: "Dynamic modeling of sail mounted hydroplanes system-part II: Hydroelastic behavior and the impact of structural parameters and free-play on flutter", 《OCEAN ENGINEERING》 *
张恩阳 等: "具有间隙非线性的全动舵系统的颤振分析", 《兵器装备工程学报》 *
肖清 等: "舵系统的颤振计算与分析", 《中国舰船研究》 *
蔡文佳: "亚音速全动控制面气动弹性分析", 《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111597715A (zh) * 2020-05-18 2020-08-28 中国核动力研究设计院 一种变截面管路自由振动特性的分析方法
CN112818580A (zh) * 2021-02-07 2021-05-18 上海机电工程研究所 基于扩充模态矩阵的间隙结构动力学模型降阶方法及系统
CN112818580B (zh) * 2021-02-07 2022-08-16 上海机电工程研究所 基于扩充模态矩阵的间隙结构动力学模型降阶方法及系统
CN113050596A (zh) * 2021-03-12 2021-06-29 北京强度环境研究所 一种随机激励下空气舵模态参数准确获取方法
CN113408072A (zh) * 2021-06-25 2021-09-17 扬州大学 风力机柔塔系统固有振动特性快速建模及仿真方法
CN113408072B (zh) * 2021-06-25 2023-10-13 扬州大学 风力机柔塔系统固有振动特性快速建模及仿真方法
CN113591178A (zh) * 2021-07-07 2021-11-02 南京航空航天大学 变高度非对称梁纵-弯耦合振动的动力学建模方法
CN113591178B (zh) * 2021-07-07 2024-03-29 南京航空航天大学 变高度非对称梁纵-弯耦合振动的动力学建模方法

Also Published As

Publication number Publication date
CN110889169B (zh) 2020-10-16

Similar Documents

Publication Publication Date Title
CN110889169B (zh) 基于多体系统传递矩阵法的舵面系统非线性颤振模型建模方法
CN110837677A (zh) 二元翼型非线性颤振时域模型的建模方法
CN110837678A (zh) 基于多体系统传递矩阵法的二元翼型频域颤振模型建模方法
Raveh CFD-based models of aerodynamic gust response
CN113806871B (zh) 一种考虑结构非线性的柔性飞行动力学建模方法
Xie et al. Static aeroelastic analysis including geometric nonlinearities based on reduced order model
CN113799136B (zh) 一种基于全状态反馈的机器人关节高精度控制系统及方法
CN105653781A (zh) 一种复合材料螺旋桨空泡性能的计算方法
CN110837676A (zh) 一种基于多体系统传递矩阵法的舵系统振动特性预测方法
CN110162826B (zh) 薄壁结构热气动弹性动响应分析方法
Lazarus et al. Multivariable active lifting surface control using strain actuation: analytical and experimental results
CN110287505B (zh) 飞行器稳定性分析方法
Kantor et al. Nonlinear structural, nonlinear aerodynamic model for static aeroelastic problems
Qian et al. Active flutter suppression of a multiple-actuated-wing wind tunnel model
CN109815637A (zh) 一种计算全柔机械臂动力学响应的仿真方法
CN113919081B (zh) 一种考虑惯性耦合的柔性飞行动力学建模与分析方法
Ritter et al. Enhanced modal approach for free-flight nonlinear aeroelastic simulation of very flexible aircraft
Liu et al. Continuous dynamic simulation for morphing wing aeroelasticity
CN110929336B (zh) 基于多体系统传递矩阵法求解三维机翼线性颤振速度的方法
Alsaif et al. Analytical and experimental aeroelastic wing flutter analysis and suppression
Cramer et al. Development of an aeroservoelastic model for gust load alleviation of the NASA common research model wind tunnel experiment
Merz A linear state-space model of an offshore wind turbine, implemented in the STAS wind power plant analysis program
Nguyen et al. Nonlinear Aeroelasticity of Flexible Wing Structure Coupled with Aircraft Flight Dynamics
Gao et al. Hopf bifurcation and feedback control of self-excited torsion vibration in the drive system
Marzocca et al. Unsteady aerodynamics in various flight speed regimes for flutter/dynamic response analyses

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant