磁浮飞行风洞动模型气动结构耦合仿真评估

高兴龙，王超，符澄,*，孙运强，寇杰

1.空气动力学国家重点实验室，绵阳 621000

2.中国空气动力研究与发展中心 设备设计与测试技术研究所，绵阳 621000

0 引 言

在航空航天领域，各类先进军/民用航空航天飞行器的发展必然给国防安全和社会经济带来巨大变化，同时也带来了大量新的空气动力学问题[1-2]，如高速气动力精确预测及一体化设计、湍流减阻、边界层转捩、地面效应等基础和前沿问题，这些问题的研究都强烈依赖风洞试验，且对风洞流场动态特性、低扰动特性及特种模拟能力要求很高，传统风洞设备越来越难以满足日益特殊的空气动力试验需求[3-4]。在先进轨道交通领域，磁浮真空管道超高速列车等先进轨道交通技术被纳入我国“十三五”重点研发计划。真空管道列车面临着复杂的激波边界层干扰、超声速气动减阻、管道激波反射、活塞效应等一系列复杂的空气动力学问题，必须在风洞中加以解决[5]。目前，传统跨超声速风洞存在试验段尺寸（长度和截面积）和气流相对运动带来的模拟真实性不足等问题，难以在受限空间内开展航空航天飞行器和超高速列车空气动力特性及力/热/结构/控制耦合问题研究[6-9]。

“飞行风洞”的概念最早由美国研究者于20 世纪90 年代提出，具备“体动风静”的特殊运行方式和性能优势。NASA 的兰利研究中心在较早前就开始关注这一全新概念的地面试验设备，称之为“高升力飞行风洞（high-lift flight tunnel）”[10]，并开展了3 期关键技术研究和论证规划工作。由于该设备的关键技术攻关难度较高，至今仍未取得较明显的进展。“磁浮飞行风洞”是利用真空管道列车概念结合动模型试验技术提出的一种新概念风洞设备，其原理是在一段封闭的直线长管道内安装磁浮模型驱动机构，利用电磁悬浮、牵引和导向技术驱动模型高速运动，模拟各类飞行器及高速列车运动的物理过程，构建接近真实飞行环境和运动特点的“体动风静”试验状态，可以满足航空航天飞行器、超高速列车等宽马赫数范围、宽雷诺数范围、低噪声、低湍流度、高真空度（高空）、特殊气体介质、受限空间条件下的空气动力学及其交叉学科地面试验需求[3,5]。

磁浮飞行风洞采用高速磁浮技术驱动试验模型在等截面、直线、密闭管道内做加速、匀速和减速运动，可以模拟高加速/高减速（高过载）过程及速度突变时的气动现象。试验过程中，风洞内模型高速运动所产生的气动力和声波的传播十分复杂，同时，受扰动的流体与结构之间还存在相互作用，属于气动结构交叉学科之间的耦合问题[11]。本文提到的气动结构耦合仿真，与传统流固耦合存在一定区别。传统流固耦合概念为固体力学研究领域所提出，更侧重于结构的动态响应问题，而本文则更侧重于管道内模型高加速/高减速和匀速运行过程的非定常气动特性求解，并不存在结构变形问题；但同时也必须考虑动模型的高加速/高减速情况与流场之间的耦合作用，传统CFD 仿真对于高加速/高减速的流场动态响应求解还存在一定局限。此外，本文研究的方法更多是为后期更深入地研究高过载条件下磁浮平台的结构振动和运动部件的功能及安全性评估等问题提供理论基础。传统CFD 仿真方法在此类动模型高速运动耦合问题的处理上还存在一定局限性，须结合动网格等技术，利用新型气动结构耦合仿真方法对动模型试验过程进行仿真，从气动结构耦合的角度研究模型在密闭管道内高速运动过程的气动特性问题，为磁浮飞行风洞流动控制系统的设计提供技术支撑。

