一种采用双弯管储气段布局的高超声速Ludwieg管设计

2022-08-23赵家权司马学昊黄冉冉熊有德

空气动力学学报 2022年4期

赵家权，司马学昊，黄冉冉，熊有德，余涛，吴杰

(华中科技大学，武汉 430074)

0 引言

高超声速空气动力学是先进高超声速飞行器设计的基础，但目前仍有诸多难题尚未完全解决，如高超声速边界层层流-湍流转捩[1-5]、激波-边界层干扰[6]、真实气体效应和多物理场耦合[7]等。为了突破高超声速空气动力学难题，提升先进高超声速飞行器的精细化设计能力，需要持续开展高超声速空气动力学基础研究。通常而言，开展高超声速空气动力学研究主要有三大手段，分别是风洞试验、数值模拟和飞行试验。近年来，高精度数值模拟和飞行试验都取得了较大的进步，风洞试验得益于丰富的地面测试手段，可以以较低成本反映流场物理机理，也得到了科研人员的广泛认可。但是，常规高超声速风洞的建设和运行成本昂贵，测试仪器的价格往往也不菲，一定程度上阻碍了高超声速试验空气动力学的基础研究的发展。为了便于基础科研机构开展高超声速空气动力学试验研究，Ludwieg在1955年首次提出了Ludwieg管的设计概念[8]。Ludwieg管结构简单，无须复杂的压力控制系统，尤其在使用快开主控阀门替代传统膜片后，其建造与运行成本相比常规下吹式高超声速风洞均大幅下降[9]。目前，高超声速Ludwieg管风洞已经得到了广泛应用，如德国不伦瑞克工业大学马赫数3和马赫数6 Ludwieg管风洞[9-10]、美国普渡大学马赫数6 Ludwieg管静音风洞[11]、美国空军实验室马赫数6 Ludwieg管风洞[12]、圣母大学马赫数6 Ludwieg管静音风洞[13]、我国中国空气动力研究与发展中心的马赫数6 Ludwieg管静音风洞等。根据德国不伦瑞克工业大学Radespiel等[9]的估算，一座口径为Φ0.5 m的马赫数6 Ludwieg管风洞整体建造费用约为100万欧元，对应的每车次运行成本约为1欧元。另一方面，受限于Ludwieg管风洞的运行原理，以上提及的高超声速Ludwieg管风洞通常采用直管作为核心部件储气段的布局设计。以运行有效时间为100 ms的Ludwieg管风洞为例，其储气段的长度需要达到17 m以上，对实验室的空间与场地要求较高。又如德国不伦瑞克工业大学、美国普渡大学、圣母大学以及我国中国空气动力研究与发展中心的Ludwieg管风洞均采用长直储气段布局，长度往往达到20 m以上。对于储气段采用弯管气动布局这种形式，尽管Koppenwallner[14]、Schrijer等[15]及美国空军实验室[12]都进行了尝试，但是弯管的参数选择更多来自经验，并未对其具体气动影响开展研究。针对Ludwieg管储气段采用弯管布局产生的气动影响不确定的问题，作者团队进行了基于双弯管储气段布局的高超声速Ludwieg管气动设计，并重点探究了双弯管储气段设计的气动影响。文章首先介绍高超声速Ludwieg管的运行原理；之后对Laval喷管型线与风洞的非定常启动过程进数值分析；而后基于建造的Φ0.25 m马赫数6 Ludwieg管风洞，对运行状态进行初步测量，并将试验结果与数值预测进行对比分析。

1 Ludwieg管理论设计

