基于几何精确应变梁的柔性机翼气动弹性建模与仿真

1月前浏览2889

导读：柔性机翼在气动载荷作用下产生较大变形，几何非线性因素不容忽视。利用机翼柔性特点，通常采用梁模型进行结构建模。从几何精确梁理论出发，结合Hamilton原理推导了几何非线性梁的动力学平衡方程。不同于经典的位移基有限元，采用梁广义应变作为插值变量，得到广义质量阵、广义阻尼阵、刚度阵及载荷列向量，建立非线性应变梁模型。结合 Newmark 数值算法和牛顿-拉夫森（Newdon-Raphson）迭代法建立了动力学方程求解算法。针对典型算例，开展静、动力学分析，并分别与有限元软件的仿真结果进行对比。

结果表明：在相当计算精度下，所建模型收敛特性更好。为进一步验证所建模型在工程实际应用中的精度和有效性，开展针对典型大展弦比机翼的地面静力试验。试验表明：仿真结果与激光位移计和光纤传感设备测得的变形值具有较高的一致性，验证了所建模型具有较高精度。

本文节选自北京航空航天大学-航空科学与工程学院许秋怡，孟杨和李书共同完成，公开发表在北京航空航天大学学报的学术论文《基于应变的几何非线性梁建模与分析》，欢迎专家、学者和工程师批评指正。诚邀读者关注10月13日20时，2024航空航天设计仿真专题报告第六期，将邀请本文作者孟杨老师做《低空经济之无人机结构强度仿真与测试》线上讲座。详情见后文。

扫码观看和回放

写在文前

本文通过引入空间位置矢量和局部坐标基向量构成的节点构形向量，结合 Cayley-Hamilton 定理构建显式的构形-应变关系，简化构形向量与应变向量间的雅可比关系推导，最终得到考虑剪切效应的几何非线性应变梁模型。给出考虑端部定向载荷及随动载荷作用下的大变形分析求解方法，开展典型对象的静、动力学分析，并分别与有限元软件的仿真计算结果进行对比。此外，开展针对大展弦比机翼的地面静力试验，进一步验证本文模型的有效性和计算精度。

一、理论建模

首先明确所涉及的梁模型具有以下 2 条基本假设：①应力-应变本构关系保持线性；②梁截面为刚性剖面，即忽略翘曲变形。图1为梁空间构形示意图。几何精确梁理论中，梁模型采用梁中心参考线和过参考线上每一点的梁截面表征。参考线并没有严格界定，可沿截面形心、质心、剪心等，选取方式会导致截面刚度矩阵发生变化。如此，梁的变形就可以被描述成梁参考线（一条空间曲线）的变形和横截面的空间平移及旋转运动。

1、构形向量定义及运动关系描述

定义参考坐标系 B，基向量为。该研究关注固支边界，参考坐标系 B 始终保持不变。给定沿梁中心线的弧长参数， L 为梁的长度，梁在任意时刻的空间构形可由 2 方面决定：

1）梁中心参考线的空间位置向量：。

2）梁截面局部坐标系G的基向量：，。基向量分别为s处横截面内正交的 2 个单位向量（通常沿截面主轴方向），则由右手法则定义：。局部坐标系G的基向量均表示在参考坐标系B中。注意，未变形时与梁中心参考线相切，即

考虑剪切变形，式 (1) 在变形后可能不再满足，即不一定与梁中心线相切。根据图1，梁内任意一点 P 的空间位置矢量可表示为：

式中：为向量从局部坐标系 G 到参考坐标系B中的坐标变换阵；向量为s处中心参考线到梁截面上点的位置向量，其表示在局部坐标系 G 中的分量列阵为。

为了方便后续推导，定义描述梁弹性变形的12×1 的列向量为

该向量为构形向量，为 (s，t) 的函数，为使表达式简洁，后文表述中将省略 (s，t)。任意点 P 的空间位置可表示为

