风洞试验雷诺数/静气动弹性效应分离方法

郭秋亭，孙岩，郭正，刘光远

1. 国防科技大学 空天科学学院，长沙 410073

2. 中国空气动力研究与发展中心 高速空气动力研究所，绵阳 622762

3. 中国空气动力研究与发展中心 计算空气动力研究所，绵阳 621000

一个多世纪以来，风洞试验伴随着飞行器的诞生，在飞行器研制中发挥了推动性的作用，目前仍然是开展飞行器设计、性能评估和流动机理研究的重要手段。高质量的风洞试验数据能够有效降低设计风险、缩短研制周期和挖掘飞行器设计潜能。但受设备物理条件的限制，风洞试验难以对全部的流场相似参数进行模拟。例如，高速风洞试验中，由于试验段尺寸和模型缩比的影响，试验雷诺数通常要比真实飞行雷诺数低1～3个量级[1]。

雷诺数是描述流动尤其是壁面附近边界层内流动演化的重要相似参数[2]。早期由于条件限制，飞行器的风洞试验主要模拟马赫数，较少考虑雷诺数的影响，为此付出过沉重的代价。20世纪50年代，美国空军研制的C-141运输机在设计时没有妥善修正雷诺数对气动性能的影响，焦点位置预判失准，险些酿成重大的飞行事故[3]。

现代大型运输类飞行器为了提高飞行品质、降低油耗和增大航程，广泛采用跨声速性能更优异的超临界翼型。但这类翼型对雷诺数变化十分敏感，真实飞行雷诺数与试验雷诺数下的气动性能差异明显[4]。为了探究飞行雷诺数下的流动机理、推动先进大型运输类飞行器的研制，美国和欧洲先后建成了具有飞行雷诺数模拟能力的低温跨声速风洞，分别为国家跨声速风洞设备[5](National Transonic Facility, NTF)和欧洲跨声速风洞[6](European Transonic Windtunnel, ETW)，能够准确预测雷诺数和静气动弹性影响。

低温高雷诺数风洞主要通过降低试验段气流总温与增大气流总压相结合的方式来提高风洞试验的雷诺数。但增压的同时也增大了速压，增加了风洞模型受到的气动载荷，导致模型出现较大的结构变形，引入了不可忽视的数据误差。数值计算和风洞试验研究结果均表明，模型结构变形与雷诺数变化引起的气动载荷变化量相近，且多数情况下二者的作用相反，使得模型静弹性变形掩盖了雷诺数对气动特性的影响，造成伪雷诺数效应[7-9]。低温高雷诺数风洞的总温、总压可以独立控制，能够保持速压不变，开展不同雷诺数下的风洞试验，单独研究雷诺数对气动特性的影响；也可以保持雷诺数不变，开展不同速压下的风洞试验，单独研究模型静弹性变形对气动特性的影响，从而将雷诺数与静弹性变形的影响进行分离[10]。

国内尚未建成总温、总压可独立控制的风洞设备，高雷诺数试验主要通过增压方式进行[11]，结构静弹性变形影响难以从试验数据分离。目前，对于增压试验中雷诺数/静气弹效应分离问题，尚无有效的解决手段。基于数值模拟或变形测量的模型变形影响修正技术在复杂外形上，仍然存在操作上的复杂性和结果的不确定性[12-14]。

为此，通过分析飞行器气动特性随雷诺数和静气弹变形的变化规律，提出一种基于风洞试验数据的雷诺数/静气弹效应快速分离方法，为风洞增压试验提供一种数据修正技术。

1 数值计算方法

1.1 NNW-FSI软件

模型绕流流场和气动弹性数值模拟均在NNW-FSI软件平台[15]开展。对湍流的模拟采用有限体积方法离散求解雷诺平均奈维尔-斯托克斯方程，该方程中的雷诺应力项通过SST(Shear Stress Transport) 涡粘模型[16]求解。气动/结构耦合数值模拟中结构静变形计算采用柔度矩阵方法，耦合界面数据传递采用薄板样条插值方法，网格运动采用径向基函数-超限插值复合方法实现。

