一种跨声速风洞透气壁模型的数值模拟研究

牟斌，王建涛，黄勇，肖中云，刘大伟

(1.中国空气动力研究与发展中心 计算空气动力研究所， 绵阳 621000) (2.中国空气动力研究与发展中心 高速空气动力研究所， 绵阳 621000)

0 引 言

宽体客机研制对模型风洞试验数据的精准度要求很高，需要消除或修正试验中的各种干扰因素。在跨声速流动试验时，风洞试验段流场中存在亚/跨/超声速的混合流动，有激波、激波/边界层干扰及洞壁波反射干扰，还有试验段开孔或开槽壁的复杂横流流动，流场十分复杂，传统线性修正方法难以应用，故厘清跨声速洞壁干扰有重要意义。

跨声速洞壁干扰影响修正一直是气动试验技术的一项研究主题。20世纪90年代末，范召林等发展了壁压信息法进行洞壁干扰修正，以FL-21风洞中实测的透气壁附近的压力分布作为风洞流场的边界条件，应用Navier-Stokes(N-S)方程，数值模拟风洞中模型的绕流流场，通过和无洞壁自由流场模拟结果的对比，得到洞壁干扰对模型气动力的影响；2016年，李鸿岩等针对小展弦比飞翼标模，利用FL-2风洞实测壁压信息构造开孔壁边界条件，通过求解N-S方程建立跨声速洞壁干扰的非线性修正方法；同年，钟世东等通过小扰动壁压信息法、全速势位流方法和基于RANS的壁压信息法，预测洞壁干扰修正量，以大风洞FL-24的试验结果为标准，对小展弦比飞翼标模确认了这种修正的有效性。壁压信息法适用范围较广，但是需要在试验的同时测量洞壁的静压分布或速度分布，而实际试验中，测压管数量可能偏少或分布不合理，得到的壁面压力采样精准度有限；同时，测压管本身也会带来干扰。

为了弥补壁压信息法的不足，人们转向研究通过建立透气壁壁面边界模型，再进行数值模拟修正洞壁干扰的新方法。日本在透气壁洞壁干扰修正方法的研究上开展了一些计算与试验研究。2006年，J.Akatsuka等改进了Harloff方法，在透气壁壁面边界模型中引入了边界层位移厚度的影响，更好地预测了从驻室流向主流时的质量流量；2010年，T.Nambu等在JAXA的2 m×2 m跨声速风洞的不同参数条件下，利用数值计算方法模拟单孔流动，得出流过孔的质量流量与孔两侧压差的关系，提出了一种新的适用于跨声速风洞开孔壁的壁面流动模型；A.Hashimoto等应用T.Nambu的模型，对ONERA-M5模型在JTWT跨声速风洞中的洞壁干扰和支架干扰进行了数值模拟，计算结果与试验结果吻合很好；应用该孔壁模型，T.Nambu等还对二维翼型试验的跨声速风洞洞壁干扰进行了数值模拟研究，分析了开孔壁洞壁下洗和阻塞效应。

近年来，国内针对透气壁的研究逐渐增加。2015年，曹世坤对单孔在FL-3风洞中的流动进行了数值模拟，根据对开孔壁壁面流动细节的分析，结合有关理论建立了开孔壁的壁面流动模型，用作模拟跨声速风洞洞壁的边界条件；2017年，金佳林等建立了多孔板模型，发展了一种研究跨声速孔壁风洞的孔壁效应的数值方法；同年，国立强应用孔壁模型的数值模拟方法，进行了开孔壁洞壁干扰的修正技术研究。

透气壁的透气特性与其几何形状(例如开闭比、孔径、壁厚等参数)有关，同时还受到试验马赫数、雷诺数、洞壁边界层厚度的影响；实际上，洞壁附近的气流偏角变化也有重要影响，这取决于试验模型的布局构型、缩比大小、飞行姿态以及与模型的距离。因此，建立透气壁壁面边界模型要想尽可能逼真地反映透气特性，必须考虑具体的风洞、具体的试验模型以及具体的试验内容，才可能较好地预测洞壁干扰。

本文针对FL-26风洞的洞壁情况，首先建立单孔流动模型，得到流量与压差的曲线关系；其次将该拟合曲线关系式作为透气壁模型，用作壁面边界条件，模拟宽体客机模型在风洞试验中的绕流；最后为了对比，还应用传统的壁压信息法，以实测壁压作为边界条件进行数值计算。

