变体飞机通过机翼折叠可实现机翼面积200％的变化^[1]。这种大幅度变形必然会引起机翼的气动力以及结构动力学特性的变化。因此需要对变体飞机机翼在折叠变形过程中的动力学特性进行研究。Scarlett等^[2]和Reich等^[3]结合有限元方法、多体动力学、气动力模型,建立了折叠机翼柔性多体动力学模型。韩景龙课题组^[4]通过对偶极子网格法的有理近似,并结合MSC. NASTRAN与ADAMS,建立了多体动力学的瞬态响应。此时使用的商业软件不便于结构模型与气动力模型的耦合。Selitrennik等^[5]通过引入时变的坐标变换矩阵、虚质量子结构综合方法以及CFD计算结果,形成了状态空间下的气动弹性模型,但是该模型中CFD的计算量较大。Zhao等^[6]建立了机翼在折叠变形的动力学方程,预测了机翼在不同变形速度下的动力学响应,但是该模型没有考虑气动力的作用,因此无法预测气动力对于机翼在折叠变形过程中动力学响应的影响。

本文采用多体动力学中的浮动坐标系法,结构动力学中的模态综合法,引入有理函数近似导出机翼折叠变形时非定常气动力,建立了一组微分代数方程来模拟机翼的折叠运动,通过约束违约稳定法^[7]获得了机翼折叠运动的动力学响应。该控制方程适用类似的多段机翼结构,并且可以改变约束方程,获得不同变形方式下的动力学响应。

1 机翼折叠变形过程中的动力学控制方程

机翼的几何模型如图 1 所示。该模型由机身,内翼和外翼3个子结构组成,分别记为子结构A、B、C。机翼折叠运动的瞬态响应,可以看作是刚体运动和弹性变形的叠加。而浮动坐标系法就是将柔性体的运动分解为刚性运动和相对浮动坐标系的弹性变形^[8]。另外,对于机翼这样的复杂结构,采用有限元离散方法所得的自由度庞大,不利于动力学方程的求解。而子结构模态综合法是对复杂结构进行动力学建模的有效手段。因此,多体动力学中浮动坐标系法和结构动力学中的模态综合法的结合是解决该问题比较好的方法^[6]。如图 1所示,oxyz为惯性坐标系,用于描述系统的运动,o_Ax_Ay_Az_A、o_Bx_By_Bz_B、o_Cx_Cy_Cz_C分别为浮动坐标系,用于描述柔性体的弹性变形。θ^A、θ^B、θ^C分别为子结构A、B、C的折叠角。

结合浮动坐标系法^[8]以及模态综合法^[9],机翼的子结构的动力学方程为

式中,Mⁱ为子结构的广义质量矩阵,包含柔性体平动惯量、柔性体平动与转动的惯性耦合项、柔性体转动与变形的惯性耦合项、柔性体转动惯性张量、柔性体转动与变形的惯性耦合项以及柔性体弹性变形的惯性项^[8]。Dⁱ、Kⁱ分别为广义阻尼矩阵和广义刚度矩阵。考虑气动力作用时,广义阻尼矩阵与广义刚度矩阵包含了气动阻尼和气动刚度。F_cⁱ为约束力矢量。F_gⁱ为子结构的重力矢量。F_vⁱ为与速度二次方有关的广义惯性力矢量。F_aⁱ为气动力载荷。qⁱ为子结构的广义位移,其具体的形式为

式中,Rⁱ为子结构的浮动坐标系的原点位置。q_fⁱ为模态位移。在机翼折叠变形的过程中,子结构的广义位移包含了刚体运动和弹性变形,气动阻尼和气动刚度的存在,因此广义质量矩阵、广义阻尼矩阵和广义刚度矩阵均是一个时变的矩阵。整个系统是一个时变的系统。

1.1 机翼结构的模态综合中的模态集选取

对于机翼结构,需要采用有限元方法对结构进行离散,提取模态信息,创建广义质量矩阵中柔性体弹性变形的惯性项,再根据柔性多体动力学理论建立动力学控制方程。子结构模态综合将复杂结构划分为若干子结构,建立子结构的模态集及模态坐标,利用界面的协调条件将各个子结构组装得到整体结构动力学方程。此时所选取的模态集直接决定着柔性体变形量的准确性。因此,这里详细给出模态集的确定。采用有限元离散化形成的无阻尼动力学方程为