NNW-FSI是国家数值风洞多学科耦合应用软件系统下规划的一款流固耦合模拟软件，在“in-house”求解软件TRIP[17]近20年研究和开发工作的基础上，通过模块化设计、组件封装、代码重构和功能拓展，实现了体系化设计和基础版本的开发。NNW-FSI软件包含前置处理、后置处理、核心解算和运行环境管理4个主要功能模块，能够开展定常/非定常流动、静态气动/结构耦合、颤振、抖振等流动现象的数值仿真，具有较高的可信度[18-19]。

1.2 计算准确性验证

以高雷诺数气动结构(High Reynolds Number Aero-Structural Dynamics,HIRENASD)机翼对静气动弹性耦合计算结果检验数值计算方法的准确性。HIRENASD机翼是气动弹性预测会议(Aeroelastic Predication Workshop, AePW)选择的标准外形，采用典型超临界翼型，气动外形和几何参数见图1[20]。HIRENASD机翼由3段组成，机翼的前缘后掠角为34°、半翼展b/2为1.286 m、机翼参考面积Aref为0.392 6 m2、平均气动弦长cref为0.344 5 m。

图1 HIRENASD机翼几何信息[20]

HIRENASD机翼计算网格采用粗密度多块结构对接网格，如图2所示。流场网格采用H型拓扑结构，机身和机翼周围分别包裹O型网格用于模拟边界层内部流动，机身表面第1层网格高度为6.0×10-7m，机翼表面第1层网格高度为4.4×10-7m，整个网格包含319个网格块、3 158 849 个网格点、3 088 384个网格单元。

图2 HIRENASD模型流场计算网格

图3给出了HIRENASD机翼弯曲变形计算与试验结果[21]的对比曲线。其中,Ma、Re和α分别表示马赫数、雷诺数和迎角,Ma=0.80,Re=1.4×107，α=3.0°;q/E为无量纲载荷因子，定义为速压q与结构弹性模量E的比值，用于表征气动载荷作用下结构变形的大小,q/E=4.7×10-7。图3中TE(Trailing Edge)表示机翼后缘，LE(Leading Edge)表示机翼前缘，Δy表示机翼弯曲位移，η为归一化的机翼展向坐标。结果表明，计算与试验得到的机翼变形沿展向分布从趋势与量值上都比较一致，说明采用的静气动弹性耦合计算方法能够准确预测跨声速绕流流场和模型变形情况。

图3 HIRENASD机翼弯曲变形

2 雷诺数/静气动弹性效应

在验证了数值模拟方法准确性的基础上，仍以1.2节的HIRENASD机翼为研究对象，进一步分析雷诺数/静气动弹性变形对模型气动系数的影响规律，为分离方法的建立提供理论支撑。

2.1 有限元模型与计算条件

图4(a)给出了HIRENASD机翼的结构有限元模型，共包含193 457个网格点、103 458个四面体网格单元，固支约束添加在与天平连接的端面。为了降低结构静变形计算中的柔度矩阵规模，提高耦合数据传递的效率与鲁棒性，选择有限元模型表面359个网格点构建简化柔度矩阵，如图4(b)所示。

图4 HIRENASD机翼有限元模型

流场与耦合计算条件依据总温不可调节风洞的试验条件给定，即雷诺数调节只能通过总压变化实现，雷诺数Re与载荷因子q/E之间为线性变化的关系。

计算条件采用3组不同的雷诺数Re，对应不同的载荷因子q/E，详细的计算条件参见表1。

表1 HIRENASD机翼计算条件

2.2 雷诺数/静气动弹性效应

图5给出了HIRENASD机翼气动力系数变化量ΔCL、ΔCD随迎角的变化曲线，其中的气动力系数变化量均是以Re为0.7×107时刚性外形的气动力系数为基准。

从图5可以看出，刚性机翼(Rigid)随着雷诺数Re的增加，升力系数CL增大，阻力系数CD减小。对于HIRENASD机翼所采用的超临界翼型，雷诺数Re增加会减小边界层厚度，使得机翼的等效弯度增大，从而升力系数增大；而Re增加会减小机翼受到的摩擦阻力，从而使阻力系数CD降低。但当迎角大于一定值，机翼上翼面开始出现流动分离，雷诺数Re对边界层流动的影响减弱，引起的升力系数变化量ΔCL也随着迎角增加开始减小，同时机翼的阻力主要由分离流动引起的压差阻力贡献，摩擦阻力在阻力中的占比已经很小，表现为迎角较大情况下Re对阻力系数CD的影响非常小。