式中：，零元素均代表 3 阶零阵，I 为3阶单位矩阵，该矩阵由梁初始构形决定，与结构变形无关。因此，只需确定构形向量，即可确定梁内任一点的空间位置。

式 (4) 对时间 t 求导可以得到任意点 P 的速度为

式 (5) 中各个向量均在参考坐标系中。

2、应变-构形关系

需要说明的是，本文所述“应变”并非工程应变，而是梁中心参考线的广义应变，类似曲线论中曲率概念。Reissner[5] 从梁内力、力矩平衡关系及几何关系出发，根据虚功原理，可得空间梁的应变-构形关系为

式中：为对弧长坐标s的导数；的 3 个分量分别为拉伸应变和沿方向的剪切应变，记为力应变；的 3 个分量分别为绕的扭转应变、绕和的弯曲应变，记为力矩应变，其物理含义为 t 时刻的曲率向量与初始曲率之差；；中的下标 0 表示初始时刻，为向量的反对称矩阵，满足：

由于局部坐标系 G 的基向量均在参考系 B 中，根据坐标变换阵的定义[18] 有如下关系成立：

记时刻 t 构形下梁的曲率向量为，3 个分量分别为绕的扭转曲率及绕这 2 个方向的弯曲曲率，且有。

将式 (8) 代入式 (6)，并引入式 (3) 定义的构形向量，可得应变-构形关系为

式中：为剪切应变和弯曲应变构成的12×12方阵，可表示为

考虑了剪切效应的应变阵，而在文献 [15]中，该矩阵中不包含剪切应变。

3、动力学方程

根据 Hamilton 原理可得

式中：和为任意 2 个固定时刻；T为单位长度的动能；U为单位长度的内能；为单位长度上的外力所作虚功，为变分符号算子。本节分别对和进行推导。

（1）动能变分

单位长度的动能T为

式中：ρ为密度；A(s)为弧长 s 处梁截面。将式 (5)

代入式（12）可得

求变分得

式中：为截面惯性阵，下标 CS 表示截面，有

其中：为单位长度的质量（线密度）；为在局部坐标系 G 中的截面质心位置；为截面惯量特性。

（2）内能变分

广义铁木辛柯梁（考虑剪切效应）应变能变分为

式中：分别为局部坐标系下作用在梁截面上的内力和内力矩。根据应力-应变本构关系，有

式中：均为3×3方阵，构成梁截面刚度矩阵。

（3）外力虚功

单位长度外力虚功为

式中：分别为作用在梁截面上的合力和合力矩；为虚位移；为虚转动。

对于飞行器而言，外力常包含随动力，即外力方向随结构状态发生变化，如气动力、推力等。随动力作用下系统是非保守的。通过坐标变换阵把相应的力矢从局部坐标系 G 中转换到参考坐标系B 中。假设已知随动力，需要通过坐标变换阵将其转换为。

（4）动力学方程

将动能变分、内能变分和外力虚功的具体表达形式代入到式 (11) 中，经分部积分最终可得

二、有限元推导与数值求解法

为求解式 (19)，传统方式是对空间位移矢量和转动进行插值，代入弱形式的积分方程，离散化处理后进行求解。引言中已经提到，由于转动量不可叠加，且会带来“奇异性”和“客观性”问题。本文将推导以应变为插值变量的离散动力学方程，其关键在于建立构形向量与应变量之间的雅可比关系。

由于篇幅原因，此章节内容不再赘叙，感兴趣的朋友可以在仿真秀app公众号对话框回复【孟杨】下载完整PDF。

三、算例分析

1、模型描述

为了探究本文模型的收敛特性和求解精度，将以悬臂梁模型为例进行数值仿真，并将计算结果与商业有限元软件 MSC.Nastran 的求解结果进行比对。在 Nastran 中，非线性静力学问题的求解选用SOL106 求解器，动响应求解选用 SOL400 求解器。这 2 种求解器以非线性位移基有限元为理论基础，处理非线性问题时，在每个迭代步更新切线刚度阵来实现非线性方程的求解。

