CN110837677A - 二元翼型非线性颤振时域模型的建模方法 - Google Patents

Abstract

本发明公开了一种二元翼型非线性颤振时域模型的建模方法，以一个舵叶连接一个扭簧的系统为对象，对系统进行简化，基于MSTMM推导弯扭耦合梁传递矩阵；确定各元件传递矩阵，拼装成总传递矩阵，建立系统整体动力学模型；令传递矩阵中刚心与质心重合，求解动力学模型，得到二元颤振模型的关键参数，即舵叶纯弯、纯扭频率；取3/4叶片展长处的截面，建立两元翼型非线性颤振模型；求解二元翼型非线性颤振模型，求得其时域响应；本发明采用MSTMM计算出二元颤振模型所需的关键动力学参数，即舵系统纯弯、纯扭频率，提出考虑间隙非线性和几何非线性的二元翼型颤振模型的建模思路，为舵系统非线性水弹性研究提供了一种有效的技术途径。

Description

二元翼型非线性颤振时域模型的建模方法

技术领域

本发明涉及一种建模方法，特别涉及二元翼型非线性颤振时域模型的建模方法，属于非线性振动技术领域。

背景技术

工程实际中，如水翼艇的水翼，潜艇中的舵翼等船舶构件作为动力学系统，结构参数设计不当时，在流场中高速运动会产生颤振现象导致结构破坏，或发生持续的弱振动现象诱发水噪声，降低水下航行器的隐蔽性。国内外学者对于舵系统或者水翼的水弹性研究，多数采用二元颤振模型，或者单个水翼加一根扭簧等简化模型来进行理论计算分析。对于二元颤振模型的计算中，关键在于如何获得二元颤振模型的计算参数。在早期的计算中，舵叶往往处理成刚体，采用舵轴的等效弹簧刚度、等效扭簧刚度作为二元颤振模型中的弹簧和扭簧，从而建立二元颤振模型。近期研究中，学者们逐渐将舵叶等水翼结构处理为柔性体，沿展向的所有剖面的翼型都是相同的，并假定每个翼型片条为刚体。水翼的弯曲和扭转变形分别采用两自由度水翼的沉浮和俯仰运动来模拟。对于非均直舵叶，一般取翼展方向70％-75％处的典型截面来建立二元颤振模型，通过二维模型来对弹性升力面建模。但是在这些学术论文中，几乎未提及如何获取水翼、控制面的纯弯、纯扭频率，以及如何将系统模型简化为二元颤振模型。因此，研究获取水翼纯弯、纯扭频率的方法具有一定意义。

舵系统的操纵系统部分往往由于长期工作磨损而产生间隙，这种间隙可能会引起结构持续发生微弱的、不衰减的振荡。这种现象并不会导致结构发生重大破坏，但会引起水噪声并降低水下航行器的隐蔽性。同时，舵叶严重的变形会使得结构呈现出明显的几何差异，位移和应变表现出非线性关系，即随着变形程度增大，刚度渐硬。传统线性颤振计算方法不能准确解决非线性系统的水动弹性问题。

多体系统传递矩阵法(MSTMM)由芮筱亭院士及其团队建立，用于多体动力学分析；该方法能实现一般多体系统动力学研究，且该方法可以实现复杂系统的快速建模与快速计算，其高计算效率方便运用于工程中。

发明内容

本发明的目的是提供一种二元翼型非线性颤振时域模型的建模方法，为舵系统非线性水弹性研究提供了一种有效的技术途径。

本发明的目的是这样实现的：一种二元翼型非线性颤振时域模型的建模方法，包括以下步骤：

步骤a、对系统进行简化，基于多体系统传递矩阵法推导弯扭耦合梁传递矩阵；

步骤b、确定各元件传递矩阵，拼装成总传递矩阵，建立系统多体动力学模型；