图5 HIRENASD机翼气动力系数变化量

而对于弹性机翼(Elastic)，气动力系数随雷诺数Re的变化与刚性状态下表现出完全迥异的情况。静气动弹性变形对升力系数CL的影响与雷诺数Re的影响相反，使得升力系数降低，且两者的影响量级接近，因此静气动弹性变形可能完全掩盖住Re对升力系数CL的真实影响，造成虚假的雷诺数效应。而对于阻力系数CD，静气动弹性变形对CD的影响与Re相同，使得阻力系数减小，但与雷诺数Re不同的是，静气动弹性变形主要影响的是机翼的诱导阻力，且阻力系数变化量随着迎角增加急剧增大，机翼静气动弹性变形将放大雷诺数Re对阻力系数CD的影响。

2.3 气动力系数随参数的变化规律

图6给出了几个特定迎角下刚性机翼气动力系数变化量ΔCL和ΔCD随lg(Re)的变化曲线，刚性机翼气动力系数变化量为相对于Re为0.7×107时刚性机翼气动力系数的差量。这里选择气动力系数变化量是因为变化量与绝对量之间只相差一个偏移量，二者随lg(Re)的变化规律相同，采用变化量可以将不同迎角下的曲线绘制在同一幅图内，并不影响对刚性机翼气动力系数随参数变化规律的分析。从图6可以看出，对于HIRENASD机翼，在当前的计算状态下，刚性机翼气动力系数或气动力系数变化量与雷诺数的对数lg(Re)之间表现为强线性变化关系。

图6 刚性机翼气动力系数变化量随lg(Re)的变化

图7实线部分给出了迎角α=-1.0°,1.0°,3.0°,5.0°时，弹性机翼气动力系数变化量ΔCf随载荷因子q/E的变化曲线，这里的气动力系数变化量是指机翼静气动弹性变形引起的气动力系数改变量，为每一个对应雷诺数Re和载荷因子q/E下机翼弹性外形与刚性外形气动力系数的差值。可以看出，弹性变形引起的气动力系数变化量与载荷因子之间也呈现出强烈的线性变化关系。

图7 弹性机翼气动力系数变化量随载荷因子q/E的变化

3 分离方法

在雷诺数效应风洞试验中，不考虑其它参数的变化，弹性外形下的气动力系数Cf,elastic可以表示为雷诺数Re和速压q的函数：

Cf,elastic=f(Re,q)

(1)

速压会引起机翼结构的静气动弹性变形，因此可以将气动力系数分成刚性外形气动力系数Cf,rigid与机翼静气动弹性变形引起气动力系数变化量ΔCf,elastic2个部分，其中Cf,rigid仅随Re变化，ΔCf,elastic随Re和q变化，关系式为

Cf,elastic(Re,q)=Cf,rigid(Re)+ΔCf,elastic(Re,q)

(2)

为了将雷诺数Re与速压q引起的静气动弹性变形效应分离，提出2个假设：

假设1模型静气动弹性变形引起的气动力系数变化量ΔCf,elastic随载荷因子q/E呈线性变化关系。

假设2雷诺数Re变化引起ΔCf,elastic的改变量很小，可以忽略不计。

对于假设1，机翼结构在线性小变形范围内，如果不考虑变形前后气动载荷分布的变化，静气动弹性变形引起的气动力系数变化量ΔCf,elastic正比于机翼的结构静变形量，而变形量正比于结构受到的气动载荷、反比于结构材料的弹性模量E，而气动载荷又正比于速度压力q，因此推断ΔCf,elastic与q/E表现为线性比例关系，2.3节中的图7也印证了这个观点。而当载荷因子q/E等于零时，机翼不受到气动载荷的作用，静气动弹性变形为零，因此ΔCf,elastic也等于零。

对于假设2，从雷诺数Re对升力系数的影响可以看到，Re增大会使升力增加，一般情况下增加量相对于绝对量是小量，而机翼静气动弹性变形主要由机翼升力载荷引起，因此Re增大对静气动弹性变形的影响很小，由此推断Re变化对ΔCf,elastic的影响量非常小。