高超声速Ludwieg管风洞的运行原理如图1所示。风洞的高压储气段与Laval喷管通过快速控制阀门分开。在风洞启动前，储气段内储存着高温高压空气，控制阀门下游的部分则通过真空泵抽成了真空。在开启快速控制阀门的瞬间，会产生一系列的非定常膨胀波，该膨胀波以声速向储气段的上游行进（图1的CD段分别表示该膨胀波的波头和波尾）；该膨胀波驱动管内的气体达到储气段启动马赫数。当膨胀波以当地声速到达储气段尾端后，再次被反射回来（AE和BF）。当反射膨胀波到达快速控制阀时，快速控制阀门关闭，风洞的运行结束。与此同时，在快速控制阀的下游，受压差驱动，气流在Laval喷管喉部形成声速流，并沿着Laval喷管膨胀加速，在试验段时获得对应设计马赫数的高超声速气流。

图1 Ludwieg管风洞的运行原理图[16]Fig. 1 Principle diagram of the Ludwieg tube tunnel[16]

储气段膨胀波前后的总温和总压关系可通过以下关系式获得[17]：

其中：下角标0代表风洞未启动的初始状态，（t,1）代表快速控制阀（快开阀）开启后储气段内膨胀波传播后的总状态；Ma1为Ludwieg管风洞储气段的启动马赫数；γ为理想气体的比热容比。膨胀波系在储气段往返运行一次对应管风洞一个周期，可由式（3）估算：

式中，L为储气段的总长度，a为当地声速。由式（3）可知，通过延长储气段的总长度可以有效增加Ludwieg管风洞的运行时间。为了确保Ludwieg管风洞100 ms的运行时间，储气段的长度需要达到17 m以上。受试验场地限制，本管风洞的储气段采用了双弯管的设计，具体构成包括3段6.5 m长的直管与2个U型接头，其中U型接头中心直径为0.4 m，储气段的内径为0.1 m。采用双弯管设计后，储气段的总长度为22 m。

2 数值方法

数值模拟采用SU2开源求解器[18]。该求解器采用二阶有限体积法，可进行从低速不可压流动到高超声速的计算，满足设计的定常与非定常模拟需求。为了准确模拟高超声速流动，所采用的湍流模式为标准Menter SST两方程湍流模型，空间离散采用AUSM格式。

2.1 网格结构与边界条件设置

为了获得更为准确的结果，Laval喷管建模采用1/4模型，如图2所示。该模型采用全结构化网格生成，壁面第一层网格高度y+为5。1/4模型的左右两侧为对称边界条件，入口和出口分别参考该Ludwieg管风洞实际运行工况，其中入口定义为压力入口边界，出口为压力出口边界。

图2 Laval喷管网格拓扑结构Fig. 2 Mesh topology of the Laval nozzle

由于该Laval喷管为中心轴对称体，网格独立性分析选择了以二维中心轴对称边界的情况，喷管周向的网格分布将参考Wu等之前已经开展过的研究[19]。在二维网格独立性分析中，分别比较了372×45、442×55、497×65、502×75四种不同网格密度下Laval喷管出口处的马赫数分布。由图3可知，网格密度为497×65时Laval喷管出口处的马赫数分布与密度为502×75的分布几乎一致，因此选442×55作为Laval喷管型线优化的参考网格密度。

图3 网格无关性分析Fig. 3 Mesh independent test

此外，为了开展风洞的启动过程以及储气段双弯管设计影响分析，需要进行全风洞的非定常数值模拟。但是采用双弯管储气段布局的模型非中心轴对称体，故采用了风洞的半模进行建模，其网格拓扑如图4所示。模型同样采用全结构化网格生成，总网格量为550万。

2.2 求解器验证

尽管SU2开源求解器具有求解高超声速流场的能力[18]，但是相关工作仍然比较少见。为了验证数值模拟的可靠性，基于德国不伦瑞克工业大学Ludwieg管风洞马赫数6的Laval喷管的型线进行了求解器的验证。在采用同样的数值方法、边界条件和计算网格的前提下，将SU2计算的结果与其他求解器进行了比较。图5显示了分别使用DLR-TAU Code、直接数值模拟[20]和SU2计算的马赫数6的 Laval喷管中心轴上的马赫数分布，三者吻合良好，表明SU2求解器可用于本文高超声速管道流动的模拟。