步骤c、令传递矩阵中刚心与质心重合，求解多体动力学模型，得到二元颤振模型的关键参数，即舵叶纯弯、纯扭频率；

步骤d、建立两元翼型非线性颤振时域分析模型；

步骤e、求解二元翼型非线性颤振模型，求得其时域响应，完成建模。

作为本发明的进一步限定，所述步骤a还包括：

根据舵叶的变形特性，将舵叶简化为多段不同参数的弯扭耦合梁，以舵叶弹性轴位置为x轴，以与扭簧相接的翼型弦线为z轴，以x轴、z轴交点为原点，y轴为经过原点垂直于xz平面的一条线，建立坐标系。

作为本发明的进一步限定，所述步骤a推导弯扭耦合梁传递矩阵具体包括：某段等截面舵叶的弯扭耦合梁模型中舵叶的长度为L，x_α是质心到弹性轴的距离，质心在Z轴正方向，则x_α为正；以弯曲和扭转振动为主，忽略其剪切效应，推导舵叶弯扭耦合梁传递矩阵，建立弯扭耦合梁的振动微分方程：

式中，m为单位长度质量，I_α为单位长度转动惯量，EI为弯曲刚度，GJ为扭转刚度，x为舵叶轴向位移，y为舵叶弯曲方向位移，θ_x为扭转角度，t为时间，b为舵叶弦长的一半，x_α质心到弹性轴的距离，且当质心在Z正方向，x_α为正；

由模态坐标转换，令

y(x,t)＝Y(x)sinωt,θ_x(x,t)＝Θ_x(x)sinωt (2)

式中，ω为圆频率；

将(2)带入(1)中，(1)式变为：

消去式(3)中的Y(x)或者Θ_x(x)，得到

式中，

W＝Y或者Θ (5)

引入无量纲长度

ξ＝x/L (6)

式(6)中L为划分后单元弯扭耦合梁长度；

式(4)可改写成无量纲形式

(D⁶+aD⁴-bD²-abc)W＝0 (7)

其中，

六阶微分方程(7)的通解可以表示为：

W(ξ)＝C₁coshαξ+C₂sinhαξ+C₃cosβξ+C₄sinβξ+C₅cosγξ+C₆sinγξ (9)

式中C₁-C₆是常数，并且

q＝b+a²/3

式(9)中的W(ξ)代表弯曲位移Y和扭转角度Θ_x在不同常数下的解；因此，

Y(ξ)＝A₁coshαξ+A₂sinhαξ+A₃cosβξ+A₄sinβξ+A₅cosγξ+A₆sinγξ (11)

Θ_x(ξ)＝B₁coshαξ+B₂sinhαξ+B₃cosβξ+B₄sinβξ+B₅cosγξ+B₆sinγξ (12)

式中A₁-A₆和B₁-B₆是两组不同的常数；

将式(11)和(12)代入式(3)可以确定常数有如下规律：

式中，

从式(11)、(12)可以得到弯曲角度Θ_z(ξ)，弯矩M_z(ξ)，剪切力Q_y(ξ)和扭矩M_x(ξ)的表达式：

因此可得

即Z(ξ)＝B(ξ)·a,a＝[A₁,A₂,A₃,A₄,A₅,A₆]^T，其中Z(ξ)为状态矢量，即方程(19)等式左边项，B(ξ)代表方程(19)等式右边的第一个矩阵；令ξ＝0，且根据方程(19)和多体系统传递矩阵法可将方程(19)写成Z_I＝B(0)·a形式；令ξ＝1得到Z_O＝B(1)·a＝B(1)B(0)^-1Z_I＝U_iZ_I；其中Z_O代表模态坐标下元件输出端的状态矢量，Z_I代表模态坐标下元件输入端的状态矢量；

得到弯扭耦合梁的传递矩阵为

U_i＝B(1)·B^-1(0) (20)

其中，