图 3 为悬臂梁模型示意图，其横截面为矩形，主要物理属性如表 1 所示。后续仿真分析中单元划分均为 20 个，收敛精度均设置为。

2、非线性静力分析

先考虑悬臂梁端部受纯弯矩作用，梁弯曲刚度为，E为弹性模量，I 为截面惯性矩。这种情况下，梁变形的精确解为半径的圆弧，即当时，梁变形成完整的圆，圆的周长等于梁的长度。图 4 为悬臂梁受端部弯矩作用变形图，仿真结果与理论解吻合。在平面问题中，虽然梁呈现大变形特征，但由于弯矩和曲率是线性相关的，采用应变梁模型处理此类问题时相当于是求解线性问题。数值计算表明，仅需要 1 个迭代步即可得到如图 4 所示的结果。

考虑悬臂梁端部受定向载荷Fz作用，方向沿 z轴正向。本文模型与 Nastran 在分析中均将梁划分为 20 个单元。梁端部轴向位移（沿 x 正向）和端部垂向位移（沿 z 正向）随端部受载大小的变化如图 5所示。图中 Nastran 结果为 SOL106 求解器计算的非线性静力学结果。本文模型计算结果与 Nastran结果基本一致。进一步地，对2 种方法计算结果进行量化对比分析。表 2 和表 3 分别为端部定向载荷从 0.6 N增大到 3.0 N 过程中端部垂向位移合轴向位移的计算结果。当端部定向载荷达到 3.0 N 时，端部轴向位移已经超过梁总长度的 40%，呈现出典型的大变形特征。根据表中相对误差数据，不同载荷作用下端部垂向位移的平均相对误差为−0.062%，轴向位移的平均相对误差为−0.045%，且随着载荷增大，相对误差始终保持在这样很小的水平。表明本文模型在处理几何非线性问题时，与 Nastran 中经典的位移基非线性有限元具备相当的计算精度。

以端部载荷为例，进一步分析对比本文模型和 Nastran 中非线性静力学求解的收敛特性。图6为Nastran 在非线性静力学求解迭代收敛曲线。在 Nastran-SOL106 求解器中，设置位移收敛准则的收敛精度为，采用增量加载方式，计算过程中会根据收敛参数的变化自适应调节载荷因子大小。

针对该算例，经 450 个迭代步后达到收敛精度。图 7 为本文模型非线性静力学求解迭代收敛曲线。采用一次性全量加载方式，本文模型仅需10 个迭代步即达到收敛精度。表 4 为本文模型各迭代步端部垂向位移。在最大变形超出梁长度40% 的大变形情况下，本文模型迅速逼近收敛值。对比 Nastran 仿真结果，本文模型在相当求解精度下，收敛特性更好。

考虑悬臂梁端部受载，固定载荷大小为 2 N，分别对定向载荷作用和随动载荷作用 2 种情况进行分析。定向载荷作用方向始终沿 z 轴正向；随动载荷作用方向则始终保持与梁参考线切向垂直。图 8 包含了悬臂梁在定向和随动载荷作用下的梁变形，采用本文模型计算出来的结果与 Nastran仿真结果吻合度较好。由于 Nastran 中本身不具有直接施加随动载荷的功能，在梁端部建立大刚度梁元（其截面刚度超过原始梁刚度的 100 倍），通过FORCE1 卡片施加从端部节点指向大刚度梁元端点的载荷，这样即保证了加载方向始终沿指定方向。这种随动载荷施加方式在一定程度上会引入系统性误差，处理工程实际问题时，需要调整短梁刚度属性，保证不改变原始结构的固有振动特性和变形特性。由图 8 可以看出，定向载荷作用下梁的构型与随动载荷作用具有明显区别。

3、非线性动响应分析