3 数值结果与分析

3.1 Laval喷管优化设计

Laval喷管是Ludwieg管风洞的主要部件，决定了风洞试验段的流场均匀性。典型的Laval喷管通常由收缩段、喉道以及扩张段组成。为了确保快开阀与收缩段的型面啮合，此次Laval喷管的收缩段设计采用三次样条曲线。在收缩比确定的情况下，开展了不同收缩段长度对喷管出口处流场影响的研究。喷管喉道采用圆弧设计，确保亚声速流动到超声速的均匀过渡。喷管的扩张段则采用了经典的Sivells方法[21]。由于此次Ludwieg管的设计定位是常规噪声型高超声速风洞，喷管的型线设计需要重点满足出口处核心区域马赫数分布的均匀性，而Laval喷管的长度极大程度上决定了风洞的建设成本，因此Laval喷管的设计重点放在了收缩段与扩张段的长度优化上，在流场品质达标的前提下尽量确保缩短Laval喷管的长度。

在收缩段的设计优化中，其长度变化范围为34.4～102.4 mm。为了便于比较，将收缩段长度为34.4 mm的情况定义为参考长度（contraction length of benchmark，CLB），不同收缩段长度下流场马赫数分布比较见图6。

图6 不同收缩段长度下Laval喷管马赫数云图比较Fig. 6 Mach number contour comparison for Laval nozzels with different contraction lengths

由马赫数云图可知，在参考收缩段长度下，虽然Laval喷管存在马赫数分布梯度，但是出口区域已经获得了较为均匀的高超声速流动。当收缩段长度增加到1.5倍参考长度时，Laval喷管的出口区域的均匀度进一步提升，之前存在的马赫数分布梯度消失。当收缩段长度为1.8倍以及2.5倍参考长度时，喷管出口处的马赫数分布保持在均匀分布状态。但是当收缩段长度为3倍参考长度时，喷管出口处的马赫数云图反而呈现不均匀分布。

进一步，比较了Laval喷管中心轴以及出口处的马赫数分布，如图7和图8所示。整体而言，不同收缩段长度下马赫数的分布均相似。但是，将局部马赫数的分布放大后可以区别不同收缩段长度下流场的差异。图7中的局部放大图是显示Laval喷管中马赫数在x/D= 3.0处获得马赫数峰值后下降（注：D表示Laval喷管的出口直径），在喷管核心区域x/D= 3.5～5.9区间呈现振荡分布。当收缩段长度为34.4 mm时，x/D= 3.5～5.9区间流向马赫数的最大相对偏差为0.16%；当收缩段长度增加到1.5倍参考长度时，马赫数的最大相对偏差下降到了0.107%；随着收缩段长度的进一步增加，喷管核心流动区域的马赫数最大相对偏差反而增加到0.13%。喷管出口处的纵向马赫数分布如图8所示。通过图8中的局部放大图可知，喷管收缩段的长度影响了喷管出口处的马赫数幅值。当收缩段长度为34.4 mm时，喷管出口处核心流区域马赫数最大相对偏差为0.034%（注明：分析区域y/D为−0.23～0.23）；当收缩段长度增加到1.5倍参考长度时，马赫数的最大相对偏差下降到了0.061%；随着收缩段长度的进一步增加，喷管出口处核心流动区域马赫数分布最大相对偏差保持在0.054%。综合比较不同收缩段长度下Laval喷管核心区域流向和纵向马赫数分布的最大相对偏差，确定此次Laval喷管收缩段的长度为1.5倍参考长度。

图7 不同收缩段长度下Laval喷管中心轴马赫数分布比较Fig. 7 Comparison of Mach number distributions along the centerline of the Laval nozzle with different contraction lengths