作为本发明的进一步限定，所述步骤b具体包括：基于MSTMM，状态矢量写为[X,Y,Θ_z,M_z,Q_x,Q_y,Θ_x,M_x]^T，X为沿坐标轴x位移对应的模态坐标列阵，Y为沿坐标轴y位移对应的模态坐标列阵，Θ_z为该点处相对于平衡位置相对于在z轴的角位移对应的模态坐标列阵，M_z为沿坐标轴z内力矩对应的模态坐标列阵，Q_x为沿坐标轴x内力对应的模态坐标列阵，Q_y为沿坐标轴y内力对应的模态坐标列阵，Θ_x为该点处相对于平衡位置相对于x轴的角位移对应的模态坐标列阵，M_x为沿坐标轴x内力矩对应的模态坐标列阵；

MSTMM中，总传递矩阵设置如下：

Z_n,n+1＝U_allZ_0,1 (21)

式中U_all＝U_n···U₃U₂U₁，U₁,U₂,U₃···U_n分别代表系统各元件的传递矩阵，Z_all＝[Z_n,n+1,···Z_1,2,Z_0,1]，Z_0,1,Z_1,2···Z_n,n+1分别代表系统各元件端点的状态矢量；

将扭簧定义为输入端，将离扭簧最远的弯扭耦合梁定义为输出端，扭簧到舵叶为传递方向；根据系统两端边界点，确定边界条件为；根据各段弯扭耦合梁参数确定其传递矩阵U_i(i＝2,3,...,8)，根据扭簧具体参数确定其传递矩阵U₁，建立系统总传递方程。

作为本发明的进一步限定，所述步骤c具体包括：将系统总传递方程中质心到弹性轴的距离都改为0(x_α≈0)，得到

和

即将弯扭耦合梁中的耦合项去掉，求解修改过后的动力学模型：

求解式(22)即可得到舵叶纯弯、纯扭频率。

作为本发明的进一步限定，所述步骤d具体包括：考虑舵叶的间隙非线性和几何非线性，根据模型建立两自由度翼型运动控制方程：

式中，m为单位展长舵叶的质量；b为舵叶弦长的一半；

为相对于弹性轴的单位展长转动惯量；r_α是水翼对刚心的回转半径；c_h和c_α分别是沉浮和俯仰阻尼；α表示俯仰运动，h表示沉浮运动；L_h为升力，取向上的方向为正方向；对弹性轴的水动力矩T_α则以抬头为正；考虑间隙非线性和几何非线性的非线性项F(h)、G(h)分别表示回复力和回复力矩，F(h)、G(α)为位移h和α的函数，其表达式为：

式中，k_h表示液压弹簧刚度系数，k_α表示扭转刚度系数；

非线性沉浮刚度系数，

非线性俯仰刚度系数；

ω_h和ω_α分别上一步骤中求出的纯弯、纯扭频率。

与现有技术相比，本发明所达到的有益效果：

(1)相比前人文献中均未给出二元翼型颤振模型中如何获得其纯弯、纯扭频率的方法，本发明提出基于多体系统传递矩阵法计算得到舵叶纯弯、纯扭频率的方法；

(2)本发明具有建模快速简单、计算速度快的特点，能快速计算出结果，解决了工程上仿真计算量过大、计算时间长的问题，为工程上类似的二元翼型颤振模型的快速建模与仿真提供参考；

(3)本发明在二元翼型颤振模型的建模中，同时考虑了间隙非线性和几何非线性，使颤振模型更接近实际，解决了传统线性颤振计算方法不能准确计算非线性系统的水动弹性问题。

附图说明

图1是本发明实例提供的一种二元翼型非线性颤振时域模型的建模流程示意图。

图2是舵叶模型简化示意图。

图3是弯扭耦合梁模型示意图。

图4是二元翼型非线性颤振模型示意图。

图5是二元翼型沉浮振幅随速度变化图。

图6是二元翼型俯仰振幅随速度变化图。

