1、http:/DOI:10.13700/j.bh.1001-5965.2021.0439基于动响应数据的大柔性机翼结构降阶方法谢长川，张铎耀，安朝*(北京航空航天大学航空科学与工程学院，北京100191)摘要：现代飞行器机翼柔性大，几何非线性问题不可忽略。基于动响应数据样本，基于谐波平衡和快速 Fourier 变换对结构动力学方程中的非线性刚度系数进行识别，建立非线性结构降阶模型。引入位移残量基模态，进行柔性机翼大变形的位移恢复。结合曲面涡格法和三维曲面插值方法搭建大柔性机翼几何非线性气动弹性分析框架。相比传统基于静力学数据回归分析的几何非线性结构降阶方法，该方法需要的载荷集数目小，提高了分析效
2、率。计算结果表明：与非线性有限元方法相比，非线性结构降阶模型准确度高，能够有效应用于大柔性机翼几何非线性静气动弹性分析，而传统的线性计算方法与非线性方法相比结果差异较大。关键词：几何非线性；快速傅里叶变换；参数识别；降阶模型；静气动弹性中图分类号：V211.47；V214.1文献标志码：A文章编号：1001-5965（2023）06-1319-12以 Helios“太阳神”无人机为代表的高空长航时无人机，无论是在军用方面还是在民用方面都有着非常良好的应用前景。通常该类型飞行器机翼需要利用大展弦比及高复合材料使用率来满足大升阻比气动特性要求，导致机翼刚度低变形大，几何非线性气动弹性问题严重1。结
3、构建模是气动弹性分析中的核心问题，准确的大变形结构建模方法是几何非线性气动弹性分析的基础。国内外研究人员针对气动弹性应用开发了不同的大变形结构建模方法。传统的位移基有限元方法模型适用度高，软件代码成熟，但自由度数高，求解效率较低。Hodges 精确梁模型2-3，应变基有限元方法4及有限段模型5等计算效率高，但不便于处理复杂模型，使用范围受限。结构降阶模型可以有效降低模型自由度规模，提高计算效率，同时保持对复杂模型的适用性6。结构动力学分析中常用的降阶方法包括静力凝聚（Guyan 缩聚）、动力凝聚及模态降阶等。线性动力学分析中的模态分析方法无法直接在几何非线性分析中应用，需要针对非线性问题特征进
4、行修改。Mcewan 等7在 2001 年给出了适用于大变形结构计算的非线性降阶方程形式及非线性刚度系数求解方法，Mignolet 和 Soize8随后补充推导了非线性项的来源。该方法选择线性结构模态作为降阶基函数，将结构动力学方程在模态空间中进行表达，忽略非定常项，利用非线性有限元分析中得到的静变形及载荷作为样本，采用回归分析方法来得到未知的非线性刚度系数。与非线性位移基有限元方法比较，可以大幅降低问题自由度数目，并可以在几乎所有的商用软件上不加修改的应用，保持了复杂模型适用性。针对大变形结构面内位移（展向位移）无法直接被低阶线性模态表征的缺陷，Hollamp 等9对该结构降阶方法进行了扩展
5、，利用将面内位移对刚度的影响缩减至弯曲模态中而并不求解直接求解面内位移的思路，发展了隐式缩减（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 fo
6、r 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此模型基础上利用模态能量准则进行了降阶基函数挑选，在结构动响应分析中取得了良好的效果。基于结
7、构降阶模型，Harmin 和 Cooper12利用偶极子格网法建立了几何非线性气动弹性分析模型，分析了大柔性机翼非线性颤振问题。An 等13-14考虑气动力随动效应，改进结构降阶模型回归分析方法，结合曲面涡格法分析了大柔性机翼几何非线性静/动气动弹性响应问题。Cesnik 等15结构降阶模型集成到基于 CFD 气动力程序的气动弹性框架中，进一步提高了分析精度。上述研究中，非线性刚度系数通过在方程静力形式上进行回归分析获得，需要大量的载荷及其静变形分析结果作为样本，过程繁琐，计算样本的时间成本较高。同时，忽略非定常项会使非线性刚度系数辨识过程丢失部分信息。在非线性系统参数辨识中，离散Fourie
8、r 变换，谐波平衡方法等频域法16-17及直接参数估计等时域法18应用广泛，基于结构动响应数据可以准确获得非线性系统信息。但相关方法在几何非线性结构刚度系数识别中应用较少。给出结构几何非线性动力学方程，基于动响应数据及谐波平衡法求解非线性刚度系数，建立大柔性机翼结构降阶模型。引入位移残差基函数恢复结构展向位移。与非线性有限元方法分析结果相比，计算准确度满足要求。结合曲面涡格法和曲面样条插值法分析大柔性机翼静气动弹性变形响应，验证结构降阶模型在几何非线性气动弹性分析中的适用性。1结构降阶模型几何非线性问题是典型的大位移-小应变-小应力问题，结构位移与应变之间是非线性关系而应力应变本构关系仍然是线
9、性的。在不考虑结构阻尼的情况下，利用动量守恒关系，结构域内的平衡方程可写为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）中的和对应于施
10、加在变形构型上并转换到参考构型上的体作用力和面力。式（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式中：为基函数分量；为待定参数即广义坐标。
11、利用 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)rX
12、pdXK(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令为结构自然模态集，基函1
13、320北 京 航 空 航 天 大 学 学 报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非线
14、性刚度系数识别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(
15、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)=
16、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参数辨识