图8 不同收缩段长度下Laval喷管出口处马赫数分布比较Fig. 8 Comparison of Mach number distributions at the exit of the Laval nozzle with different contraction lengths

扩张段作为Laval喷管的关键设计部分，直接影响了高超声速气流的品质。通常喷管长度越长，高超声速流动的膨胀越缓慢，喷管出口处的高超声速流动的均匀性也越好。但是，喷管长度越长，其加工成本越高，型面精度控制难度也越大，风洞试验段的有效口径也更小。因此，Laval喷管扩张段的优化需要在喷管的长度与流场品质上进行折中。选择了对5组不同扩张段长度进行全Laval喷管数值模拟，马赫数云图比较见图9。当扩张段长度为L/D= 6.1时，喷管出口处形成了均匀的高超声速流场。当扩张段长度变为L/D= 4.8时，扩张段的流动加速更快，但是喷管出口处流场仍呈现均匀分布。随着扩张段长度进一步缩短，扩张段的流动加速更加剧烈，在喷管出口处可见明显的马赫数梯度分布。

图9 不同扩张段长度下Laval喷管马赫数云图比较Fig. 9 Comparison of Mach number contours for the Laval nozzels with different expansion lengths

采用同样的分析方法，将不同扩张段长度下Laval喷管的中心轴和出口处马赫数分布进行提取与比较，分析喷管核心流动区域马赫数的最大相对偏差，如图10和图11所示。由图10可知，喷管扩张段长度对流场马赫数的分布有明显影响—喷管越长，高超声速气流加速越平缓，喷管核心流动区域的马赫数分布也更为平坦。通过定量分析发现，当喷管扩张段长度由L/D= 6.1变为L/D= 4.8时，喷管核心流动区域流向马赫数的最大相对偏差由0.15%上升到0.3%。当扩张段的长度进一步缩短后，喷管核心流动区域流向马赫数的最大相对偏差迅速上升到0.64%。

图10 不同扩张段长度下Laval喷管中心轴马赫数分布比较Fig. 10 Comparison of Mach number distributions along the centerline of the Laval nozzle with different expansion lengths

图11 不同扩张段长度下Laval喷管出口处马赫数分布比较Fig. 11 Comparison of Mach number distributions at the exit of the Laval nozzle with different expansion lengths

不同扩张段长度下喷管出口处马赫数分布比较见图11。由该图可知，当喷管扩张段长度为L/D=6.1时，马赫数的分布最为平坦；与此同时，由于喷管的长度最长，喷管出口处的有效直径也最小。随着喷管扩张段长度缩短，喷管出口处马赫数最大相对偏差逐渐增加。在喷管扩张段长度为L/D= 4.1时，喷管出口处马赫数最大相对偏差为0.166%，达到国军标GJB 4399—2002马赫数最大相对偏差要求。综合考虑流场品质以及建造成本，该Laval喷管的扩张段长度选择了L/D= 4.8。

3.2 快开阀对流场影响分析

快开阀主控阀门作为本风洞的设计特点之一，大幅提升了Ludwieg管风洞的运行效率。但是，快开阀位于Laval喷管喉道前，其存在对试验段流场的影响需要明晰。在此分别针对Laval喷管有无快开阀的情况进行了流场模拟，对应的马赫数云图如图12所示。

由图12可见，快开阀首先影响了流动在扩张段的加速收缩，而后受影响的流动经过扩张段膨胀，导致喷管核心区域流向的马赫数分布均匀性有一定程度的降低，但是核心区域马赫数整体仍呈现均匀分布特征。将图12中不同情况下，喷管出口处的马赫数进行提取与比较（图13），可见马赫数的分布几乎一致，但是采用快开阀的情况Laval喷管的边界层厚度增加了约50%，该现象与Kozulovic等之前发现的现象相同[22]。

图12 Laval喷管有无快开阀情况下的马赫数云图对比Fig. 12 Mach number contour comparison for the Laval nozzel with or without the fast-acting valve