根据假设1与假设2，式(2)可以变化为

Cf,elastic(Re,λ)=Cf,rigid(Re)+k2λ

(3)

式中：λ=107×q/E；k2为待确定的系数。由式(3) 可以看出，雷诺数/静气动弹性效应分离的关键是确定系数k2，进而确定ΔCf,elastic。但通过式(3)无法直接确定比例系数k2，因此提出第3个假设，通过给定Cf,rigid随雷诺数Re的函数关系确定系数k2。

假设3刚性外形气动力系数Cf,rigid随雷诺数的对数lg(Re)呈线性变化关系。

飞行器气动力特性的雷诺数效应研究结果表明[22]，对于给定的飞行器外形和绕流条件，存在一个雷诺数的自准区，在该范围内，气动力特性与雷诺数的对数呈现近似线性变化的关系，第2节的数值模拟结果(图6)也十分吻合这种规律。风洞试验中，为了外插获取飞行雷诺数下的气动数据，试验条件通常会根据经验选择自准区范围内的雷诺数。假设3用于确定式(3)中的系数k2，但在有些情况下假设3可能并不严格满足，后文将对不满足情况的影响做进一步详细的分析。

基于假设3，式(3)可以变换为

Cf,elastic(ξ,λ)=k1ξ+k2λ+k3

(4)

式中：ξ=lg(Re/Re0)；Re0=106；k1、k3为待确定的系数。利用3组不同雷诺数Re与载荷因子q/E组合下的气动力试验数据即可求解线性方程组确定系数k1、k2、k3为

(5)

当不同雷诺数Re与载荷因子q/E的气动力数据超过3组，可以利用最小二乘方法计算待定系数k1、k2、k3。

待定系数k1、k2、k3确定后，利用式(4)即可得到气动力系数随雷诺数与载荷因子的变化特性，其中式(4)右端第1项表示随Re变化的气动力系数，右端第2项表示模型静气动弹性变形引起的气动力系数变化量，第3项k3表示雷诺数Re等于Re0时的刚性气动力系数，即

(6)

4 分离方法测试

采用大展弦比翼/身/平尾组合体风洞模型的试验数据对同步分离方法进行测试，并通过不同载荷因子外插得到的刚性外形风洞试验数据对分离方法获得的结果进行验证和评估。

4.1 测试对象及数据

测试对象选择NASA设计的通用研究模型(Common Research Model,CRM)的翼/身/平尾组合体构型CRM-WBT0，外形如图8所示。CRM模型代表了一类典型的现代宽体客机外形，其设计巡航马赫数Ma为0.85，巡航升力系数为0.50，是第4～6届阻力预测会议(Drag Prediction Workshop,DPW)选择的标准外形。CRM模型的详细几何参数见文献[23]。

图8 CRM-WBT0模型外形

CRM模型先后在Langley研究中心的NTF风洞、Ames研究中心的11 ft (1 ft=0.305 m)跨声速风洞以及欧洲跨声速风洞ETW开展过风洞试验。其中，CRM模型在ETW开展的风洞试验得到了欧洲战略风洞研究潜能改进项目(European Strategic Wind Tunnels Improved Research Potential,ESWIRP)的资助，获得了丰富的试验数据。因此，采用CRM-WBT0模型在互联网上公开发布的ETW风洞试验数据开展分离方法的测试与验证。

表2给出了测试使用的3组试验数据条件，可以利用车次号在ETW发布的CRM模型试验数据中获取对应的气动力系数。

表2 试验参数

此外，为了对气动力系数分离效果进行评估，还采用了1组与表2中第2组具有不同载荷因子的试验数据，用于修正得到刚性外形的气动力系数数据。这组试验数据的载荷因子q/E为5.058×10-7，其他参数条件与第2组相同，对应的ETW试验的车次号为213。

4.2 气动力系数分离结果