式中,分别为子结构的质量矩阵、刚度矩阵、界面处的力矢量。u为节点位移矢量。固定界面模态综合法是常用的综合方法^[9]。首先将子结构的自由度分为2个集合,即为结构内部节点位移u_I和边界节点位移u_J。将子结构的无阻尼动力学方程按该集合进行分块

式中,为边界处的内力矢量。引入模态坐标pⁱ以及模态矩阵Φⁱ

式中,Φ为固定界面的主模态矩阵,Φ_c为约束模态矩阵。此时模态矩阵包含子结构的6个刚体模态,不能直接应用于机翼折叠变形的瞬态响应分析。因此,需要进行矩阵变换。为了求得变换矩阵,对下列公式进行特征值求解

式中,,ω为结构的频率。计算出特征向量按质量矩阵进行正则化,得到特征向量φ_k。由这些特征向量得到变换矩阵T,正交化后的模态矩阵为

此时刚体模态就会显式地表达出来,剔除这些刚体模态,便可以建立机翼的结构动力学方程。在工程中,阻尼矩阵可以通过模态阻尼比给出。此时,引入阻尼矩阵,机翼子结构结构动力学方程为

式中,Φ^iT·uⁱ。

1.2 机翼折叠变形过程中的时变气动载荷

采用偶极子法可以获得频域下的非定常气动力,通过有理近似可以导出任意运动情况下的非定常气动力。因此给时域非定常气动力的模拟提供了思路。作为机翼折叠变形中的气动载荷初步研究,本文采用偶极子网格法获得了机翼的气动力。根据虚功原理,可以得到作用在结构网格节点的等效力

式中,为广义气动力影响系数矩阵,k为减缩频率,b为参考长度。ρ为密度,V为来流速度。D_aic为非定常气动力影响系数矩阵。G_s1和G_s2为位移插值转换矩阵,G_s为力插值转换矩阵,S为气动网格的面积矩阵。因此,F_a与子结构的广义位移qⁱ有关,也是一个时变的矢量。

1.3 机翼折叠变形的动力学控制方程

引入作用在结构网格节点的等效力

此时Q为频域气动力。为了进行机翼折叠变形时的瞬态响应分析,这里采用最小状态法进行频域气动力的有理函数近似获得时域内的气动力。此时有

式中,A₀、A₁、A₂、D_s和E_s,均为最小状态法近似中的未知矩阵,可以通过最小二乘法求得,R_s为气动力滞后矩阵。q_a为气动力状态变量。

对于机翼的折叠变形运动,存在着运动约束和驱动约束,可以表达成如下形式

式中

f_d^B 表示期望的折叠角。此时,子结构A与C均保持水平的状态。

考虑到约束方程,引入拉格朗日乘子λ,此时约束力F_c可以通过-CT_qλ来表达,C_q为约束方程的雅克比矩阵。最终机翼折叠变形的控制方程为 式中,M=blockdiag(M^A,M^B,M^C),D=blockdiag(D^A,D^B,D^C),K=blockdiag(K^A,K^B,K^C)。可以看出该动力学控制方程是微分代数方程,对于该方程可以采用约束违约稳定法来求解^[8]。将约束方程对时间进行2次求导,可以得到约束加速度方程,此时可以得到

式中,。α和β为正数,对于违约进行修正^[8]。

2 机翼折叠变形中的瞬态动力学响应分析

为了预测机翼折叠变形的动力学响应,令折叠角θ^B的变化符合以下规律

式中,120°代表机翼折叠变形结束时,子结构B所要达到的角度。T_s代表变形结束时所需要的时间。

本文将各个子结构简化为等厚度的铝板,在MSC.NASTRAN中被离散成四节点的CQUAD4单元,其单元属性为PSHELL,这样单元的每个节点有6个自由度,即3个平动自由度u_x、u_y、u_z,以及3个转动自由度θ_x、θ_y、θ_z。机翼的有限元模型如图 2所示。通过DMAP语言来获得子结构的质量矩阵和刚度矩阵,对于每个子结构均保留6阶模态,每阶所对应的模态阻尼比为0.01。因此最终得到机翼折叠变形的动力学控制方程,即公式(14)有46个方程,采用约束违约稳定法解进行求解,本文中取α=5,β=10,整个程序的编制均在MATLAB中完成。

2.1 柔性机翼在重力作用下瞬态响应

首先,不考虑气动力的作用。即在公式(14)中忽略与气动力有关的项。翼尖法向位移响应如图 3所示。

