基于动响应数据的大柔性机翼结构降阶方法.pdf

http:/DOI:10.13700/j.bh.1001-5965.2021.0439基于动响应数据的大柔性机翼结构降阶方法谢长川，张铎耀，安朝*(北京航空航天大学航空科学与工程学院，北京100191)摘要：现代飞行器机翼柔性大，几何非线性问题不可忽略。基于动响应数据样本，基于谐波平衡和快速 Fourier 变换对结构动力学方程中的非线性刚度系数进行识别，建立非线性结构降阶模型。引入位移残量基模态，进行柔性机翼大变形的位移恢复。结合曲面涡格法和三维曲面插值方法搭建大柔性机翼几何非线性气动弹性分析框架。相比传统基于静力学数据回归分析的几何非线性结构降阶方法，该方法需要的载荷集数目小，提高了分析效率。计算结果表明：与非线性有限元方法相比，非线性结构降阶模型准确度高，能够有效应用于大柔性机翼几何非线性静气动弹性分析，而传统的线性计算方法与非线性方法相比结果差异较大。关键词：几何非线性；快速傅里叶变换；参数识别；降阶模型；静气动弹性中图分类号：V211.47；V214.1文献标志码：A文章编号：1001-5965（2023）06-1319-12以 Helios“太阳神”无人机为代表的高空长航时无人机，无论是在军用方面还是在民用方面都有着非常良好的应用前景。通常该类型飞行器机翼需要利用大展弦比及高复合材料使用率来满足大升阻比气动特性要求，导致机翼刚度低变形大，几何非线性气动弹性问题严重1。结构建模是气动弹性分析中的核心问题，准确的大变形结构建模方法是几何非线性气动弹性分析的基础。国内外研究人员针对气动弹性应用开发了不同的大变形结构建模方法。传统的位移基有限元方法模型适用度高，软件代码成熟，但自由度数高，求解效率较低。Hodges 精确梁模型2-3，应变基有限元方法4及有限段模型5等计算效率高，但不便于处理复杂模型，使用范围受限。结构降阶模型可以有效降低模型自由度规模，提高计算效率，同时保持对复杂模型的适用性6。结构动力学分析中常用的降阶方法包括静力凝聚（Guyan 缩聚）、动力凝聚及模态降阶等。线性动力学分析中的模态分析方法无法直接在几何非线性分析中应用，需要针对非线性问题特征进行修改。Mcewan 等7在 2001 年给出了适用于大变形结构计算的非线性降阶方程形式及非线性刚度系数求解方法，Mignolet 和 Soize8随后补充推导了非线性项的来源。该方法选择线性结构模态作为降阶基函数，将结构动力学方程在模态空间中进行表达，忽略非定常项，利用非线性有限元分析中得到的静变形及载荷作为样本，采用回归分析方法来得到未知的非线性刚度系数。与非线性位移基有限元方法比较，可以大幅降低问题自由度数目，并可以在几乎所有的商用软件上不加修改的应用，保持了复杂模型适用性。针对大变形结构面内位移（展向位移）无法直接被低阶线性模态表征的缺陷，Hollamp 等9对该结构降阶方法进行了扩展，利用将面内位移对刚度的影响缩减至弯曲模态中而并不求解直接求解面内位移的思路，发展了隐式缩减（implicitcondensation,IC）方法，随后又进一步发展了 扩 展 隐 式 减 缩（implicitcondensationexpansion,ICE）方法来表征面内位移10。Guo 和 Przekop11在收稿日期：2021-08-03；录用日期：2021-08-27；网络出版时间：2021-09-1411：09网络出版地址： J.北京航空航天大学学报，2023，49（6）：1319-1330.XIE C C，ZHANG D Y，AN C.Reduced order method for large flexible wing structure based on dynamic response dataJ.Journal ofBeijing University of Aeronautics and Astronautics，2023，49（6）：1319-1330（in Chinese）.2023年6月北京航空航天大学学报June2023第49卷第6期JournalofBeijingUniversityofAeronauticsandAstronauticsVol.49No.6此模型基础上利用模态能量准则进行了降阶基函数挑选，在结构动响应分析中取得了良好的效果。基于结构降阶模型，Harmin 和 Cooper12利用偶极子格网法建立了几何非线性气动弹性分析模型，分析了大柔性机翼非线性颤振问题。An 等13-14考虑气动力随动效应，改进结构降阶模型回归分析方法，结合曲面涡格法分析了大柔性机翼几何非线性静/动气动弹性响应问题。Cesnik 等15结构降阶模型集成到基于 CFD 气动力程序的气动弹性框架中，进一步提高了分析精度。上述研究中，非线性刚度系数通过在方程静力形式上进行回归分析获得，需要大量的载荷及其静变形分析结果作为样本，过程繁琐，计算样本的时间成本较高。同时，忽略非定常项会使非线性刚度系数辨识过程丢失部分信息。在非线性系统参数辨识中，离散Fourier 变换，谐波平衡方法等频域法16-17及直接参数估计等时域法18应用广泛，基于结构动响应数据可以准确获得非线性系统信息。但相关方法在几何非线性结构刚度系数识别中应用较少。给出结构几何非线性动力学方程，基于动响应数据及谐波平衡法求解非线性刚度系数，建立大柔性机翼结构降阶模型。引入位移残差基函数恢复结构展向位移。与非线性有限元方法分析结果相比，计算准确度满足要求。结合曲面涡格法和曲面样条插值法分析大柔性机翼静气动弹性变形响应，验证结构降阶模型在几何非线性气动弹性分析中的适用性。1结构降阶模型几何非线性问题是典型的大位移-小应变-小应力问题，结构位移与应变之间是非线性关系而应力应变本构关系仍然是线性的。在不考虑结构阻尼的情况下，利用动量守恒关系，结构域内的平衡方程可写为19Xk(FijSjk)+0b0i=0 uiX 0（1）0b0iFijSjkuiXk00t0t0u0式中：为参考构型下的结构密度；为单元上的体作用力；为变形梯度；为第 2 类 P-K 应力张量；为结点变形分量；为结点参考位置坐标分量；表示结构域，式（1）使用了 Einstein 求和约定。在结构域边界上，面力在边界条件上，而指定边界位移，边界条件可表示为FijSjkn0k=t0iX t0（2）u=0X u0（3）n0t0u式中：为垂直于边界向外的单位法向量；为b0t0结构域的位移矢量。注意式（1）及式（2）中的和对应于施加在变形构型上并转换到参考构型上的体作用力和面力。式（1）式（3）给出了确定结构任意变形位置应力场和位移场的控制方程，式（3）表示参考构型下边界位移为 0。vi=vi(X)vi=vi(X)给定满足微分方程，即式（1）的函数以及边界条件，即式（2）和式（3）的函数，式（1）等效积分的弱形式可表示为w00vi uidX+w0viXk(FijSjk)dX=w00vib0idX+wt0vit0ids（4）uui(X,t)假定位移场 分量满足：ui(X,t)=Mn=1qn(t)U(n)i(X)（5）U(n)i(X)qn(t)vi=U(n)i(X)n=1,2,MM式中：为基函数分量；为待定参数即广义坐标。利用 Galerkin 方法求方程的近似解，给定，为所取基函数数量，即模型阶数。将式（5）代入式（4）中，经过代数运算可得到降阶后的结构动力学方程，其张量形式可表达为Mmn qn+K(1)mnqn+K(2)mnlqnql+K(3)mnlpqnqlqp=fm（6）式中：Mmn=w00U(m)iU(n)idX（7）K(1)mn=w0U(m)iXkCiklpU(n)lXpdX（8）K(2)mnl=12(K(2)mnl+K(2)lmn+K(2)nlm)（9）K(2)mnl=w0U(m)iXjCijkpU(n)rXkU(l)rXpdXK(2)lmn=w0U(l)iXjCijkpU(m)rXkU(n)rXpdXK(2)nlm=w0U(n)iXjCijkpU(l)rXkU(m)rXpdX（10）K(3)mslp=12w0U(m)iXjU(n)iXkCjkswU(l)rXsU(p)rXwdX（11）fm=w00U(m)ib0idX+wt0U(m)it0ids（12）Ciklp、Cijkp、Cjksw式中：为 4 阶材料本构张量，在线性本构关系下为常数。在上述降阶方程推导中，没有限制基函数的形式。结构线性模态即结构自然模态满足 Galerkin方法中对于基函数的要求，同时具有在结构分析中容易求解得到的特点，因此采用线性模态作为式（5）中的基函数是自然合理的选择。=1,2,M令为结构自然模态集，基函1320北 京 航 空 航 天 大 学 学 报2023年数可表示为U(m)=m（13）MmnK(1)mn代入式（7）、式（8）中并考虑结构自然模态正交性，可得广义质量项及广义刚度项为Mmn=Mmm=nMmn=0m,n（14）K(1)mn=Kmm=nK(1)mn=0m,n（15）MmKm式中：和分别为广义质量项和广义刚度项。物理空间位移与广义坐标关系可表示为qm=TmMuMm（16）M式中：为结构物理质量阵。式（6）可简化为Mm qm+Kmqm+K(2)mnlqnql+K(3)mnlpqnqlqp=fm（17）K(2)mnlK(3)mnlp式中：非线性刚度系数及为结构降阶模型中的待识别参数。2非线性刚度系数识别K(2)mnlK(3)mnlpqmfm采用基于谐波平衡法的多自由度非线性系统频域参数识别法识别式（17）中的非线性刚度系数项和。将式（17）中的广义坐标响应及广义力以等时间间隔进行离散，可以得到Mm qm(i)+Kmnqm(i)+K(2)mnlqn(i)ql(i)+K(3)mnlpqn(i)ql(i)qp(i)=fm(i)（18）i=1,2,NNqmQm(k)k=1,2,LL式中：，为时间离散总数。对广义坐标响应进行离散 Fourier 变换，其变换序列表示为，为频率离散总数，变换序列间满足：Qm(k)=N1i=0qm(i)ej2ikN（19）j式中：为虚数单位。qm(i)qn(i)ql(i)qn(i)ql(i)qp(i)fm(i)同理，将式（18）中的、非线性广义坐标二次 项、三 次 项及 模 态 力也进行离散 Fourier 变换可以得到以下频域序列：Qm(k)=N1i=0 qm(i)ej2ikN（20）Q(2)nl(k)=N1i=0qn(i)ql(i)ej2ikN（21）Q(3)nlp(k)=N1i=0qn(i)ql(i)qp(i)ej2ikN（22）Fm(k)=N1i=0fm(i)ej2ikN（23）经过离散 Fourier 变换后频域结构动力学方程可表示为MmQm(k)+KmQm(k)+K(2)mnlQ(2)nl(k)+K(3)mnlpQ(3)nlp(k)=Fm(k)k=0,1,L1（24）u qmQm(k)Qm(k)Q(2)nl(k)Q(3)nlp(k)Fm(k)在给定形式及大小的测试载荷下，通过非线性有限元软件计算一定时间内的结构时域响应数据样本，利用式（16）计算时域响应广义坐标并进行等时间间隔离散得到离散广义坐标响应样本，利用式（19）变换后可以得到频域序列样本，结合式（19）式（23），计算广义加速度、非线性广义坐标二次项、三次项及模态力频域序列样本分别为、及。将含有未知系数的部分保留在方程左侧，将已知量移至方程右侧，式（24）可写为回归问题形式：K(2)mnlQ(2)nl(k)+K(3)mnlpQ(3)nlp(k)=Fm(k)MmQm(k)KmQm(k)k=1,2,L（25）LLK(2)mnlK(3)mnlp式（25）一共包含 个线性方程，选择合适的值，求解式（25）的最小二乘解，即可识别出非线性刚度系数与。参数辨识流程如图 1 所示。给定外载荷非线性有限元软件计算时域动响应动响应数据样本结构模型等时间间隔离散频域样本结构动力学方程回归分析非线性刚度系数Qm,Fm,Qnl,Qnlp(2)(3)离散Fourier变换图1参数辨识流程Fig.1Flowchartofparameteridentification3结构大变形的位移恢复一般来说为了满足计算资源消耗的要求，结构第6期谢长川，等：基于动响应数据的大柔性机翼结构降阶方法1321v模型计算所选择的模态阶数有限。对于小变形问题，低阶结构模态往往可以满足结构分析需要，但对于大变形问题，结构模态的模态分量不能完备地描述结构在大变形情况下的所有位移分量。例如大柔性机翼结构在弯曲变形较大时，会产生明显的展向位移，该位移分量无法使用线性模态进行恢复。如图 2 所示，为较大的弯曲变形所导致的展向位移，低阶结构模态难以描述这一位移分量。类似的情况还包括因扭转变形较大所产生的弦向位移。v升力初始结构非线性变形线性模态图2弯曲大变形导致的展向位移Fig.2Spanwisedisplacementinducedbylargebendingdeformation引入位移残量基模态，将结构大变形表示为结构线性模态与位移残量基模态的叠加。将结构线性模态无法表征的位移残量表示为R=X q（26）RX=u(1),u(2),u(N)q=q(1),q(2),q(N)T式中：为残量矩阵；为物理空间下的等时间间隔离散动响应样本位移矩阵，每列由各离散时间点下的物理坐标响应组成；为等时间间隔离散的广义坐标响应样本，每列由各离散时间点下的广义坐标响应组成：q(i)=q1(i),q2(i),qM(i)T i=1,2,N（27）给定表征位移残量基模态的广义坐标向量形式为20Q(i)T=qs,1(i),qs,2(i),qs,G(i)=q21(i),q1(i)q2(i),q2M(i)（28）Q=Q(1),Q(2),Q(N)（29）qs,j(i)(j=1,2,G)G式中：为位移残量基模态广义坐标，由结构线性模态广义坐标的二阶多项式组合构成，为位移残量基模态广义坐标数目。位移残量基模态可表示为=RQ+（30）式中：上标+表示广义逆。该模态可以表征结构因为大变形导致展向位移及弦向位移等。大柔性结构的变形可表征为u=q+Q（31）位移恢复流程如图 3 所示。动响应样本位移矩阵X动响应样本广义坐标q残量矩阵R=Xq位移残量基模态位移残量基模态广义坐标Q位移恢复u=q+Q图3位移恢复流程Fig.3Flowchartofdisplacementrecovery4静气动弹性分析传统的气动弹性分析中一般采用平板气动力模型，不考虑机翼变形时升力面的曲面效应，当机翼变形较小时精度较高。对于具有几何非线性效应的大变形机翼情况，忽略曲面效应的平板气动力模型计算的结果则将对气动力分析造成较大影响1,13。采用曲面涡格法21-22作为气动力模型，结合非线性结构降阶模型进行大柔性机翼静气动弹性分析。4.1曲面涡格法xyz曲面涡格法以机翼中弧面作为气动面基准，在机翼的展向和弦向上划分面元，并在面元上取1/4 弦线的中点为气动力作用点，3/4 弦线中点处设定为边界条件控制点控制点，过气动力作用点布置涡环基本解，进行气动力求解。涡环单元由 4 段等强度直线涡首尾相接而成，自由涡由后缘涡格拖出，平行于来流方向。在控制点处满足不可穿透条件。机翼涡格划分及涡环布置情况如图 4 和图 5所示。曲面涡格法气动力求解坐标系为 轴沿气流来流方向，轴水平向右，轴由右手定则确定。zyOx图4机翼涡格划分Fig.4Vortexlatticedivisionofwing对于非后缘涡环网格而言，其涡环对空间任意一点的诱导速度可看作 4 段直线涡在该点诱导速度的线性叠加。位于后缘的涡环网格对空间任意一点的诱导速度可看作 3 段直线涡及两段半无限长涡线在该点诱导速度的线性叠加22。根据 Biot-Savart1322北 京 航 空 航 天 大 学 学 报2023年定理，所有涡格控制点处的诱导速度可表示为Vx=WxVy=WyVz=Wz（32）VxVyVzxyzWxWyWzxyz式中：、分别为所有涡格控制点处沿、轴方向的诱导速度分量的列向量；、分别为涡格沿、轴方向的诱导速度影响系数矩阵；为涡格强度向量。i根据 Neumann 边界条件，第 个涡格控制点处的不可穿透边界条件可以表达为(V+Vi)ni=0（33）VViinii式中：为无穷远处来流速度矢量；为所有涡环在第 个涡格控制点产生的诱导速度矢量；为第个涡格控制点处的法向矢量。i根据 Kutta-Jukowski 理论，第 个涡格气动网格力作用点的气动力为fAi=VFi（34）Fii式中：为大气密度；为第 个涡格力作用点（1/4 弦线中心点处）的总涡强。4.2结构/气动插值在结构气动耦合的过程中，位移信息及力信息需要准确地在气动和结构界面之间传递。二维插值方法不再适用于大变形气动弹性分析，需要采用三维曲面插值方法进行结构/气动的信息传递21,23。USXAUA结构结点变形矢量及气动网格节点变形矢量之间的关系可写为UA=GUS（35）G式中：为结构节点与气动节点之间的位移插值矩阵。FAFS根据虚位移的任意性，气动力与结构载荷的插值关系式为FS=GTFA（36）4.3静气动弹性分析方法大柔性机翼几何非线性静气动弹性分析中，气动载荷和结构运动与时间无关。静气动弹性分析中，首先给出计算工况及计算模型，采用非线性结构降阶模型及定常曲面涡格法，结合三维曲面插值方法及松耦合迭代流程计算机翼静气动弹性变形，最终收敛得到静气动弹性分析结果。位移收敛条件为?Ui+1StipUiStip?（37）UiStipi式中：为第 个迭代步中翼尖处结构结点位移；为收敛准则。大柔性机翼静气动弹性分析流程如图 6 所示。计算模型计算工况定常曲面涡格法计算曲面气动力力插值非线性结构降阶模型计算结果变形位移插值变形收敛?更新气动模型静气动弹性结果YN图6大柔性机翼静气动弹性分析流程Fig.6Staticaeroelasticanalysisflowchartoflargeflexiblewing5算例分析5.1大展弦比机翼模型xy以某单梁式平直机翼为算例模型，研究非线性结构降阶模型准确性。翼根固支取为坐标原点，以机翼前缘到机翼后缘方向为 轴正方向，翼根到翼尖方向为 轴正方向，z 轴满足右手定则。为了调整机翼的结构变形，在翼尖设置配重杆。机翼有限元模型如图 7 所示，机翼模型设计参数如表 1 所示。模型线性结构模态分析结果如表 2 及图 8 所示。该模型第 1 阶模态即一阶弯曲模态频率很低，接近 1Hz，模型柔性很大，同时由于翼尖配重杆的图5涡环布置情况Fig.5Vortexlatticearrangementzxy图7有限元模型Fig.7Finiteelementmodel第6期谢长川，等：基于动响应数据的大柔性机翼结构降阶方法1323存在，即使主梁截面扁平，一阶扭转模态频率及水平一弯模态频率均低于 30Hz，可以预见模型在气动力作用下将出现较大变形。5.2测试载荷的选取t选取适当的测试载荷作为非线性有限元软件输入数据，激发出结构非线性特性，计算时域动响应，提供辨识非线性刚度系数所需要的样本数据。采用“3211”多级方波信号形式的测试载荷作为输入。“3211”多级方波输入相比于其他输入，优势在于可以通过调整特征时间将频带移动到所希望激发的频带上，其数学模型为=0t 0,t1),t t1+7t,+)t t1,t1+3t)t t1+3t,t1+5t)t t1+5t,t1+6t)t t1+6t,t1+7t)（38）tt1t=0.1 s式中：为时间；为加载起始时间；为加载载荷幅值。选择机翼模型在 16m/s 来流风速，3迎角下的气动力作为加载载荷幅值。本文中，令，将“3211”形式测试载荷的有效频带调整至覆盖结构垂直一弯模态固有频率上，测试载荷对应垂直一弯模态的广义力的功率谱密度及方波输入如图 9所示，其频带较宽，低频弯曲模态频率附近能量大，满足动力学响应辨识需求。0.90.70.50.30.10.1226101418频率/Hz功率谱密度/(WHz1)(a)广义力功率谱密度642042610.10.30.701237 3 1 5 45678时间/s广义力/N(b)广义力时域变化图9垂直一弯模态的广义力功率谱密度及方波输入Fig.9Powerspectraldensityandsquarewaveinputofgeneralizedforceof1stbendingmode利用 MSC.Nastran 软件进行非线性时域动响应计算。时间步长 0.001s，共计算 7s。大柔性机翼模型翼尖处垂直位移响应如图 10 所示，机翼变形较大，超过机翼展长的 40%，测试载荷及动响应样本数据能够反映模型的非线性特性。5.3静力学验证首先对非线性结构降阶模型进行静力学变形计算正确性验证。选择机翼模型在 3迎角下，来流风速分别为 10m/s、16m/s 和 22m/s 时的气动力做表1机翼模型设计参数Table1Designparametersofwingmodel模型参数取值半翼展/mm1000弦长/mm100弹性轴位置50%机翼弦长主梁密度/(kgm3)7.75103主梁截面形状35mm1.5mm(矩形)配重杆长度/mm200配重杆质量/g62表2前六阶模态Table2Firstsixmodes模态编号模态振型模态频率/Hz1垂直一弯1.1792垂直二弯7.7243垂直三弯22.194一阶扭转22.955水平一弯27.476垂直四弯44.27zxyzxyzxyzxyzxyzxy(a)垂直一弯(b)垂直二弯(c)垂直三弯(d)一阶扭转(e)水平一弯(f)垂直四弯图8前六阶模态振型图Fig.8Vibrationmodediagramoffirstsixmodes1324北 京 航 空 航 天 大 学 学 报2023年为验证载荷。利用 MSC.Nastran 软件进行非线性静变形计算，与所建立的非线性结构降阶模型计算结果进行对比。机翼主梁翼尖垂直位移计算对比结果如表 3 所示。表3翼尖垂直位移对比Table3Comparisonoftipverticaldisplacement来流风速/(ms1)垂直位移/mm相对误差/%非线性有限元 非线性结构降阶模型1053.453.30.216153.7151.81.222313.0304.32.8从计算结果可以看出，非线性结构降阶模型准确度较高，相对误差在 1%左右。主梁垂直位移计算结果如图 11 所示。当计算载荷较小时，非线性方法计算结果与线性计算结果接近，小变形假设可以满足。而当载荷增大后，结构变形逐渐增大，非线性方法计算结果与线性计算结果之间的不同逐渐增大。采用小变形假设时结构变形计算不准确。非线性结构降阶模型计算结果与非线性有限元计算结果一致性较高。通过计算结果也可以发现，在垂直位移计算中，非线性降阶方法得到的非线性刚度相较于有限元方法略大，这与基于静力学响应数据建立降阶模型方法的结论是一致的7-8。主梁展向位移计算结果如图 12 所示。线性有限元方法并不能准确恢复机翼模型因为大弯曲变形导致的展向位移，而非线性有限元方法及非线性结构降阶模型计算结果一致性较高。当机翼模型垂直变形达到展长的 30%附近时，展向位移约为6%，已经能够对气动力计算产生明显影响20，使用非线性方法计算机翼变形对于气动弹性等分析具有重要意义。5.4动力学验证对非线性结构降阶模型进行动力学变形计算正确性验证。选择“3211”多级方波信号输入形式的验证载荷，采用机翼模型在 3迎角下，来流风速分别为10m/s 和12m/s 的气动力做为验证载荷幅值。利用有限元软件进行非线性动响应计算，与所建立的非线性结构降阶模型计算结果进行对比。机翼主梁翼尖垂直位移时域响应计算对比结果如图 13 所示。从计算结果可以发现，非线性结构降阶模型与非线性有限元方法计算结果一致性较好。机翼主梁翼尖展向位移时域响应计算对比结果如图14 所示，与静力学计算相同，线性方法不能计算展向时域响应，60040020002004006001012345678垂直位移/mm时间/s图10翼尖处垂直位移响应Fig.10Responseoftipverticaldisplacement非线性结构降价模型非线性有限元线性有限元非线性结构降价模型非线性有限元线性有限元非线性结构降价模型非线性有限元线性有限元60504030201001020002004006008001 000 1 200垂直位移/mm180140100602020垂直位移/mm40030035020025010015050050垂直位移/mm展向坐标/mm(a)来流风速为10 m/s20002004006008001 000 1 200展向坐标/mm(b)来流风速为16 m/s20002004006008001 000 1 200展向坐标/mm(c)来流风速为22 m/s图11不同来流风速下的主梁垂直位移Fig.11Verticaldisplacementofmainbeamwithdifferentflowvelocities第6期谢长川，等：基于动响应数据的大柔性机翼结构降阶方法1325非线性结构降阶模型计算结果准确度能够满足要求。5.5静气动弹性分析=0.5 mm?Ui+1StipUiStip?在涡格法计算中，沿展向和弦向划分 40 4 气动网格，计算 3迎角时机翼在不同来流风速下的静气动弹性变形，最终收敛准则条件为，即。作为对比，分别计算了基于非线性有限元方法和曲面涡格法的非线性静气动弹性分析方法及线性气动弹性计算方法在相同工况下的机翼静气动弹性变形。机翼主梁翼尖处各来流风速下的静气动弹性垂直位移计算结果如图 15 及表 4所示。线性方法非线性结构降价模型非线性有限元线性有限元非线性结构降价模型非线性有限元线性有限元非线性结构降价模型非线性有限元线性有限元00.40.81.21.62.020002004006008001 000 1 200展向位移/mm0481216展向位移/mm10100302050406070展向位移/mm展向坐标/mm(a)来流风速为10 m/s20002004006008001 000 1 200展向坐标/mm(b)来流风速为16 m/s20002004006008001 000 1 200展向坐标/mm(c)来流风速为22 m/s图12不同来流风速下的主梁展向位移Fig.12Spanwisedisplacementofmainbeamwithdifferentflowvelocities2001501005001005015020025015050501502501012345678垂直位移/mm时间/s1012345678垂直位移/mm时间/s(a)来流风速为10 m/s(b)来流风速为12 m/s非线性结构降价模型非线性有限元图13不同来流风速下的翼尖垂直位移Fig.13Verticaldisplacementofwingtipwithdifferentflowvelocities2024610812505101520251012345678展向位移/mm时间/s1012345678展向位移/mm时间/s(a)来流风速为10 m/s(b)来流风速为12 m/s非线性结构降价模型非线性有限元图14不同来流风速下的翼尖展向位移Fig.14Spanwisedisplacementofwingtipwithdifferentflowvelocities1326北 京 航 空 航 天 大 学 学 报2023年计算结果与非线性方法计算结果存在显著差异。基于非线性结构降阶模型的计算方法能够与基于非线性有限元方法计算结果保持很高的一致性，最大相对误差不超过 5%。机翼主梁翼尖处各风速下的静气动弹性展向位移计算结果如图 16 及表 5 所示。线性方法不能计算机翼展向位移，而基于非线性结构降阶模型的气动弹性分析方法能够给出准确计算结果，与基于非线性有限元方法的气动弹性分析方法相比，最大相对误差不超过 7%。对于气动弹性分析而言，机翼扭转角计算同样对气动载荷预测具有重要影响。机翼主梁翼尖处各来流风速下的静气动弹性扭转角计算结果如图 17 及表 6 所示。计算结果表明，与基于非线性表4静气动弹性垂直位移对比Table4Comparisonofstaticaeroelasticverticaldisplacement来流风速/(ms1)垂直位移/mm非线性有限元非线性结构降阶模型线性气动弹性分析1053.453.354.316153.7149.9164.022313.0296.6418.260504030201001020002004006008001 000 1 200垂直位移/mm18016014012010060408020020垂直位移/mm5003004001002000100垂直位移/mm展向坐标/mm(a)来流风速为10 m/s20002004006008001 000 1 200展向坐标/mm(b)来流风速为16 m/s20002004006008001 000 1 200展向坐标/mm(c)来流风速为22 m/s非线性结构降阶模型+曲面涡格法非线性有限元+曲面涡格法线性气动弹性分析非线性结构降阶模型+曲面涡格法非线性有限元+曲面涡格法线性气动弹性分析非线性结构降阶模型+曲面涡格法非线性有限元+曲面涡格法线性气动弹性分析图15不同来流风速下静气动弹性垂直位移Fig.15Staticaeroelasticverticaldisplacementwithdifferentflowvelocities00.40.81.21.62.020002004006008001 000 1 200展向位移/mm20246101281416展向位移/mm10100503020406070展向位移/mm展向坐标/mm(a)来流风速为10 m/s20002004006008001 000 1 200展向坐标/mm(b)来流风速为16 m/s20002004006008001 000 1 200展向坐标/mm(c)来流风速为22 m/s非线性结构降阶模型+曲面涡格法非线性有限元+曲面涡格法线性气动弹性分析非线性结构降阶模型+曲面涡格法非线性有限元+曲面涡格法线性气动弹性分析非线性结构降阶模型+曲面涡格法非线性有限元+曲面涡格法线性气动弹性分析图16不同来流风速下的静气动弹性展向位移Fig.16Staticaeroelasticspanwisedisplacementwithdifferentflowvelocities第6期谢长川，等：基于动响应数据的大柔性机翼结构降阶方法1327有限元方法的气动弹性分析方法相比，基于非线性降阶模型的气动弹性分析方法结果相对误差均小于线性方法计算结果。同时，从翼根到翼尖处，机翼扭转角计算误差逐渐增大，线性结果增加更加明显。机翼主梁翼尖静气动弹性位移随来流风速变化情况如图 18 所示。与基于非线性有限元方法的气动弹性分析方法相比，基于非线性结构降阶模型的气动弹性分析方法计算结果一致性较高，而线性气动弹性分析结果与非线性方法分析结果存在显著差异。结构降阶模型能够有效应用于机翼结构变形预测及非线性气动弹性分析。45015030008101214161820222481012141618202224垂直位移/mm10010203050604070展向位移/mm来流风速/(ms1)(a)垂直位移来流风速/(ms1)(b)展向位移非线性结构降阶模型+曲面涡格法非线性有限元+曲面涡格法线性气动弹性分析非线性结构降阶模型+曲面涡格法非线性有限元+曲面涡格法线性气动弹性分析图18不同来流风速下的翼尖静气动弹性位移Fig.18Staticaeroelastictipdisplacementwithdifferentflowvelocities6结论1)建立大柔性结构非线性动力学方程及降阶模型，基于动响应数据样本，采用谐波平衡法和快速 Fourier 变换的参数识别频域方法计算了降阶模表5静气动弹性展向位移对比Table5Comparisonofstaticaeroelasticspanwisedisplacement来流风速/(ms1)展向位移/mm非线性有限元 非线性结构降阶模型 线性气动弹性分析101.61.701613.613.702257.953.50表6静气动弹性扭转角对比Table6Comparisonofstaticaeroelastictwist来流风速/(ms1)扭转角/()非线性有限元 非线性结构降阶模型 线性气动弹性分析100.3430.3410.348160.9900.9391.045222.0781.9172.6450.0200.0150.0100.00500.0050.0200.0150.0100.00500.00520002004006008001 000 1 200扭转角/rad扭转角/rad0.050.010.030.040.0200.01扭转角/rad展向坐标/mm(a)来流风速为10 m/s20002004006008001 000 1 200展向坐标/mm(b)来流风速为16 m/s20002004006008001 000 1 200展向坐标/mm(c)来流风速为22 m/s非线性结构降阶模型+曲面涡格法非线性有限元+曲面涡格法线性气动弹性分析非线性结构降阶模型+曲面涡格法非线性有限元+曲面涡格法线性气动弹性分析非线性结构降阶模型+曲面涡格法非线性有限元+曲面涡格法线性气动弹性分析图17不同来流风速下的静气动弹性扭转角Fig.17Staticaeroelastictwistwithdifferentflowvelocities1328北 京 航 空 航 天 大 学 学 报2023年型中的非线性刚度。引入位移残量基模态，考虑大变形情况下位移恢复问题，能够计算因为弯曲及扭转大变形导致的机翼展向及垂直位移。2)对非线性结构降阶模型进行了静力学验证及动力学验证。与非线性有限元方法相比，垂直位移计算精度很高，在大变形情况下相对误差最大，为2.8%；展向位移计算精度较好，其最大误差为展长的0.17%。非线性结构降阶模型辨识得到的结构刚度略大。作为对比，线性有限元方法在大变形情况下并不能给出准确的垂直位移计算结果，在垂直位移上相对误差最大为 6.29%，并且无法计算结构展向位移。3)结合非线性结构降阶模型及曲面涡格法进行大柔性机翼几何非线性静气动弹性分析。与基于非线性有限元方法的气动弹性分析方法相比，计算精度较高，其中垂直位移方向上最大误差为5.4%。线性气动弹性方法计算结果与非线性方法相比差异较大，在垂直位移方向上最大误差为 33.6%。参考文献（References）谢长川,吴志刚,杨超.大展弦比柔性机翼的气动弹性分析J.北京航空航天大学学报,2003,29(12):1087-1090.XIECC,WUZG,YANGC.AeroelasticanalysisofflexiblelargeaspectratiowingJ.JournalofBeijingUniversityofAeronauticsandAstronautics,2003,29(12):1087-1090（inChinese）.1HODGESDH,DOWELLEH.Nonlinearequationsofmotionfortheelasticbendingandtorsionoftwistednonuniformrotorblades:NASA-TN-D-7818R.Washington,D.C.:NASA,1974.2HODGESDH.ReviewofcompositerotorblademodelingJ.AI-AAJournal,1990,28(3):561-565.3SIMOJC.Afinitestrainbeamformulation.Thethree-dimension-aldynamicproblem.PartIJ.ComputerMethodsinAppliedMech-anicsandEngineering,1985,49(1):55-70.4CAST