0 引言

飞机在飞行过程中，无时无刻不承受着载荷，准确预测飞机载荷是研制新机的关键技术^[1]，也是飞机结构强度设计的基础。由于真实飞机为连续弹性体，分布载荷的作用必然引起结构变形，使得气动力分布发生变化。且随着飞机设计技术的不断发展，大量复合材料以及大展弦比设计应用到机翼中，考虑飞机的弹性变形对载荷的影响越来越重要。同时，民用飞机适航条款也规定：如果载荷作用下的变形会显著改变外部载荷或内部载荷的分布，则必须考虑载荷分布变化的影响^[2]。可见建立有效的民用飞机静气弹载荷修正方法是十分必要的。另一方面，民用飞机运行场景复杂、工况数量庞大，这对修正方法的计算效率也提出了很高的要求。

截至目前，国内外学者已经建立了多种有效的静气弹载荷修正方法以考虑载荷分布变化的影响。MSC.Flightloads平台^[3]利用线性结构有限元方法考虑弹性变形，采用偶极子格网法计算气动力，通过求解气动弹性线性方程组得到载荷重分布结果。另外，采用计算结构动力学（computational structural dynamics，简称CSD）与计算流体力学（computational fluent dynamics，简称CFD）耦合计算进行静气弹载荷修正也成为了近年的研究趋势^[4]，该耦合方法一般采用模态法^[5]进行结构变形计算，基于动网格技术^[6]进行迭代计算，计算精度较高。但存在计算量大、计算周期长的问题，不能满足大量工况的修正需求。

目前，在工程中应用得更多的是基于NACA-TN3030报告^[7]所发展起来的一种快速修正方法，该方法的理论基础是工程梁法^[8]和片条气动力理论^[8]。它的气动模型建模简单，能与结构模型高效配合，计算速度快，能够满足大量载荷工况的修正要求。刘亚美等人^[9]将此方法应用到某型号飞机分布载荷计算中，结果表明，该方法是一种实用且可靠的线性方法。瞿进^[10]结合飞行载荷设计的特点，对该方法的设计细节进行了完善。郭天天等人^[11]用亚音速线性化理论改进了该方法中的二元升力线气理论，并应用于疲劳载荷计算，扩展了该方法的应用范围。尽管该方法已得到了一定程度的应用，但对其结果的试验验证工作还鲜有报道^[12]，其在解决实际工程问题中的精确程度依然未知。

本文对现有的基于工程梁和片条气动力理论的静气弹载荷修正方法进行了补充和完善，建立了一套完整的设计流程和计算方法，并结合风洞试验数据对该方法进行了验证。

1 修正流程

首先建立静气动弹性载荷修正流程，以机翼的几何外形、刚度以及刚性载荷作为输入，以弹性修正后的分布载荷作为输出。

修正流程主要包括以下4个步骤：

a）刚性载荷计算

刚性载荷的数据来源于风洞试验中刚性模型测压点的压力系数。因为输出的分布载荷将直接转换为有限元节点载荷，所以对网格密度有较高的要求。但测压数据点分布相对稀疏，所以需要先基于气动外形建立致密的面元网格，然后将测压点数据插值到面元节点上。

b）机翼变形计算

计算外载荷引起的机翼变形时，采用一根变刚度的工程梁来描述机翼，机翼变形用梁的变形代替，分为弯曲变形和扭转变形。对于后掠机翼，顺气流剖面的迎角改变量由这两部分变形共同决定。

c）各站位气动力变化计算

机翼变形会引起各站位气动力发生改变，采用马蹄涡影响系数法可计算出各站位气动力增量，由气动力增量计算出各站位由弹性变形引起的迎角改变量，循环迭代可得到各站位收敛后的气动力增量。

d）弹性修正后的分布载荷计算

上述计算步骤获得了弹性修正后各站位的迎角改变量以及气动力增量，该步骤进一步将它们转化为弹性修正后的分布载荷并完成输出。

详细的计算流程如图1所示。

需要说明的是，在各站位气动力分布变化计算过程中，需将计算出的各站位气动力分布变化量作为输入进行下一次迭代计算，以获取各站位附加迎角变化量。根据各站位附加迎角变化量是否为0来判断计算是否收敛，如果还未收敛，将更新各站位迎角，进入下一次迭代运算，直至收敛为止。

2 修正方法

2.1 刚性分布载荷计算