图7是当间隙为α_s＝0.0065,ξ_s＝0.0005，来流速度为V_non＝0.01的相轨图。

图8是当间隙为α_s＝0.0065,ξ_s＝0.0005，来流速度为V_non＝0.01的频谱分析图。

图9是来流速度下响应频率的分布图。

具体实施方式

下面结合附图对本发明作进一步描述；以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

本发明的一种二元翼型非线性颤振时域模型的建模方法，参见图1所示。现以某水下舵叶为例，根据一根扭簧连接一个舵叶的系统建立二元翼型非线性颤振模型，进行时域颤振分析。系统结构如图2所示。

系统具体参数如下：扭簧刚度为K_α＝1.23×10⁷N·m/rad，质心到弹性轴的距离x_α＝0.2889，水翼对刚心的回转半径r_α ²＝0.405，舵叶弦长的一半b＝0.9m，质量比μ＝0.403，刚心距离弦线中点的距离与半弦长的比值a＝-0.48，沉浮间隙无量纲量ξ_s＝0～0.0005，俯仰间隙无量纲量a_s＝0～0.0065rad。

具体实施步骤如下：

步骤一：对系统进行简化，基于多体系统传递矩阵法推导弯扭耦合梁传递矩阵。

根据舵叶的变形特性，将舵叶简化为7段不同参数的弯扭耦合梁，如图2所示，以舵叶弹性轴位置为x轴，以与扭簧相接的翼型弦线为z轴，以x轴、z轴交点为原点，y轴为经过原点垂直于xz平面的一条线，建立坐标系。如图3所示为某段等截面舵叶的弯扭耦合梁模型，该段舵叶的长度为L，x_α是质心到弹性轴的距离，质心在Z轴正方向，则x_α为正。

舵叶主要以弯曲和扭转振动为主，忽略其剪切效应，推导舵叶弯扭耦合梁传递矩阵，建立弯扭耦合梁的振动微分方程：

式中，m为单位长度质量，I_α为单位长度转动惯量，EI为弯曲刚度，GJ为扭转刚度，x为舵叶轴向位移，y为舵叶弯曲方向位移，θ_x为扭转角度，t为时间，b为舵叶弦长的一半，x_α质心到弹性轴的距离，且当质心在Z正方向，x_α为正；

由模态坐标转换，令

y(x,t)＝Y(x)sinωt,θ_x(x,t)＝Θ_x(x)sinωt (2)

式中，ω为圆频率；

将(2)带入(1)中，(1)式变为：

消去式(3)中的Y(x)或者Θ_x(x)，得到

式中，

W＝Y或者Θ (5)

引入无量纲长度

ξ＝x/L (6)

式(6)中L为划分后单元弯扭耦合梁长度；

式(4)可改写成无量纲形式

(D⁶+aD⁴-bD²-abc)W＝0 (7)

其中，

六阶微分方程(7)的通解可以表示为

W(ξ)＝C₁coshαξ+C₂sinhαξ+C₃cosβξ+C₄sinβξ+C₅cosγξ+C₆sinγξ (9)

式中C₁-C₆是常数，并且

q＝b+a²/3

式(9)中的W(ξ)代表弯曲位移Y和扭转角度Θ_x在不同常数下的解；因此，

Y(ξ)＝A₁coshαξ+A₂sinhαξ+A₃cosβξ+A₄sinβξ+A₅cosγξ+A₆sinγξ (11)

Θ_x(ξ)＝B₁coshαξ+B₂sinhαξ+B₃cosβξ+B₄sinβξ+B₅cosγξ+B₆sinγξ (12)

式中A₁-A₆和B₁-B₆是两组不同的常数；

将式(11)和(12)代入式(3)可以确定常数有如下规律：

式中，

从式(11)、(12)可以得到弯曲角度Θ_z(ξ)，弯矩M_z(ξ)，剪切力Q_y(ξ)和扭矩M_x(ξ)的表达式：