先考虑定向载荷作用。在悬臂梁端部施加动载荷，载荷方向始终沿z向，载荷大小定义为

Nastran 中采用 SOL400 求解器进行动力学求解，时间步长为0.002 s，时间步为2 000 步，求解0～10 s的动响应；应变梁中时间步长同样设置为 0.002 s。图9 为端部垂向位移和轴向位移的响应曲线。前 3 s时间内，两者响应曲线几乎完全重合；随着时间推移，两者在一些波峰和波谷处出现一定偏差，这主要是数值误差积累所致；整体来看，本文模型与Nastran 仿真结果一致性较好。

考虑随动载荷情况，在端部施加与式 (57) 载荷大小一致的动载荷，区别在于载荷方向随着梁变形始终保持与梁中心参考线切向垂直。图 10 为随动正弦载荷作用下位移响应。整体来看，本文模型与Nastran 仿真具有较好的一致性，但相较于定向载荷情况，两者的相对误差有所增大。除了数值误差之外，这种相对误差的增大主要源于前文提到的Nastran 中随动载荷施加方式导致的系统性误差。

四、试验验证

为进一步验证本文模型在工程实际应用中的精度和有效性，开展针对典型大展弦比机翼的地面静力试验。

1、试验模型描述

图 11 为用于地面试验的机翼模型，其主要参数如表 5 所示，其中，机翼翼型为 NACA0015，弹性轴位置在 50% 弦长处。机翼主梁为一根矩形截面钢制梁，其物理属性与表 1 所列悬臂梁一致。沿机翼展向共有 8 个翼盒，由木制翼肋、蒙板和热缩蒙皮制成。翼盒仅中点处与主梁表面胶接，且每两段翼盒间保留 3 mm 间隙，从而避免翼盒带来附加刚度，确保机翼刚度仅由主梁提供。机翼端部布置60g的配重锤，便于在重力作用下使机翼发生大变形。

仿真分析中将机翼等效为梁模型。刚度特性由主梁决定，质量特性由机翼主梁、翼盒和配重锤三部分构成。称重后，将翼盒和配重锤分散为集中质量点进行建模。

图 12 为地面静力试验示意图，机翼根部与基座固连，考察重力作用下机翼变形情况。试验通过2 种方式获取机翼变形：①激光位移计直接测量机翼变形值；②光纤传感设备获取机翼主梁表面应变，再利用文献 [19] 所述应变-位移恢复算法计算静变形。

图 13 为光纤光栅传感器黏贴示意图，传感器为方法②中黏贴于机翼主梁表面的光纤光栅（fiber grating，FBG）传感器。沿主梁上表面平行布置 2 根光纤，分别记为 s1 和 s2，每根光纤串接 5 个FBG 传感器。图中括号内的数字代表传感器位置与根部距离，两根光纤间距 27 mm。

2、试验结果分析

图 14 为机翼主梁表面弯曲应变分布。本文模型仿真分析时将梁划分为 25 个单元。由于单元内部应变为常值，图中应变梁对应的应变曲线呈阶梯状。光纤 s1 和 s2 对应的应变曲线是对离散测点数据样条插值得到的。其中，判定光纤s2 的第 2 个传感器 s2-2 采集的应变值为异常数据，插值拟合时不予考虑。后续计算机翼弯曲变形时，采用 2 根光纤应变数据的平均值作为实测应变值。总体而言，本文模型仿真结果比实测应变值小，但相差不大。

图 15 为自重作用下机翼垂向变形。图中依次绘制了激光位移计、光纤传感测试得到的实测值及本文模型、Nastran 计算得到的仿真值。机翼尖部变形值超过翼展的 30%，呈现典型大变形特征。4 条曲线呈现较好的一致性。为进一步量化分析应变梁的计算精度，沿机翼展向取 5 个站位，评估变形相对误差，如表 6 所示，为了提高试验值的可靠性，表中垂向变形的试验值取激光位移计和光纤传感测试值的平均值。