1 数值模拟方法

1.1 控制方程及其离散

1.1.1 控制方程

控制方程选取任意坐标系下守恒型的雷诺平均N-S方程组(RANS)，应用有限体积法离散。

(1)

式中：

Ω

为控制体；∂

Ω

为控制体边界；为控制体边界外法向单位向量；d

V

是体积微元；d

S

是面积微元；为面上的通量，可采用各种格式离散得到；为守恒变量。

在格心格式中，控制体取为网格单元，守恒变量存储在网格中心，其定义为

=[

ρ

ρu

ρu

ρu

e

]通量张量可分为无黏项和黏性项，具体形式详见文献[17]。

1.1.2 控制方程离散

控制方程无黏项的空间离散采用Roe通量差分格式(flux-difference splitting scheme)，黏性项采用中心差分格式离散。时间离散主要应用交替方向隐式近似分解(AF-ADI)及对称高斯-赛德尔近似LU分解(LU-SGS)。流场收敛加速措施应用局部时间步长、多重网格等方法。计算分别采用S-A一方程湍流模型和

k

-

ω

SST二方程湍流模型。

1.2 边界条件

所用到的边界条件分为基于一维流动特征关系式的无反射远场边界条件、固壁无滑移边界条件、对称边界条件三类。对黏性固壁来说，流体运动速度等于物面的运动速度；黏性计算中，固壁外虚拟两排点后采用统一的格式在全场计算，统一计算界面通量后，将物面无黏通量修正到只有压力贡献上。

1.2.1 试验段上游入流边界条件

入流边界条件需要指定总压、总温。首先计算内点的马赫数及流动角：

(2)

边界面上的压力和密度由式(3)得到：

(3)

温度由状态方程得到，可由此得到声速。结合外插得到的马赫数和流动角，可得边界面上的三个速度分量。对于试验段入口，流动沿

x

向正向。

1.2.2 壁压法透气壁边界条件

应用壁压信息法，壁面压强取风洞试验测压结果。由于风洞测压数据存在跳动，扣除空风洞壁压数据后，还存在一定程度的散布，可能是测压孔加工质量等原因所致。可采用三次样条曲线或“五点三次平滑方法”，把误差光顺掉，再将壁压数据使用样条插值法插值到风洞壁面上。虚拟点速度直接取内点外插。

1.2.3 透气壁模型壁面边界条件

应用透气壁模型时，壁面界面附近法向速度

V

由透气壁模型给出，压强

p

及密度

ρ

均采用内点(下标i)外插：

(4)

其中，切向速度(

u

v

w

)由内点速度矢量减去法向速度分量得到。

1.2.4 单孔流动模拟用的驻室外边界出入流条件

模拟试验段壁板的单孔流动时，设置驻室的计算外边界，必须特别谨慎。小孔没有流动时，驻室流动静止；有流动时，驻室外边界的流动可能“流出”，也可能“流入”。

对于“流出”情况，按照亚声速出流处理即可。

(5)

对于“流入”情况，假定总压为驻室静压，总密度为来流密度，即

(6)

外边界的虚拟点马赫数及流动方向均由内点外插，即

(7)

根据一维等熵流动关系式，采用式(6)～式(7)得到：

(8)

根据式(8)求得的密度、压强得到声速

a

后，得到边界点的速度

V

，进而得到各分量：

(9)

1.3 并行重叠网格技术

网格生成在CFD数值模拟中占据重要地位。宽体客机标模试验需要模拟的部件包括：飞机主体、支杆、弯刀、洞壁等，这些部件之间的距离很近，生成多块对接的网格难度很大。并且风洞洞壁与飞机存在相对运动，对接网格无法满足要求。本文应用重叠网格技术模拟飞机与风洞洞壁流场。

对于静态计算，重叠网格插值及单元信息预先处理并存储，并行计算剖分网格时，根据网格变化信息将重叠网格相关信息进行变换即可，不需要重新“挖洞”及搜索贡献单元。

2 单孔流动建模

2.1 透气壁模型法简介