因此可得

即Z(ξ)＝B(ξ)·a,a＝[A₁,A₂,A₃,A₄,A₅,A₆]^T，其中Z(ξ)为状态矢量，即方程(19)等式左边项，B(ξ)代表方程(19)等式右边的第一个矩阵；令ξ＝0，且根据方程(19)和多体系统传递矩阵法可将方程(19)写成Z_I＝B(0)·a形式；令ξ＝1得到Z_O＝B(1)·a＝B(1)B(0)^-1Z_I＝U_iZ_I；其中Z_O代表模态坐标下元件输出端的状态矢量，Z_I代表模态坐标下元件输入端的状态矢量；

所以弯扭耦合梁的传递矩阵为

U_i＝B(1)·B^-1(0) (20)

其中，

步骤二：确定各元件传递矩阵，拼装成总传递矩阵，建立系统多体动力学模型；

基于MSTMM，状态矢量写为[X,Y,Θ_z,M_z,Q_x,Q_y,Θ_x,M_x]^T，X为沿坐标轴x位移对应的模态坐标列阵，Y为沿坐标轴y位移对应的模态坐标列阵，Θ_z为该点处相对于平衡位置相对于在z轴的角位移对应的模态坐标列阵，M_z为沿坐标轴z内力矩对应的模态坐标列阵，Q_x为沿坐标轴x内力对应的模态坐标列阵，Q_y为沿坐标轴y内力对应的模态坐标列阵，Θ_x为该点处相对于平衡位置相对于x轴的角位移对应的模态坐标列阵，M_x为沿坐标轴x内力矩对应的模态坐标列阵；

MSTMM中，总传递矩阵设置如下：

Z_n,n+1＝U_allZ_0,1 (21)

式中U_all＝U_n···U₃U₂U₁，U₁,U₂,U₃···U_n分别代表系统各元件的传递矩阵，Z_all＝[Z_n,n+1,···Z_1,2,Z_0,1]，Z_0,1,Z_1,2···Z_n,n+1分别代表系统各元件端点的状态矢量；

将扭簧定义为输入端，将离扭簧最远的弯扭耦合梁定义为输出端，扭簧到舵叶为传递方向。根据系统两端边界点，确定边界条件为：输入端Z_1,0＝[0,0,0,M_z1,Q_x1,Q_y1,0,M_x1]^T，输出端Z_9,0＝[X₉,Y₉,Θ_z9,0,0,0,Θ_x9,0]^T。Z_1,0为模态坐标下系统输入端的状态矢量，Z_9,0为模态坐标下系统输出端的状态矢量，

根据各段弯扭耦合梁参数确定其传递矩阵U_i(i＝2,3,...,8)，根据扭簧具体参数确定其传递矩阵U₁，建立系统总传递方程为：

步骤三：令传递矩阵中刚心与质心重合，求解多体动力学模型，得到二元颤振模型的关键参数，即舵叶纯弯、纯扭频率；

将系统总传递方程中质心到弹性轴的距离都改为0(x_α≈0)，得到

和

即将弯扭耦合梁中的耦合项去掉，求解修改过后的动力学模型，即可得到舵叶纯弯、纯扭频率。

具体求解方法如下：

将边界条件代入舵系统总传递方程可得到系统特征方程：

式中

为去掉零元素后得到的矩阵，为消去

中与

中零元素对应列得到的方阵；

由于

和U_all都只与系统的结构参数和固有频率有关，对于实际的线性振动系统，式(23)必有非零解，所以

求解式(24)即可求出舵叶的纯弯频率ω_h、纯扭频率ω_α。

步骤四：建立两元翼型非线性颤振时域分析模型；