时空守恒元和解元方法（Space Time Conservation Element and Solution Element Method,CE/SE）是近年来兴起的一种全新的高分辨率守恒型方程数值积分方法，由Chang[12]于1995 年提出。该方法将时间和空间统一起来同等对待，利用守恒型积分方程，通过定义守恒元（Conservation Element,CE）和解元（Solution Element,SE）使局部和整体都严格满足守恒律[13]。经过开发，该方法已成为一种新型的可求解二维和三维非稳态欧拉方程的可压缩体求解方法，并被应用于冲击波、气动噪声、磁流体力学等典型多物理场耦合问题[14-19]。

与传统“体静风动”风洞试验技术相比，“体动风静”风洞试验技术主要基于动模型试验技术原理，依靠动力驱动模型在滑轨上运动，从而开展相关的空气动力学测试[3]。国内外在高速动模型试验技术与计算模拟方面开展了大量的研究。日本Kobayasi 物理研究所建立了小型高速列车压力波测量装置，德国宇航中心的哥廷根隧道模拟试验平台可以开展静态和瞬态列车空气动力学研究，韩国建筑技术研究院建立了电动机带传动弹射的小型列车动模型试验平台。国内开展动模型试验技术研究的主要有中南大学轨道交通安全重点实验室、中国科学院先进轨道交通力学研究中心和西南交通大学试验中心等。在气动结构计算方面，中南大学王乐卿和高广军等[3]通过仿真研究了高速列车车身风阻制动板气动外形设计问题；西南交通大学李一凡等[4]利用CFD方法研究了导流装置对高速磁浮列车气动特性的影响。

本文基于CE/SE 方法和气动结构分区耦合技术，根据磁浮飞行风洞的实际试验需求，合理建立仿真工况，从气动结构耦合的角度对动模型在磁浮飞行风洞内的加速、匀速和减速全过程进行气动特性仿真评估，分析动模型在密闭管道内高速运动的流固耦合问题，获得受扰动流场各参数的时间历程曲线，为磁浮飞行风洞流动控制系统和消波措施等的方案设计提供支撑。

1 问题描述与方法

对磁浮飞行风洞的磁浮平台和轨道模型进行适当简化，计算运行过程中平台及模型所受到的气动力。简化后的管道、轨道横截面和平台的模型如图1 所示，其中，z 为模型运动方向，x 为管道横截面的水平方向，y 为管道横截面垂直方向；管道截面可以沿轴向直接拉伸形成三维管道模型，管道的轴向总长度为磁浮飞行风洞的长度。

图1 磁浮飞行风洞磁浮平台及管道截面组合模型图Fig.1 Combined model between maglev platform and tunnel section of maglev flight tunnel

1.1 流体求解方法

首先给出一维的波传播方程：

式中：a ＞0，为对流项常数；u 为波传播速度；t 为时间。该方程令单元空间离散，在不同空间位置和给定时间产生离散点。CE/SE 方法通过构建二维欧拉空间域 E2将时间作为附加空间坐标考虑，在 E2内的时间和空间上应用高斯散度定理，从而将微分形式的式（1）改写为积分形式[13]：

式中：S(V)为 E2内任意时间-空间域 V 的边界；=(au,u)；=dσ，dσ和分别为S(V)表面单元外法向的面积和单位向量。在 E2内建立一个单元体积，该体积内的单元在时间和空间上守恒且被统一处理，这是构建守恒元CE 的关键。之后建立SE，保证SE 单元内的变量足够小，从而可以在 E2内坐标上某点中心的时空区域内部边界附近对流体变量进行近似泰勒级数展开：

式中：u*为SE 单元内部流体变量；uq为时空区域守恒元q 点的流体变量；x 为空间坐标，xq为q 点的空间坐标，tq为q 点的时间坐标。通过流体的对流扩散方程将时间和空间导数相关，从而仅存在uq(x,t)和它的空间导数这两个未知量。

为了令系统封闭，定义2 个CE 的方程，之后获取积分形式的对流方程，如公式（2）所示，则CE/SE格式可以保证SE 在时间和空间上统一守恒。为保证流体守恒，沿CE 所形成的直线域进行时间和空间积分：