目前的理论和试验研究均表明开孔壁的边界条件是十分复杂的。开孔壁的透气特性与洞壁几何参数(如开闭比、孔径、壁厚等)有关，还受试验马赫数、雷诺数、洞壁边界层厚度以及洞壁附近的气流偏角变化的影响。洞壁附近的气流偏角取决于试验对象，即与模型的形状、大小、姿态以及距离模型的远近有关。因此，应当针对每一座风洞建立透气壁模型，针对具体的试验计算洞壁干扰。开孔壁上透气孔的数量众多，不可能对每个孔都进行数值模拟，因为不仅网格划分困难而且计算机资源也难以满足要求，建立的物理模型需要进行合理的简化。研究表明，在整个开孔壁上应用单孔流动模型来预测洞壁干扰是可行的。因此本文仅对开孔洞壁上的一个单孔流动进行建模。

对于指定风洞来说，影响开孔流动的主要参数包括：孔两侧压差、来流马赫数、试验段边界层位移厚度。考虑到FL-26风洞试验段距离稳定段出口为11.5 m，边界层已经充分发展，在模型长度变化范围内，边界层的变化不大，因此，在建模中没有考虑边界层位移厚度变化的影响。针对每个马赫数建模：指定孔两侧压差，计算通过孔的流量，用多项式拟合，得到压差与流量关系的拟合公式。

单孔流动建模完成后，在实际应用中，风洞的孔壁作为透气壁处理。以近壁第一层网格单元的压强和驻室压强差为变量，应用单孔模型得到流量，再根据壁面网格面积，得到法向速度，而后利用1.2.3节透气壁模型壁面边界条件处理即可。

2.2 单孔计算网格及边界条件

单孔几何模型的建立主要采用FL-26跨声速风洞的试验段及其开孔壁板的孔尺寸。该风洞是一座半回流、暂冲引射式跨声速风洞，其试验段横截面尺寸为2.4 m×2.4 m。目前有四个试验段：全模试验段、半模试验段、张线试验段和内埋试验段。其中，全模试验段四壁开孔，开孔率为4.3%。孔径为24 mm，左右壁板孔深20 mm，上下壁板定板孔深20 mm，动板孔深15 mm，顺流向60°，稳定段到开孔壁起始距离11.5 m。

FL-26跨声速风洞主要技术性能指标如下：

试验段尺寸：2.4 m(宽)×2.4 m(高)×7 m(长)；

Ma

范围：0.3～1.2，1.4；

总压范围：1.1×10～4.5×10Pa；

Re

范围：

Re

=1.76×10～17.00×10(

c

=0.24 m)；

迎角范围：-22°～22°，-10°～34°，(12°双转轴机构)；

侧滑角范围：-12°～12°(12°双转轴机构)。

比较试验段、驻室和孔径尺寸可以看出：透气孔相对于风洞洞壁、驻室尺寸是很小的。驻室内流动速度极小，可以看做是静止来流。三维计算模型简化后如图1所示，

x

向为流向，

y

向为风洞顶壁方向，

z

向为侧壁方向；在

y

向上带孔的洞壁把计算区域分成驻室和试验段两部分。

(a) 单孔截面图

(b) 单孔立体示意图

初场设置：驻室及小孔内流动速度为0，压强为驻室压强；试验段流场设置为自由来流。试验段来流压强

p

与驻室压强

p

的关系为

(10)

边界条件设置：

x

方向驻室前后以及

y

方向驻室上面，均设置为1.2.4节的驻室出入流条件，驻室下壁面设置为无黏固壁；试验段上壁面为无黏固壁；

x

方向上游设置为自由来流，

x

方向下游设置为亚声速出流；

z

方向上，驻室和试验段的左右面均设置为对称边界条件。

单孔计算网格如图2所示。

(a) 单孔对称面截面网格

(b) 单孔物面网格

网格共14块，网格格子数约200万，网格法向最小间距为2×10，

y

≈0.2。驻室及试验段相对于孔尺寸很大，满足黏性计算要求。单孔计算网格采用“O”型网格结构，消除奇性轴影响，孔附近网格密度、均匀性也满足喷流计算要求。

2.3 单孔流动计算结果

以宽体客机风洞试验模型为对象，研究0.8≤

Ma

≤0.92范围内的洞壁干扰效应。计算状态依据试验选取4个马赫数：0.80，0.85，0.89，0.92。计算主要考虑中小压差的影响，即-1.0<Δ

C