当T_s=2 s,T_s=6 s时,分别在1.27 s、3.81 s附近,θ^B=90°,此时子结构B与子结构A与C垂直,翼尖法向位移达到最大值。明显地,柔性对翼尖法向位移响应有较大影响。在整个变形中,翼尖的位移衰减收敛,并且变形结束时,位移稳定值要比刚性机翼时的位移稳定值要小。由于在公式(14)中引入了结构阻尼,位移响应呈收敛的趋势。T_s越大,意味着变形越缓慢,翼尖的法向位移变化也较缓慢。

对于刚性机翼,作用在子结构B上合力在机翼折叠变形开始和结束时的值均等于子结构B的重力,即4.89 N。而作用在子结构B上的合力矩,在变形运动开始时大约为341.36 Nmm,在变形结束时保持在-157.18 Nmm。初始值和稳定值均等于子结构B的重力产生的力矩。作用在子结构C上的合力,初值和稳定值基本为子结构C的重力,即11.564 N。作用在子结构C上的力矩,其初值和稳定值基本一样,均为2 885 Nmm。这是因为机翼在折叠变形的过程中,子结构C始终保持水平,所以子结构C的重心相对于其浮动坐标系没有很大的变化。在变形的过程中,无论是作用在子结构B上的合力、合力矩,还是作用在子结构C的力、力矩,均发生变化。

柔性对作用在子结构B上的合力、合力矩,以及作用在子结构C上的力、力矩有较大的影响。与刚性相比,作用在柔性子结构B上的合力有较大振荡,已经改变了在刚性状态下的变化趋势,但是随着时间的增加,合力呈收敛的趋势。对于合力矩而言,也有一定的震荡。作用在柔性子结构C上的合力、合力矩有较大的振荡,均呈衰减收敛的趋势。重力是一种体力,柔性体重力产生的广义力与弹性变形有关。所以在折叠变形时弹性变形发生变化,另外,弹性力与阻尼力也会对合力和合力矩产生一定的影响。

2.2 柔性机翼在重力和气动力作用下瞬态响应

其次,考虑在重力和气动力同时作用下柔性机翼的瞬态响应,在公式(14)中,ρ=1.226 kg／m³,M_∞=0.0,分别计算了不同来流速度下的瞬态响应如图 5和图 6所示。

不同来流速度下翼尖法向位移响应如图 5所示。当T_s=2 s时,在1.27s附近,θ^B=90°,翼尖法向位移达到最大值。与柔性机翼的位移响应有些类似,在不同的来流速度下翼尖的法向位移有一定的振荡。但从局部的放大图中可以看出,随着来流速度的增加,折叠变形结束时的振荡减小。这主要是由于气动阻尼的存在,也表明气动力有辅助折叠变形的作用。

从图 6中可以看出,在气动力的作用下,随着来流速度的增加,作用在子结构B上的合力、合力矩,以及作用在子结构C上的力、力矩均呈衰减振荡。考虑气动力时,气动力与弹性位移相关,且是分布力,因此相对于重力产生的力矩而言,是较小的,因此与柔性机翼的结果类似。由于气动阻尼的存在,力与力矩均有所减小。在一定的来流速度范围内,气动力有助于机翼的折叠变形。

3 结 论

结合多体动力学中的浮动坐标系法与结构动力学中模态综合法,本文建立了机翼折叠变形过程中的动力学控制方程,该控制方程考虑了一些非线性项,如与速度二次方有关的广义惯性力,可以用于分析机翼在缓慢变形和快速折叠变形的瞬态响应。分析了机翼柔性对瞬态响应的影响。将有理函数近似导出的折叠变形运动中的非定常气动力近似表达式引入到动力学控制方程中,分析了柔性机翼的瞬态气动弹性响应。可以得出以下结论:

1) 运用浮动坐标系法与模态综合法建立的折叠机翼动力学控制方程。该控制方程适用类似的多段机翼结构,并且可以改变约束方程,可以获得不同变形方式下的瞬态响应;

2) 对机翼折叠变形过程越缓慢,得到的瞬态响应越平稳;对于柔性折叠机翼而言,柔性对机翼在折叠变形过程的瞬态响应影响较大;

3) 考虑气动力时,在一定的流速范围内,来流速度越大,在变形结束后翼尖位移响应的振荡越小。因此气动力在一定的程度上有助于机翼的折叠变形运动。