图13 喷管有无快开阀条件下出口处马赫数纵向分布比较Fig. 13 Comparison of Mach number distributions at the exit of the nozzel with or without the fast-acting valve

3.3 Ludwieg管风洞启动与运行过程分析

图14显示了Ludwieg管风洞的启动过程，分别给出了风洞启动过程中不同时间点的马赫数云图。由图可知，风洞在2 ms内完成了启动过程，在喷管出口处形成了马赫数6高超声速流动。之后，高超声速流动将持续直至反射的膨胀波达到喷管的收缩段。

图14 Ludwieg管风洞启动过程马赫数云图变化（P0 = 10 bar，T0 = 300 K）Fig. 14 Mach number contour variation during the starting process of the Ludwieg tube tunnel (P0 = 10 bar, T0 = 300 K)

图15 展示了风洞运行过程中不同时间节点膨胀波系在储气段的演化特征。t= 2.23 ms时，Laval喷管的出口处已经形成了高超声速流动，膨胀波系处于快开阀上游，波系整体强度较高；随着风洞的稳定运行，膨胀波系朝向储气段尾部行进，在t= 22.23 ms时到达储气段的第一个U形接头，此时膨胀波系的强度减弱；在t= 68.23 ms时，膨胀波系到达储气段的尾部，发生反射，转而朝向快开阀行进；在t= 121.23 ms时，膨胀波系到达了Laval喷管的收缩段，膨胀波的强度大幅减弱；在t= 127.23 ms时，通过膨胀波波头强度可以得知膨胀波已经发生反射，风洞运行处于第二个运行周期。在Ludwieg管风洞单个运行周期中，Laval喷管段的流动处于稳定状态，而膨胀波波系处于持续减弱状态。

图15 Ludwieg管风洞运行过程密度梯度云图变化（P0 = 10 bar，T0 = 300 K）Fig. 15 Density gradient contour variation during the operation process of the Ludwieg tube tunnel(P0 = 10 bar, T0 = 300 K)

更进一步，分析了Ludwieg管风洞运行过程中快开阀上游位置x= –1.5D处纵向马赫数和压力随时间的分布变化特征，如图16和图17所示。由图16可见，在t= 2 ms时，膨胀波还未达到x= –1.5D处，该处流动尚未启动，尽管Laval喷管已经形成了高超声速气流；在t= 2.23 ms时，膨胀波到达该位置，并形成了启动马赫数0.058。随着风洞的稳定运行，膨胀波系沿着储气段发展，在x= –1.5D处的边界层厚度也由t= 2.23 ms时储气段直径的7%增长到了t= 121.23 ms时储气段直径的67%。值得注意的是，在该过程中由于边界层的厚度持续增长，x= –1.5D处的马赫数始终处于增长趋势。在t= 123.23 ms时，x= –1.5D处的马赫数开始下降，表明Ludwieg管风洞的一个稳定运行周期已结束。

图16 x = –1.5D位置不同时刻马赫数分布比较Fig. 16 Comparison of Mach number distributions at x = –1.5D for different time instances

图17显示了x= –1.5D处不同时刻总压分布的变化特征。由该图可知，随着膨胀波系的运动，总压由10 bar变为了9.25 bar，满足式（2）；此外，在Ludwieg管风洞的一个运行周期内，即从t= 2.23～121.23 ms内，该处的总压几乎保持不变。当一个运行周期结束后，膨胀波系重复之前的运动，形成台阶状的总压分布特征。

图17 x = –1.5D位置不同时刻总压分布比较Fig. 17 Comparison of total pressure distributions at x = –1.5D for different time instances

3.4 弯管对Ludwieg管流场的影响分析