取3/4叶片展长处的截面，考虑间隙非线性和几何非线性，建立二元水翼模型，如图4所示。在水翼弹性轴处固定一根刚度为K_h弹簧以及一根刚度为K_α的扭簧，弹性轴在翼弦中点前ab处，质心到弹性轴的距离为x_αb，水翼弦长为2b。水翼有两个自由度，即随弹性轴的沉浮运动h(向下为正)和绕弹性轴的俯仰转动α(抬头为正)；图4中，h_s是沉浮间隙，α_s是俯仰间隙。

根据模型建立两自由度翼型运动控制方程：

式中，m为单位展长舵叶的质量；b为舵叶弦长的一半；

为相对于弹性轴的单位展长转动惯量；r_α是水翼对刚心的回转半径；c_h和c_α分别是沉浮和俯仰阻尼；α表示俯仰运动，h表示沉浮运动；L_h(t)为升力，取向上的方向为正方向；对弹性轴的水动力矩T_α(t)则以抬头为正。考虑间隙非线性和几何非线性的非线性项F(h)、G(h)分别表示回复力和回复力矩，F(h)、G(α)为位移h和α的函数，其表达式为：

式中，k_h表示液压弹簧刚度系数，k_α表示扭转刚度系数；

非线性沉浮刚度系数，

非线性俯仰刚度系数。

根据Theodorsen非定常任意水动力流体理论，在不可压缩流中的二元水翼，水翼单位展长的升力L_h(t)(向上为正)和对刚心的俯仰力矩T_α(t)(迎风抬头为正)分别为：

式中，ρ为不可压缩流体的密度，σ为表征时间的量，

为无量纲时间，a为刚心距离弦线中点的距离与半弦长的比值，并在中点后为正；

为了方便推导非线性颤振方程，将升力和力矩整理为如下简洁形式：

式中，

步骤五：求解二元翼型非线性颤振模型，求得其时域响应；

由于式(31)中存在积分项，引入新的状态变化量如下：

根据式(25)、(26)、(27)、(30)、(31)和式(32)可以推导出状态空间中两自由度二元非线性水翼无量纲形式的水弹性方程为

式中，

其中，

式中，ω_h和ω_α分别上一步骤中求出的无耦合沉浮和俯仰固有频率，

为频率比，

为翼型绕弹性轴的无量纲回转半径，x_α为翼型质量中心到弹性轴的无量纲距离，ξ_ξ和ξ_α分别为沉浮和俯仰运动的阻尼比，μ为质量比，U为来流速度，V_non为无量纲来流速度，

为无量纲时间，R_ξ为非线性沉浮刚度系数无量纲量，R_α是非线性俯仰刚度系数无量纲量，ξ_s是沉浮间隙无量纲量，α_s是俯仰间隙无量纲量。

采用龙格库塔法求解式(33)及给定初始条件ξ(0)＝0，α(0)＝0.01rad(一般给α一个较小的初值)可以得到二元水翼的时域响应。

沉浮和俯仰振幅随着V_non变化图如图5和6所示；从图中可以看出，当俯仰间隙和沉浮间隙的值为α_s＝0.0065,ξ_s＝0.0005时，V_non＜0.07舵系统发生极限循环振动(LCO)现象，这种持续的振动可能诱发水噪声，降低水下航行器的隐蔽性。当V_non＞0.07时，系统只是发生了静变形，可能会降低舵面的操纵效率。在V_non≈0.08左右时，系统的静变形发生跳跃现象。因此，系统在包含间隙非线性的情况下，可能发生LCO和静变形等现象。

从图7可以看出，当间隙为α_s＝0.0065,ξ_s＝0.0005，来流速度为V_non＝0.01，舵系统发生了LCO现象，对应的频谱分析图如图8所示，在频率为10.09Hz左右有一个离散峰，离散峰处的能量相比其它处要高，表明此处容易激发出水噪声。

对图5、6工况下计算出的响应做频谱分析，获得响应频率，绘制出不同来流速度下响应频率的分布图，如图9所示。从图中可以看出，随着速度的增加，系统的响应频率逐渐降低。