首先，对机翼进行网格划分，将风洞试验数据库中各个剖面的压力系数通过双线性插值的方法插值到网格单元点上，得到每个网格单元的压力系数，如图2所示。已知第i个网格单元中心点O=（x，z），与其相邻的四个风洞数据测压点为Q₁₁=（x₁，z₁）、Q₁₂=（x₁，z₂）、Q₂₁=（x₂，z₁）、Q₂₂=（x₂，z₂），它们对应的风洞测压点的压力系数分别为cp（Q₁₁）、cp（Q₁₂）、cp（Q₂₁）和cp（Q₂₂）。

其次，在X方向进行线性插值得到R₁以及R₂的压力系数cp（R₁）、cp（R₂）；再对Z方向进行线性插值，通过上一步计算出的R₁与R₂在z方向上插值计算出cp（O），则第i个网格单元的压力系数cp（i）=cp（O）。计算公式如下：

通过公式（4），可以计算得到每个网格单元所受气动力，从而计算每个刚轴点所受的累积剪力F、弯矩M和扭矩T。

式中：s（i）为网格单元面积；$\overrightarrow{v}（i）$为网格单元方向向量；q为动压。

如图3所示，P（x′，y′，z′）为任一网格单元的中点坐标，G（x，y，z）为任一刚轴点的坐标，累积载荷从翼尖往翼根累积，通过z′（i）＜z（i）判断共有n个网格单元在此刚轴点的翼尖侧，则此刚轴点的累积载荷为：

式中:f_x（i）、f_y（i）以及f_z（i）为网格单元气动力分别在x、y以及z轴上的分量，求得弯矩和扭矩都是对x轴和z轴取矩，若要对刚轴取矩，需进行以下转化（Λ为机翼后掠角）:

2.2 机翼变形计算

在机翼变形计算中，假设机翼为弹性悬臂梁，通过求解机翼刚轴上的弯矩M和扭矩T来求解顺气流剖面迎角变化量（如图4所示）^[7]：

式中，s为沿刚轴方向自翼根到计算点的距离，Λ为机翼后掠角，EI为抗弯刚度，GJ为抗扭刚度。

将机翼从翼尖到翼根离散成n段，可将弯矩和扭矩写成列向量形式：

由图4中的几何关系可得：

其中，Δx为划分片条宽度。

计算迎角变化量时从翼根到翼尖累积，则可以求得由每个刚轴站位上的刚性分布载荷引起的顺气流剖面的迎角变化{Δα_s0}。机翼的弹性变形会带来机翼载荷变化{ΔM₁}和{ΔT₁}，此部分载荷变化引起{Δα_s1}，以此迭代计算，直到翼尖处顺气流剖面迎角变化Δα_sk（1）=0，最终得到顺气流剖面迎角变化为:

2.3 影响系数计算

首先在机翼上布置如图5所示的马蹄涡，从翼尖至翼根编号1到n。马蹄涡的参考点布置在机翼的1/4弦线处，为满足后缘边界条件在3/4弦线上布置控制点，由这个马蹄涡系构造马蹄涡影响系数矩阵[K]，其元素值K_ij表示第i个控制点受到同侧翼面上（图中左翼）Γ_j以及另一侧翼面相同位置的马蹄涡的共同影响。其中，Γ_j为第j个马蹄涡的强度，由毕奥-萨瓦公式可计算各个涡段对第i个控制点的影响系数^[13]。

2.4 气动分布变化计算

此部分计算中不考虑机翼剖面弯度变形及剖面弦向压心的变化，对机翼剖面的顺流迎角增量采用线性处理。按前文机翼上马蹄涡的布置，在刚性机翼计算中，矩阵方程见公式（15）：

式中:[K]为影响系数矩阵; {Γ}为环量列阵; V₀为自由来流速度; {α}为迎角分布列阵。

根据库塔-儒科夫斯基定理有：

式中:{e}为单位展长的气动分布增量; ρ为空气密度。

联立公式可得：

2.5 弹性修正后的分布载荷计算

通过2.2节方法求得第i个测压剖面弹性修正后的迎角α_t，风洞试验数据库中的数据包含不同迎角下各个测压点的压力系数。通过${\alpha }_{q_1}$＜α_t＜${\alpha }_{q_2}$判断出修正后的迎角α_t在风洞试验数据库两个相邻迎角${\alpha }_{q_1}$和${\alpha }_{q_2}$的范围内，则可以得此剖面第i个测压点弹性修正后的压力系数：