1.2 气动结构耦合平台

对于气动结构耦合问题，须将流体求解器和结构求解器进行一致求解，寻找信息传递变量，建立耦合平台。这里可以将CE/SE 方法和结构动力学方法进行耦合，以实现高速运动过程的气动结构耦合计算仿真。上述2 个方法可以独立进行计算，因此可以分别对流场和结构进行网格划分和时间步长设置，耦合系统自动追踪2 个域中的最小时间步长进行计算。气动结构耦合须将拉格朗日结构嵌入到流体域中，并在每个时间步内将节点位移和速度信息作为向流体求解器传递的变量。在初始求解时间步内，CE/SE 方法首先用初始速度和压力求解可压缩流场，获得界面压力；之后将压力作为边界条件施加到结构域上，对结构域进行求解后，得到节点的位移和速度；再将边界节点上的位移和速度返回到流体场，从而使流体和结构进入收敛迭代的过程；当满足收敛条件后，再进入下一个时间步的计算。整个系统的流固耦合迭代过程如图2 和3 所示。

图2 流固耦合平台Fig.2 Fluid structure interaction platform

图3 基于CE/SE 算法的流固耦合求解流程图[20]Fig.3 Solving flowchart of FSI by CE/SE algorithm[20]

1.3 仿真模型及参数设置

结构有限元模型主要由磁浮飞行风洞磁浮平台和试验模型组成。流体模型如图4 所示。流体网格域的长度为1 000 m。根据CE/SE 方法的特点，流场采用结构化网格，最小网格尺寸为0.01 m，结构和流体网格共计约1 200 万单元。结构域和流体域网格分别如图5 和6 所示。

图4 有限元流体模型示意图Fig.4 Layout of finite element flow model

图5 磁浮平台和试验模型结构域网格图Fig.5 Structure mesh of maglev platform and test model

图6 管道横截面流体域网格图Fig.6 Flow mesh of tunnel section

仿真模型采用的长度单位为m，时间单位为s，质量单位为kg。平台和结构模型假设为不变形，材料采用刚体材料。

2 仿真结果及评估分析

2.1 仿真工况

磁浮飞行风洞管道总长度为 1000 m，分为加速段、匀速段和减速段，根据试验条件建立2 个仿真工况：工况1 的最大速度为170 m/s，工况2 的最大速度为340 m/s，加速段和减速段为匀加速和匀减速。工况参数如表1 所示，所有的仿真工况均考虑重力因素。仿真模型模拟的试验阻塞比为20%，即模型运动过程中阻力面积占管道横截面可用面积的20%。

表1 仿真工况参数Table 1 List of simulation conditions

2.2 气动结果分析评估

工况1、2 模型运行速度随时间变化曲线如图7、8 所示。取风洞管道内部纵向位置为0、450、500、600 和1 000 m 的5 个横截面，输出密闭管道内横截面所受压力随时间变化的曲线，工况1 各横截面的结果如图9 所示。从450、500 和600 m 横截面的情况可以看出：这3 个横截面均经历了2 个压力波峰值，第一个压力波峰出现的时刻略早于模型到达所选横截面的时刻，说明模型头部压缩波的运动速度明显快于模型本身；第一个压力波达到峰值后马上急剧下降，之后缓慢增大，这应该是模型头部产生的压缩波运动经过所选横截面后，尾部产生的膨胀波又出现在所选横截面导致的；第二个压力波峰出现时，模型处于减速阶段，可能是压缩波遇到管道末端壁面产生的反射波重新运动到所选横截面位置导致的。1 000 m 处管道末端横截面只产生了一个波峰，在模型减速运动阶段出现。

图7 工况1 模型运行速度随时间变化曲线Fig.7 Model moving velocity vs.time for condition 1

图8 工况2 模型运行速度随时间变化曲线Fig.8 Model moving velocity vs.time for condition 2