同样，本方法对不同叶片建立二元翼型非线性颤振模型，只需修改简化模型和系统中具体参数即可。

本发明并不局限于上述实施例，在本发明公开的技术方案的基础上，本领域的技术人员根据所公开的技术内容，不需要创造性的劳动就可以对其中的一些技术特征作出一些替换和变形，这些替换和变形均在本发明的保护范围内。

Claims (6)

1.一种二元翼型非线性颤振时域模型的建模方法，其特征是，包括以下步骤：

步骤a、对系统进行简化，基于多体系统传递矩阵法推导弯扭耦合梁传递矩阵；

步骤b、确定各元件传递矩阵，拼装成总传递矩阵，建立系统多体动力学模型；

步骤c、令传递矩阵中刚心与质心重合，求解多体动力学模型，得到二元颤振模型的关键参数，即舵叶纯弯、纯扭频率；

步骤d、建立两元翼型非线性颤振时域分析模型；

步骤e、求解二元翼型非线性颤振模型，求得其时域响应，完成建模。

2.根据权利要求1所述的二元翼型非线性颤振时域模型的建模方法，其特征是，所述步骤a还包括：

根据舵叶的变形特性，将舵叶简化为多段不同参数的弯扭耦合梁，以舵叶弹性轴位置为x轴，以与扭簧相接的翼型弦线为z轴，以x轴、z轴交点为原点，y轴为经过原点垂直于xz平面的一条线，建立坐标系。

3.根据权利要求2所述的二元翼型非线性颤振时域模型的建模方法，其特征是，所述步骤a推导弯扭耦合梁传递矩阵具体包括：某段等截面舵叶的弯扭耦合梁模型中舵叶的长度为L，x_α是质心到弹性轴的距离，质心在Z轴正方向，则x_α为正；以弯曲和扭转振动为主，忽略其剪切效应，推导舵叶弯扭耦合梁传递矩阵，建立弯扭耦合梁的振动微分方程：

式中，m为单位长度质量，I_α为单位长度转动惯量，EI为弯曲刚度，GJ为扭转刚度，x为舵叶轴向位移，y为舵叶弯曲方向位移，θ_x为扭转角度，t为时间，b为舵叶弦长的一半，x_α质心到弹性轴的距离，且当质心在Z正方向，x_α为正；

由模态坐标转换，令

y(x,t)＝Y(x)sinωt,θ_x(x,t)＝Θ_x(x)sinωt (2)

式中，ω为圆频率；

将(2)带入(1)中，(1)式变为：

消去式(3)中的Y(x)或者Θ_x(x)，得到

式中，

W＝Y或者Θ (5)

引入无量纲长度

ξ＝x/L (6)

式(6)中L为划分后单元弯扭耦合梁长度；

式(4)可改写成无量纲形式

(D⁶+aD⁴-bD²-abc)W＝0 (7)

其中，

六阶微分方程(7)的通解可以表示为：

W(ξ)＝C₁coshαξ+C₂sinhαξ+C₃cosβξ+C₄sinβξ+C₅cosγξ+C₆sinγξ (9)

式中C₁-C₆是常数，并且

q＝b+a²/3

式(9)中的W(ξ)代表弯曲位移Y和扭转角度Θ_x在不同常数下的解；因此，

Y(ξ)＝A₁ coshαξ+A₂ sinhαξ+A₃ cosβξ+A₄ sinβξ+A₅ cosγξ+A₆ sinγξ (11)

Θ_x(ξ)＝B₁ coshαξ+B₂ sinhαξ+B₃ cosβξ+B₄ sinβξ+B₅ cosγξ+B₆ sinγξ (12)

式中A₁-A₆和B₁-B₆是两组不同的常数；

将式(11)和(12)代入式(3)可以确定常数有如下规律：

式中，

从式(11)、(12)可以得到弯曲角度Θ_z(ξ)，弯矩M_z(ξ)，剪切力Q_y(ξ)和扭矩M_x(ξ)的表达式：