<1.0。

在[-0.1,0.1]区间以间距0.01划分，其余部分以0.1划分区间。通过指定不同的压差，数值计算积分得到通过小孔的质量流量，绘制压差和流量的关系曲线，如图3所示。利用多项式拟合曲线，得到拟合公式，即最后应用到全机计算的透气壁模型。

图3 Ma=0.92流量随压差变化比较曲线

拟合压差与流量关系曲线，采用分段三次多项式曲线，具体形式为

y

(

x

)=

a

+

bx

+

cx

+

dx

(11)

曲线拟合时分三段：[-1,-0.09)，[-0.09,0.09)，[0.09,1]。

从图3可以看出：马赫数0.92情况下，拟合曲线与实际曲线误差不到0.1%，其他马赫数下的拟合曲线精度类似。

不同马赫数下流量随压差变化曲线如图4所示，可以看出：在小压差情况下，质量流量和压差的关系呈现明显的线性关系，随着压差进一步增大，流量增加变缓，非线性现象增强；随着马赫数增大，通过孔的流量增大；在正压差下(流出试验段情况)，随着压差进一步增大，各马赫数下的流量差别减小，在压差为1.0时，各马赫数下的小孔流量趋于同一值，但在负压差(流入试验段情况)下未出现流量趋于同一值的现象。需要注意的是，流量与压差的关系并不关于原点对称，略有偏移，这可能是由于边界层流动的影响而造成的。

图4 不同马赫数下流量随压差变化曲线

马赫数为0.85时，不同压差下马赫数云图和流线图如图5所示，可以看出：对于“流出试验段”情况，小压差和大压差流场结构相似，由于试验段流出气流的引射作用，驻室气流先缓慢向小孔汇聚，后跟随着试验段喷出气流远离小孔；对于负压差“流入试验段”情况，驻室侧气流的压力大于试验段主流，被引射流入小孔，而且压力越大，流入气体越多；在大压差情况下，流入试验段气流与主流交汇处形成分离涡，对试验段壁板边界层会有较大影响。

(a) ΔCp=0.05

(b) ΔCp=0.50

(c) ΔCp=-0.10

(d) ΔCp=-0.50

马赫数为0.92，压差为1.0时的典型流出状态流场图如图6所示，可以看出：在大压差下，试验段主流加速流出小孔，在小孔内部达到超声速，在小孔出口形成斜激波，并在小孔外形成“波节”。

(a) 压力等值线

(b) 马赫数等值线

负压差下的流入状态如图7所示。

(a) 压力等值线

(b) 马赫数等值线

从图7可以看出：吸入气体后，在小孔后会形成一个死水区，减小了流管面积，试验段主流加速。

3 宽体客机洞壁干扰数值模拟

3.1 程序验证

计算外形为宽体客机标模，不考虑风洞洞壁及支杆等。模型基本参数为

参考面积：0.286 62×10mm

力矩参考点：(844.727 mm，0 mm，0 mm)

参考长度：0.202 51×10mm

计算网格取半模，网格共454块，网格点为36 109 294，格子数为33 930 240，网格长宽均取为模型机身长度的50倍以上。模型的整体网格如图8所示，可以看出：外场采用“H”型网格，物面附近采用“O”型网格，生成边界层网格，法向网格最小间距为1.5×10，模型平均气动弦长

C

为202.51 mm，法向最小间距小于10

C

要求。对称面网格在机翼、平尾、头部均加密，整个空间网格都比较均匀。

图8 宽体客机标模的计算网格

计算网格边界条件主要包含物面、对称面、远场。来流马赫数依据试验给定，雷诺数在4个马赫数下变化不大，统一取为

Re

=5.1×10。计算格式选Roe格式，限制器取光滑限制器(Albada型)，三阶插值，湍流模型分别采用S-A和

k

-

ω

SST模型，时间离散选LU-SGS离散，计算中应用三重多重网格加速计算。为了验证本文方法的可靠性，以

Ma

=0.85为典型状态，采用自研软件PMB3D与美国CFL3D V6.4软件进行对比计算。马赫数为0.85、攻角为6.84°的典型流场如图9所示。

(a) 物面压力等值线

(b) 截面马赫数云图和流线