图9 工况1 管道内各横截面压力变化曲线Fig.9 Pressure vs.time curves of tunnel sections for condition 1

工况2 各横截面的结果如图10 所示。从图10可以看出：450、500 和600 m 横截面处的情况与工况1 类似，均出现了了模型前端的压缩波和反射波导致的2 次峰值。但工况2 管道末端1 000 m 横截面处出现了多个峰值，应该是压缩波在末端壁面多次反射的结果。

图10 工况2 管道内各横截面的压力变化曲线Fig.10 Pressure vs.time curves of tunnel sections for condition 2

将2 个工况压力随时间变化的曲线与模型运动的速度变化曲线进行对比，可以看出：模型匀速运行过程中，试验段横截面上的压力波均出现了波动峰值，不利于气动试验的开展，须采取一定的消波措施削弱压力波。

图11 为工况1 模型运动过程中所受气动阻力随时间变化的曲线。从图中可以看出：模型所受气动阻力在加速段与匀速段相差不大，在减速段达到最大，最大阻力约为初始阻力（t=0 时）的1.1 倍。图12 为工况2 模型运动过程中所受气动阻力随时间的变化曲线。从图中可以看出：在加速段接近匀速段时，模型所受的气动阻力开始逐渐减小，在减速段，气动阻力达到最大，气动阻力峰值约为加速段所受气动阻力均值的2 倍。

图11 工况1 模型运动过程气动阻力随时间变化曲线Fig.11 Aerodynamic drag vs.time curve during model moving for condition 1

图12 工况2 模型运动过程气动阻力随时间变化曲线Fig.12 Aerodynamic drag vs.time curve during model moving for condition 2

2.3 管道内波系传播结果

图13 和14 为模型在直线密闭管道内加速、匀速和减速运动过程中周围非定常流场演化的三维速度和Schlieren 数云图。从图中可以看出：波阵面的前移速度明显快于模型的运动速度，在模型未到达管道末端时，波阵面已先于模型到达管道末端壁面并产生反射。

图13 模型运动不同时刻速度云图Fig.13 Velocity contours for model moving at different times

图14 模型运动不同时刻Schlieren 数分布云图Fig.14 Schlieren contours for model moving at different times

2.4 管道末端压力分布

为开展消波措施设计，通过仿真计算获得管道末端的压力波变化情况。在管道模型的末端截面不同位置选取若干点（图15），并输出这些点的流场参数时间历程曲线，如图16 所示。

图15 管道末端截面选取单元位置示意图Fig.15 Layout for choosing elements’ positions at end wall

图16 工况1 管道末端截面选取单元的压力变化曲线Fig.16 Pressure vs.time of choosing elements for condition 1

表2 为本文气动结构耦合仿真方法（FSI）与传统CFD 方法得到的风洞沿程压力和气动阻力对比。从表2 可以看出：2 种方法的结果差异基本在5%以内[21]，说明了本文方法的可行性；而整体上CFD 仿真结果小于FSI 仿真结果，可能是CFD 方法未考虑气动与结构耦合效应的影响。

表2 FSI 与CFD 仿真结果[21]对比Table 2 Comparison results of FSI and CFD[21] simulation

3 管道消波措施分析

为确定多孔介质参数，通过数值仿真进行分析。仿真模型如图17 所示，给定二维模型，划分为高压区、低压区和多孔介质区，长度取1 m；高、低压区压差设置为5 000 Pa，通过高、低压区压差形成压缩波，压缩波强度为2 500 Pa，向右传播；在低压区中间平面处设置压力监测点，监测面平均压力随时间的变化规律。

图17 多孔介质参数数值计算模型Fig.17 Numerical calculation model of parameters for porous medium

根据一维等熵理论，压缩波在传播过程中遇到压力间断面会产生反射压缩波。本文计算了不同黏性系数（即不同流阻率）多孔介质条件下反射压缩波的强度。为了便于对比，首先计算了多孔介质黏性系数为0（即无多孔介质）时的压缩波反射工况，结果如图18 所示，压缩波与右端壁面接触后，向左方向反射时，无强度和能量损失。