因此可得

即Z(ξ)＝B(ξ)·a,a＝[A₁,A₂,A₃,A₄,A₅,A₆]^T，其中Z(ξ)为状态矢量，即方程(19)等式左边项，B(ξ)代表方程(19)等式右边的第一个矩阵；令ξ＝0，且根据方程(19)和多体系统传递矩阵法可将方程(19)写成Z_I＝B(0)·a形式；令ξ＝1得到Z_O＝B(1)·a＝B(1)B(0)^-1Z_I＝U_iZ_I；其中Z_O代表模态坐标下元件输出端的状态矢量，Z_I代表模态坐标下元件输入端的状态矢量；

得到弯扭耦合梁的传递矩阵为

U_i＝B(1)·B^-1(0) (20)

其中，

4.根据权利要求3所述的二元翼型非线性颤振时域模型的建模方法，其特征是，所述步骤b具体包括：基于MSTMM，状态矢量写为[X,Y,Θ_z,M_z,Q_x,Q_y,Θ_x,M_x]^T，X为沿坐标轴x位移对应的模态坐标列阵，Y为沿坐标轴y位移对应的模态坐标列阵，Θ_z为该点处相对于平衡位置相对于在z轴的角位移对应的模态坐标列阵，M_z为沿坐标轴z内力矩对应的模态坐标列阵，Q_x为沿坐标轴x内力对应的模态坐标列阵，Q_y为沿坐标轴y内力对应的模态坐标列阵，Θ_x为该点处相对于平衡位置相对于x轴的角位移对应的模态坐标列阵，M_x为沿坐标轴x内力矩对应的模态坐标列阵；

MSTMM中，总传递矩阵设置如下：

Z_n,n+1＝U_allZ_0,1 (21)

式中U_all＝U_n···U₃U₂U₁，U₁,U₂,U₃···U_n分别代表系统各元件的传递矩阵，Z_all＝[Z_n,n+1,···Z_1,2,Z_0,1]，Z_0,1,Z_1,2···Z_n,n+1分别代表系统各元件端点的状态矢量；

将扭簧定义为输入端，将离扭簧最远的弯扭耦合梁定义为输出端，扭簧到舵叶为传递方向；根据系统两端边界点，确定边界条件为；根据各段弯扭耦合梁参数确定其传递矩阵U_i(i＝2,3,...,8)，根据扭簧具体参数确定其传递矩阵U₁，建立系统总传递方程。

5.根据权利要求4所述的二元翼型非线性颤振时域模型的建模方法，其特征是，所述步骤c具体包括：将系统总传递方程中质心到弹性轴的距离都改为0(x_α≈0)，得到

和

即将弯扭耦合梁中的耦合项去掉，求解修改过后的动力学模型：

求解式(22)即可得到舵叶纯弯、纯扭频率。

6.根据权利要求5所述的二元翼型非线性颤振时域模型的建模方法，其特征是，所述步骤d具体包括：考虑舵叶的间隙非线性和几何非线性，根据模型建立两自由度翼型运动控制方程：

式中，m为单位展长舵叶的质量；b为舵叶弦长的一半；

为相对于弹性轴的单位展长转动惯量；r_α是水翼对刚心的回转半径；c_h和c_α分别是沉浮和俯仰阻尼；α表示俯仰运动，h表示沉浮运动；L_h为升力，取向上的方向为正方向；对弹性轴的水动力矩T_α则以抬头为正；考虑间隙非线性和几何非线性的非线性项F(h)、G(h)分别表示回复力和回复力矩，F(h)、G(α)为位移h和α的函数，其表达式为：

式中，k_h表示液压弹簧刚度系数，k_α表示扭转刚度系数；

非线性沉浮刚度系数，

非线性俯仰刚度系数；

ω_h和ω_α分别上一步骤中求出的纯弯、纯扭频率。

Images