式中: $c{p}_{q_1}$和$c{p}_{q_2}$为对应迎角${\alpha }_{q_1}$和${\alpha }_{q_2}$的压力系数。

最后，通过2.1节的分布载荷计算方法可得到弹性修正后的分布载荷。

3 方法验证

3.1 算例介绍

本文选取某型民用飞机机翼风洞试验模型为分析对象。风洞试验模型包含一副刚性机翼模型和一副弹性机翼模型。对于刚性模型，测量数据包含测压点处的气动压力系数以及翼根的累积载荷（弯矩、剪力和扭矩）。弹性模型的测量数据为翼根累积载荷（弯矩、剪力和扭矩）。

本文选取以下两种典型工况作为算例：工况1的计算条件为飞行马赫数0.6，迎角5°；工况2的计算条件为飞行马赫数0.8，迎角5°。

3.2 弹性载荷修正分析

根据图1的计算流程，需要先按照2.1节介绍的理论方法插值计算刚性气动力，并与风洞试验测量出的刚性模型翼根累积载荷对比。选取低、中、高网格密度（对应网格数量为190、760以及1 710），验证刚性气动力插值计算方法的准确性，这里选取工况1进行对比分析，如表1所示。

由表1的数据可知，使用中密度网格就可以达到满意的计算精度，因此在后续计算中采用中密度网格划分方法，保证计算效率。中密度网格的划分结果如图6所示，图中b为翼展，展向沿x轴，以下用x/b代表机翼展向相对位置。

根据第1节机翼弹性修正流程对机翼的弹性载荷进行计算，得到如图7所示的工况1压力分布云图。由图7可知，弹性修正后上表面前缘的负压区以及下表面的正压区相比于刚性模型有所减小，尤其是在靠近翼尖的位置变化明显。这是由于靠近翼尖位置刚度较小，由气动力造成的结构变形较大，造成分布载荷变化较大。这与试验中观察到的现象一致，说明了修正方法的有效性。

进一步对累积载荷进行求解，弯矩和扭矩是对图3中x轴和z轴取矩，结果如图8所示。在各个站位上弹性修正的累积载荷均小于刚性累积载荷，同时越靠近翼根的位置累积载荷的差值越大。这是由于机翼的气动压心在刚轴之后，使得机翼产生负扭角（如图9所示），从而使机翼上的升力减小，累积剪力同步减小。

载荷在机翼上是从翼尖向翼根逐渐累积的，因此翼根处的弹性与刚性剪力之间的差异最为显著。具体来说，翼根的弹性累积剪力比刚性累积剪力减少了13.15%。

由于机翼的弯曲变形导致压力中心向内移动，升力减小的同时，升力臂的长度也在缩短。因此，弹性累积弯矩相对于刚性弯矩较小，翼根的弹性累积弯矩比刚性弯矩减少了19.80%。

机翼产生的负扭角使得气动压力中心前移，并且升力减小，这导致弹性累积扭矩也比刚性扭矩要小。翼根的弹性累积扭矩比刚性扭矩减少了22.10%。

工况2的计算结果与工况1结果的趋势相同。

表2展示了本文方法对两个算例工况的计算时间（计算机CPU主频为2.6 GHz，运行内存为16 GB），数据表明本文所采用的方法计算速度快，能够满足大量工况的修正需求。使用单台计算机对1万个工况进行修正，耗时在30 min以内。

最后，将计算得到的结果与试验数据进行对比，如表3所示。从表3可以看出，两个工况计算的翼根剪力和弯矩与试验结果对比，误差都在10%以内，翼根扭矩误差稍大，略超10%。且本文计算结果与试验结果相比均偏保守。

这表明本文所建立的修正方法能够在保证计算效率的同时达到工程满意的精度要求。

表3（续）

4 结论

本文充分考虑现代大展弦比民用飞机设计特征和载荷修正要求，对现有的基于工程梁和片条气动力理论的静气弹载荷修正方法进行了补充和完善，建立了一套完整的设计流程和计算方法。所建立的修正方法逻辑清晰、建模过程简单、迭代收敛迅速，适合工程应用；进一步结合弹性模型风洞试验数据对该方法进行了验证。算例对比分析结果表明本文所建立的快速修正方法能够在保证计算效率的同时达到工程满意的精度要求。

参考文献