为了认识储气段采用弯管布局对Ludwieg管流场的影响，首先对膨胀波系经过第一个U形接头的过程进行了分析，如图18所示。在t= 22.23 ms时，膨胀波波头快到达U形接头，此时膨胀波为平面波。当膨胀波遇到U形接头后，膨胀波系的运动受离心力作用呈现非对称分布，在U形接头外侧膨胀波强度较高，而内侧强度较弱；即使如此，此处并未观察到明显膨胀波反射现象。当膨胀波系进一步行进，其逐渐恢复为平面波。

图18 膨胀波系经过第一个U形接头的密度梯度的云图Fig. 18 Density gradient contours of the expansion waves passing through the first U-tube joint

图19展示了采用双弯管的储气段中心轴上膨胀波波图，其中Δx-t的斜率表示膨胀波行进的速度。由所绘出的波图可知，风洞在运行过程中膨胀波在弯管部分的斜率并未发生明显变化，更进一步表明采用本弯管设计参数时膨胀波的行进并未受到明显影响。此外，该图还展现了Ludwieg管风洞一个运行周期内储气段里高低压气流间断面位置变化信息，并可以确定风洞一个有效周期约为122 ms。

图19 储气段内膨胀波系的Δx-t波图Fig. 19 Δx-t diagram of the expansion waves in the storage tube

最后，提取储气段x= –1.2D处的壁面压力，其随时间的变化如图20所示，可见风洞的有效运行时间为121.23 ms。进一步将风洞的有效运行时间进行局部放大，发现膨胀波在经过第二个弯管处（t= 43 ms和83 ms）储气段压力存在波动，相对偏差仅为0.04%，可以忽略不计；但是，膨胀波在经过第一个弯管处时储气段压力不存在类似波动。目前该现象的原因尚不清楚，需要进一步开展高精度数值模拟。此外，在单个风洞运行周期内，储气段的压力下降仅为0.64%。

图20 风洞运行过程中储气段流向位置x = −1.2D壁面压力随时间变化Fig. 20 Time variation of the wall pressure at x = −1.2D of the storage tube during the tunnel operation

4 试验验证

基于前文的气动设计，在华中科技大学建成了一座Φ0.25 m的马赫数6 Ludwieg管风洞（简称HUSTHLT），如图21所示。该Ludwieg管风洞由双弯管储气段、快开阀、Laval喷管、试验段、扩张段和真空罐组成。风洞总长约12 m，占地面积约18 m2；风洞每两车次运行间隔为10 min，每天可以运行60次车。

图21 华中科技大学高超声速Ludwieg管风洞Fig. 21 Picture of the hypersonic Ludwieg tube tunnel

在对应图20中x= −1.2D位置安装总压传感器以测量Ludwieg管风洞运行过程中的总压变化。试验段内则采用Pitot探头齐平安装XCQ-062 Kulite传感器。在与数值模拟相同工况运行风洞后，将储气段和试验段的压力数据记录下来，并与非定常数值模拟的数据进行比较，如图22所示，可知风洞的有效运行时间达116 ms，与数值预测结果吻合良好。此外，储气段总压在一个运行周期内的最大压力相对偏差为0.58%。基于Pitot-Rayleigh关系，计算得试验段的来流马赫数为6.05，与风洞设计的名义马赫数相对偏差仅为0.83%。值得留意的是，在数值模拟中观察到的压力波动并未发现，其原因可能是压力传感器的精度有限。

图22 风洞运行过程中储气段压力与试验段Pitot探头压力变化Fig. 22 Time variation of the pressure in the storage tube and the Pitot probe measured pressue in the test section

最后，采用高速纹影系统对试验段的Pitot探头进行脱体激波拍摄，以进一步确定该Ludwieg管风洞的有效运行时间与流场品质。相机的帧速设置为1000帧/秒，每两幅纹影图片之间的时间间隔是1 ms。由图23纹影序列可知，从脱体激波形成到t= 116 ms过程中，Pitot探头的脱体激波稳定，自由来流中无多余的激波与扰动，流场品质较好。