图9给出了CRM-WBT0模型的气动力系数分离结果。图中、ETW Rigid表示利用相同雷诺数不同载荷因子下试验数据外插得到的刚性外形气动力数据;ETW Elastic为ETW发布的风洞试验数据;Modeling为采用分离方法预测的刚性外形气动力系数，通过对比Modeling表征的气动力预测数据与ETW Rigid表征的修正刚性外形气动力数据，评估预测结果的准确性。具体的插值思路是：认为在载荷因子q/E为零时，机翼不受气动载荷的作用，即机翼没有变形，利用不同载荷因子数据外插得到q/E等于零的试验数据，即对应刚性外形下的气动力。对于其他雷诺数，前文分析了雷诺数对静气动弹性变形引起的气动力系数变化量影响很小，可以利用给定雷诺数下的弹性与刚性外形气动力系数差量修正得到其他雷诺数下的刚性外形气动力数据。这种基于真实的风洞试验数据修正得到刚性外形气动力系数，用于评估与验证分离方法的结果具有较高的可信度。

可以看出，在模型迎角α≤4.0°时，采用分离方法得到的刚性外形气动力系数与试验修正数据吻合较好，仅阻力系数预测趋势存在小量偏差。原因主要有2个方面。一是模型的迎角较小，机翼受载荷作用后产生的结构静弹性变形对阻力系数的影响非常小，分离方法预测偏差超出了机翼结构静弹性变形对阻力系数的影响量，图9(a)中ETW Rigid和ETW Elastic 2组阻力系数之间的最大差量仅为1.07×10-4。二是分离方法建立在3个前提假设之上，对不满足这些条件的预测结果存在一定的偏差，如图9(b)～图9(d)所示。此外，在迎角α=5.0°时，流动出现分离，原有的理论假设不再满足，升力系数CL与俯仰力矩系数Cm分离结果与ETW Rigid之间存在明显的偏差。

图9 CRM-WBT0模型气动力系数分离结果

4.3 分离方法的分析与讨论

在3个假设完全成立的前提下，建立的分离方法有着严格的数学推导，能够从含有模型静气动弹性变形试验数据中分离得到理想的刚性外形气动力结果。假设1与假设2有着理论上的依据，并不依赖于模型外形与绕流条件。但假设3的满足条件与绕流条件相关，并不是完全成立的。如图9(e)中，此时迎角较大流动出现了分离现象，升力系数CL与俯仰力矩系数Cm随雷诺数的对数并不是线性变化的，从而使得待定系数k1、k2、k3的计算出现不准确，容易造成分离方法获得的结果与刚性试验数据修正结果之间出现较大的偏离。

式(3)中，Cf,rigid随雷诺数变化的函数关系作为系数k2的定解条件，当Cf,rigid随lg(Re)满足线性变化时，利用式(5)能够得到合理的k2值，从而得到合理的刚性外形气动力系数分离结果，如图9(a)～图9(d)所示；而当Cf,rigid随lg(Re)不满足线性变化时，采用式(5)得到分离结果偏差较大，如图9(e)所示。

从上述分析可以看出，依据式(5)的分离方法能够对Cf,rigid随lg(Re)呈现强线性变化的情况得到很好的分离结果。而对于其他情况，可以将分离方法进行变换，由式(3)可以得到：

Cf,rigid(Re)=Cf,elastic(Re,λ)-k2λ

(7)

式(7)中，首先利用式(5)得到一个k2的初值，然后在初值附近变换不同的k2值，得到不同的Cf,rigid随Re的变化数据，并以此分析Cf,rigid随Re的变化规律。此外，该分离方法也可以与其他技术一起，如CFD/CSD耦合数值模拟，共同评估静气动弹性变形影响的系数k2，然后分离得到无弹性变形影响的雷诺数效应试验数据。

5 结 论

利用静气动弹性耦合数值模拟方法研究了雷诺数与静气动弹性变形对气动力系数的影响规律，基于线性理论提出了3个假设，并发展了一种雷诺数/静气动弹性变形效应分离方法，通过大展弦比风洞模型的试验数据对分离方法的准确性进行了检验。结果表明提出的雷诺数/静气动弹性变形效应分离方法能够为常规风洞的雷诺数效应试验提供一种简便快捷的数据分析技术，预估雷诺数与静气动弹性变形对气动力特性的影响，且在刚性模型气动力系数随雷诺数对数呈现近似线性变化的情况下能够得到合理的气动力系数分离结果。