图18 无多孔介质时压缩波的传播特性Fig.18 Compression wave propagation characteristics without porous medium

常规吸声棉的黏性系数为5.6 × 108m-2，给定多孔介质黏性系数，设置孔隙率为1，计算结果如图19所示：压缩波在与多孔介质接触后，形成强度约1 800 Pa 的反射压缩波，说明该多孔材料黏性系数过大，造成了压缩波在接触面的强反射；后续压力缓慢上升到4 890 Pa，说明多孔介质对压缩波有弛豫作用，且对压缩波有110 Pa 对应能量的耗散。

图19 多孔介质黏性系数为5.6 × 108 m-2 时压缩波的传播特性Fig.19 Compression wave propagation characteristics at viscous factor 5.6 × 108 m-2

降低多孔介质黏性系数至1 × 107m-2时，计算结果如图20 所示。从图中可以看出：压缩波与多孔介质接触后产生的反射波压力缓慢增大了300 Pa，压缩波与右端壁面接触后的反射压缩波通过多孔介质的强度约为1 400 Pa，且后续约有300 Pa 的压力缓升。此时壁面反射压缩波强度过高，且最后压力升高到5 000 Pa，说明此时多孔材料黏性较低，对能量耗散也较小。

图20 黏性系数1 × 107 m-2 时压缩波传播特性Fig.20 Compression wave propagation characteristics at viscous factor 1 × 107 m-2

图21 给出了无多孔介质、有高黏性多孔介质、有低黏性多孔介质时监测点的时域压力变化曲线。可以看出：多孔介质对压缩波具有显著的弛豫作用。由于压缩波强度与压力梯度成线性关系，所以多孔介质可以显著降低压缩波强度。多孔介质黏性越高，压缩波与其接触后的反射强度越高，能量耗散也越高。因此，要达到较好的降压效果，应选择黏性系数合理的多孔介质，使得压缩波尽可能地通过多孔介质，同时又有不错的弛豫作用。

图21 无多孔介质、有高/低黏性多孔介质时监测点的时域压力变化曲线Fig.21 Time-domain pressure change curve of monitoring point without porous medium and with high/low viscosity porous medium

经过计算分析，确定了一个较优的黏性系数：2 × 107m-2。如图22 所示，此时入射压缩波和反射压缩波强度分配合理，且都有不错的弛豫效果。

图22 多孔介质黏性系数为2 × 107 m-2 时的监测点压力变化曲线Fig.22 Variation curve of monitoring point ’s pressure at viscous factor 2 × 107 m-2

综上所述，尾端消波方案考虑填充黏性系数为2 × 107m-2、孔隙率≥0.95 的多孔介质。有关此类问题的研究较少，没有可以借鉴的方向，为达到最优设计目标，可通过仿真和实验方法开展关键技术研究解决以下问题：

1） 此种多孔介质非常规棉类多孔介质，须专门设计制造。

2） 考虑到压缩波的强度，对多孔介质的抗冲击性、抗疲劳性也有相应的要求。

3） 目前确定的黏性系数仅针对强度为 2 500 Pa的压缩波较优，对于其他强度的压缩波是否合适尚有待验证，还可以继续优化。

4） 尾端消波方案还可以尝试通过添加空气层、使用阻抗渐变的多孔介质来优化。

5） 多孔介质的固定结构须专门设计，应既能固定多孔介质，又不造成明显的压缩波反射。

4 结 论

本文基于CE/SE 方法对动模型管道内模型高速运动过程进行了气动结构耦合仿真，建立的每个工况均计算收敛，且可以获得模型运动过程中完整的流场结果。通过仿真结果评估了管道内压力波的传播规律和模型运动过程所受气动力的变化规律。开展了多孔介质参数的仿真设计，计算结果可为管道消波的设计方案比选提供参考。