从图9(a)可以看出：压力等值线显示在机翼前缘有比较复杂的激波-边界层干扰，物面流线清晰显示出机翼上出现分离，且分离点从翼根处后缘逐渐发展到前缘。从图9(b)可以看出：激波-边界层的相互作用在机翼中段形成大分离涡。

跨声速流动数值模拟的难点在于攻角超过一定范围会出现激波诱导分离，激波位置及分离点的准确捕捉对于准确计算气动力均非常重要。

计算的升力、阻力、俯仰力矩系数与试验值对比如图10所示。

(a) 升力系数

(b) 阻力系数

(c) 俯仰力矩系数

(d) 升阻比

从图10可以看出：PMB3D与CFL3D计算结果一致性非常好，用S-A模型的结果几乎完全重合，而SST模型在大攻角略有差异。与FL-26包含洞壁、支架干扰的试验值相比，计算升力在3°攻角以前吻合较好，3°攻角以后差别逐渐增大，推测可能是计算中模拟的分离区偏大，导致升力系数显著下降。用SST模型在大攻角时与试验值相差更大，极有可能是SST计算的分离区比S-A模型的分离区更大所致。从阻力来看，用S-A模型计算的阻力在2°攻角后与试验值吻合很好，1°和-2°攻角差别较大，计算低估了阻力。俯仰力矩系数与试验值变化规律一致，拐点接近，尽管斜率差别较大。

验证计算表明，在外形、网格一致的情况下，本文解算器PMB3D与国外知名软件CFL3D具有很好的一致性。与试验测试结果的对比表明，计算的升力、阻力系数与试验规律一致，但对这种具有强激波-边界层相互作用的绕流，在量值上还有较大差异。

3.2 洞壁干扰计算结果

计算模型为全模，周围计算网格如图11所示，共分859块，1～856块为飞机网格，857～859块为风洞洞壁网格，格子数为37 683 200。攻角变化通过旋转飞机网格用重叠网格技术实现。

图11 整体网格结构

应用透气壁模型和壁压信息两种方法模拟洞壁干扰，本文透气壁模型方法得到的洞壁压力云图如图12所示，可以看出：洞壁影响在流向上的不均匀性。

图12 洞壁压力云图(单孔模型方法)

透气壁模型法基于单孔流动的建模结果，数值模拟时在洞壁边界条件上应用；壁压信息法基于风洞洞壁上测得的压力值，应用多项式平滑后再应用“双线性插值方法”插值到数值模拟的壁面，即洞壁边界压力用插值结果，但速度及密度根据流场外插得到。

三种方法计算的升阻特性与试验值的比较如图13～图15所示，“plane”为飞机在自由大气中的计算结果；“plane+WT”为飞机在风洞中采用透气壁模型法的计算结果；“plane+WT2”为飞机在风洞中采用壁压信息法的计算结果；“Exp.”为FL-26风洞带洞壁、支架干扰的试验结果。

(a) 升力系数

(b) 阻力系数

(a) 升力系数

(b) 阻力系数

(a) 升力系数

(b) 阻力系数

从图13～图15可以看出：对于升力系数曲线，用透气壁模型模拟洞壁的计算结果相比自由大气结果有明显改善，在攻角-2°～4°间，升力系数曲线与试验值重合最好，尤其是马赫数为0.89时升力曲线在攻角-2°～6°的与试验重合；用壁压法模拟洞壁的升力系数计算结果也略有改进，但改变量相对较小。对于阻力系数曲线，透气壁模型法在攻角-2°～2°的结果也有改善，但在2°～6°计算值大于试验值；壁压法的计算结果在马赫数0.80时有改进，但在马赫数0.85和0.89时，与自由大气的结果很接近，差量较透气壁模型法要小。另外，计算中发现马赫数0.92时气动力受到支撑干扰很大，还需进一步考虑支撑干扰及洞壁干扰的耦合效应。

4 结 论

(1) 在马赫数0.80/0.85/0.89下，洞壁干扰使得升力、阻力增大，透气壁模型与壁压法定性一致。

(2) 在带洞壁干扰的计算中，采用透气壁模型得到的升力系数与试验值更为接近；阻力在攻角-2°～2°范围内有明显改善，但在2°～6°范围内的计算结果大于试验结果，需要进一步研究。

(3) 与本文壁压信息法相比，模拟跨声速洞壁干扰采用透气壁模型方法得到的修正量较大。