一种求解低速跨流域流固耦合问题的离散统一气体动理学格式

王 勇，刘 沙，2，卓丛山，2，钟诚文，2，*

(1. 西北工业大学 航空学院，西安 710072；2. 西北工业大学 翼型、叶栅空气动力学重点实验室，西安 710072)

0 引 言

运动边界和流固耦合（Fluid-structure interaction，FSI）问题是航空航天工程、船舶和海洋工程、土木和交通工程等领域经常面临的多场耦合问题[1-2]。数值模拟FSI 问题目前已发展大量方法。近些年，随着技术进步，以微纳米机电系统（micro-electro-mechanicalsystems，MEMS）和超低轨道航天器为代表的应用领域面临跨流域运动边界和流固耦合问题，如MEMS中微尺度梁的振动问题[3]、航天器返回过程[4]和超低轨道航天器与外层稀薄大气相互作用问题[5]。在这类问题中，由于特征尺度的显著减小或气体密度显著降低，流体的稀薄气体效应是这类流固耦合问题需要考虑的重要影响因素之一。根据克努森数（Kn）定义[6]，流动可划分为连续流、滑移流、过渡流和自由分子流四个流域。由于上述新问题中流体流域通常处于滑移流域或过渡流域，而传统的流固耦合求解方法中流体部分通常采用基于Navier-Stokes 方程的适用于连续流域的宏观方法，因此有必要发展可考虑稀薄气体效应影响的流固耦合求解新方法。

近年来，在跨流域流动数值模拟方法方面，离散统一气体动理学格式（discrete unified gas kinetic scheme，DUGKS）[7]这一新兴数值方法得到了广泛关注。该方法在格式通量计算部分，借鉴了格子Boltzmann 方法（lattice Boltzmann method，LBM）[8]的特征线构造思想，较统一气体动理学格式（unified gas kinetic scheme，UGKS）[9]显著简化，计算效率进一步提升。在网格实施方面，由于标准LBM 采用迁移-碰撞模式，Boltzmann 方程的对流项隐藏在分布函数沿特征线的迁移过程中[10]，因此计算网格需采用均匀笛卡尔网格；而DUGKS 在实施过程中由于采用对流项显式求解思想，尽管计算复杂度有所提升，数值耗散也有所增加，但可采用的网格也更加灵活，结构、非结构及混合网格均可，对不规则计算域的适用性显著提高。在数值耗散方面，与传统的有限体积LBM[11-12]相比，由于对流项计算过程中引入迁移-碰撞思想，因此数值耗散小，对网格尺度和时间步长约束小。目前该类方法已具有全流域流动计算能力，并在快速发展中[13-14]。考虑到运动边界问题的重要性，目前已有开展少量工作[15]，仍需进一步完善。

在运动边界数值模拟方面，网格技术通常分为基于动网格技术的任意拉格朗日-欧拉法（arbitrary Lagrangian-Eulerian，ALE）[16]和基于笛卡尔网格的浸入边界法（immersed boundary method，IBM）[17]。ALE可采用贴体网格，在模拟高雷诺数流动及结构运动模式相对简单的FSI 问题时有较大优势，计算效率较高，如航空工程中的颤振数值模拟。IBM 由于壁面附近采用笛卡尔网格，适用于雷诺数较低的流动问题，模拟大变形动边界问题有显著优势。目前，基于DUGKS的ALE[15]和IBM[18]均有发展，其中ALE-DUGKS 在模拟考虑稀薄气体效应影响的运动边界问题时有良好表现。而在低速连续流域范围，以LBM、DUGKS为代表的介观数值方法，由于无需求解宏观方法中的压力泊松方程[19]，方法实施也相对简单。ALEDUGKS 在低速连续流和稀薄流的独特优势是为进一步发展低速跨流域运动边界和流固耦合数值方法奠定基础。在流固耦合数值模拟方面，发展较为成熟的是松耦合计算方法[20]，即在每个数值模拟时间步长内采用计算流体力学（computational fluid dynamics，CFD）和计算结构动力学（computational structural dynamics，CSD）各自成熟的求解方法，并在流-固耦合界面进行数据交换实现流体域和固体域的耦合推进。

综上，本文的主要工作是针对MEMS 和航天工程面临的跨流域流动运动边界和流固耦合问题，采用新兴跨流域数值方法DUGKS，将已发展的运动边界ALE-DUGKS 扩展至低速，同时兼顾高速流动的流固耦合计算框架，利用网格的贴体性采用较少的网格量更高效地模拟边界层流动。针对连续流域圆柱绕流强迫振动和自由振动两个典型算例进行模拟，验证方法的计算能力，并进一步采用跨流域求解方法验证方法的跨流域计算能力和拓展性。

1 基于DUGKS 的流固耦合求解方法

1.1 ALE-DUGKS 格式

ALE-DUGKS 是王勇等[15]在原始DUGKS[7]基础上构造的跨流域流动运动边界求解格式。原始DUGKS 格式构造基于Boltzmann-BGK 方程：

式中，下标k为控制体界面编号，下标b为界面中点处相应变量值。式（8）和式（9）中的V和S分别为控

1.2 流固耦合计算方法

ALE-DUGKS 具备运动边界问题求解能力，耦合计算结构动力学并采用松耦合求解方法后，可将DUGKS发展成可考虑稀薄气体效应的跨流域流固耦合计算框架。在实施过程中，流体域网格变形采用Laplace光顺法[23]，结构域控制方程为：

2 算例验证

2.1 静止圆柱绕流算例

为了验证DUGKS 方法并为后续算例提供对比数据，首先进行静止圆柱绕流模拟。离散速度模型采用D2Q9 模型。参考文献[25，26]的算例设置，计算域见图1，此处D为圆柱直径，根据阻塞率5%布置上、下边界远场距离。计算域左侧边界和上、下边界设置为入口边界条件，边界外法线为nb；对于进入计算域的分布函数方向eα·nb＜0，分布函数给定为平衡态；离开计算域的分布函数方向由内场分布函数插值获得。计算域右侧边界设置为出口边界条件，分布函数由内场插值获得。图2 为计算网格，总网格单元为37 506，其中壁面附近采用四边形贴体网格，最小网格尺度为D/256，其余为三角形网格。初始化来流速度为u=U∞=0.05，v=0； 初始密度 ρ∞=1.0。雷诺数范围为Re= 60～150。升力系数CL定义为：

图1 圆柱绕流计算域和边界条件Fig. 1 The configuration and boundary conditions for the flow around a circular cylinder

图2 圆柱绕流计算网格Fig. 2 A mesh for the flow around a circular cylinder

式中f为自然涡脱落频率。本文与Prasanth[25]和He[26]的数值模拟结果以及Williamson[27]的平行涡脱落模型吻合良好。在雷诺数Re= 60～80 范围内数值结果略高于平行涡脱落模型，主要原因是随着雷诺数降低，黏性影响范围增大，远场边界条件对结果影响显著。减小阻塞率后，可明显改善数值结果和模型之间的差异。

图3 三组雷诺数下圆柱绕流升力系数时间发展历程Fig. 3 Time evolutions of lift coefficients for flows with three Reynolds numbers around a circular cylinder

图4 Re = 60～150 静止圆柱绕流St 对比Fig. 4 A comparison of St for flows around a stationary cylinder at Re = 60～150

2.2 Re = 100 圆柱横向强迫振动算例

在静止圆柱绕流算例基础上，进一步模拟圆柱横向强迫振动算例以验证算法的运动边界模拟能力。计算参数与算例1 相同，雷诺数Re= 100。圆柱y方向运动方程为：

式中，A为圆柱振幅，fy为圆柱结构频率，fRe=100为静止圆柱绕流自然涡脱落频率，A*为圆柱无量纲振幅，f*为与自然涡脱落频率的比值。在文献[28]中，Kumar 等将圆柱涡脱落分为非锁频区、过渡区和锁频区。图5 为各区分界线图，与Kumar 等结果整体十分吻合。图6 为五组A*和f*组合下的升阻力系数时间发展历程和升力系数功率谱图，其中阻力系数CD为将式（20）中FL替 换成FD。当无量纲频率f*＜1 时，非锁频区（图6（a））和锁频区（图6（b））有清晰的分界线。在非锁频区，由于流体自然涡脱落频率和结构频率接近，升阻力系数出现类似“拍”现象，功率谱图中流动自然涡脱落频率占优且包含其他无规律频率信息；在锁频区，涡脱落频率锁定在结构频率，同时其他频率与主频呈现清晰的倍频现象。当无量纲频率f*＞ 1 时，非锁频区和锁频区之间存在过渡区。在非锁频区和过渡区，升阻力系数同样出现“拍”现象（图6（c）和图6（d））。当流动自然涡脱落频率占优时，流动呈现非锁频现象；当结构自然频率占优时，流动则处在过渡区范围。对于功率谱图，非锁频区和过渡区除主频外，还包含其他无规律频率信息，锁频区其他频率信息则呈现清晰的倍频现象（图6（e））。

图5 Re = 100 圆柱横向强迫振动非锁频区、过渡区和锁频区Fig. 5 No lock-in, transitional, and lock-in regimes for flows around a forced-oscillating circular cylinder at Re = 100

图6 五组不同振幅A＊和频率f ＊下升阻力系数时间发展历程和升力系数功率谱对比Fig. 6 Time evolutions of lift and drag coefficients, and the power spectra of lift coefficient at five different amplitudes and frequencies of oscillation

2.3 Re = 100 圆柱自由振动算例

为验证ALE-DUGKS 模拟流固耦合问题的能力，首先对Re= 100 圆柱自由振动进行模拟。结构自由振动控制方程为x和y方向两自由度无阻尼振动方程，对应于式（18）中的质量矩阵M、刚度矩阵K和外力矩阵Fext分别为：时间内，圆柱处在锁定状态，流场发展成周期性流动后释放圆柱。图7 为U*= 6.25 时两个方向位移时间发展历程和收敛后的耦合位移，其中x和y为圆柱在两个方向的位移值。可以发现周期性非定常升力使得圆柱在y方向产生较大位移，非定常阻力则使圆柱在x方向运动到平衡位置后以较小振幅振动，且两方向耦合位移呈现类似“8”的结构；圆柱在其他折减速度时有类似现象。图8 为圆柱y方向最大位移和涡致振动频率对比，其中ymax/D为无量纲最大位移，St= 0.166 为静止圆柱绕流的无量纲涡脱落频率。本文与Singh[29]、Zhang[30]等数值结果吻合良好，主要差异在锁频和非锁频的过渡区间。图9 为U*= 4.0 时y方向位移时间历程，与图7 结果不同，此时流动处在非锁频状态，圆柱位移很小。图10 则为流动处在锁频和非锁频的过渡状态，需要较长的数值模拟时间获得周期性振动结果。此外，数值模拟发现，阻塞率较小时，圆柱最大位移偏高，即远场距离对圆柱位移有压缩效应。这也是本文与文献[29]相同，取5%阻塞率的原因。

图7 圆柱自由振动折减速度U＊ = 6.25 结构位移Fig. 7 Displacements of the free-oscillating circular cylinder at U＊ = 6.25

图8 Re = 100 圆柱自由振动y 方向最大位移和涡脱落频率对比Fig. 8 Comparisons of the maximum y-direction displacement and vortex shedding frequency for the flow around a freeoscillating circular cylinder at Re = 100

图9 圆柱自由振动U＊ = 4.0 时y 方向位移时间发展历程Fig. 9 Time evolution of the y-direction displacement of a freeoscillating circular cylinder at U＊ = 4.0

图10 圆柱自由振动U＊ = 7.8 时y 方向位移时间发展历程Fig. 10 Time evolution of the y-direction displacement of a freeoscillating circular cylinder at U＊ = 7.8

2.4 Re = 60～150 圆柱自由振动算例

与Re= 100 圆柱自由振动算例不同，该算例主要验证结构自然频率固定时来流风速变化对结构位移和频率的影响。结构参数仍为式（23）～式（25），此时无量纲频率给定为Fs= 16.6/Re，其中16.6 使得Re=100 时结构无量纲频率约等于静止圆柱绕流无量纲涡脱落频率。图11 为y方向最大位移和涡脱落频率对比，与Prasanth 等[25]的数值结果十分吻合，仅在锁频和非锁频的过渡区间有差异。主要原因是Prasanth 等在计算过程中通过动态调节雷诺数（增加或减小）以分析风速的动态变化对过渡区流动机理的影响；本文主要为方法验证，忽略动态影响机理分析，计算中每给定一个雷诺数获得收敛结果后完成计算。对于y方向位移和升阻力系数，在非锁频和锁频区与算例2.3 的Re= 100 自由振动算例类似，但在过渡区较为复杂。图12 为Re= 130 的时间历程，经过频谱分析，当无量纲时间小于1 800 时，y方向位移不断增加，此时涡脱落频率中圆柱的自然涡脱落频率占优；当时间大于1 800 后，涡脱落频率锁定在结构自然频率，圆柱位移也收敛到稳定的周期性振动状态。

图11 Re = 60～150 圆柱自由振动y 方向最大位移和涡脱落频率对比Fig. 11 Comparisons of the maximum y-direction displacement and vortex shedding frequency for flows around a free-oscillating circular cylinder at Re = 60～150

图12 Re = 130 圆柱自由振动y 方向位移和升阻力系数时间发展历程Fig. 12 Time evolutions of the y-direction displacement, and lift and drag coefficients for the flow around a free-oscillating circular cylinder at Re = 130

2.5 Re = 60 圆柱横向强迫振动连续流和稀薄流对比

为验证算法的扩展能力，采用跨流域计算框架对Re= 60 圆柱横向强迫振动进行数值模拟。跨流域计算中，克努森数Kn、马赫数Ma和雷诺数Re关系为[31]：

图13 Ma = 0.086 6、Re = 60 圆柱强迫振动升力系数时间发展历程Fig. 13 Time evolutions of lift coefficients for the flow around a forced-oscillating circular cylinder at Re = 60 and Ma = 0.086 6

针对弱可压缩流动，式（7）和式（27）求解见文献[15]。离散速度采用28×28 高斯点。Ma= 0.4 时圆柱自然涡脱落频率St= 0.132，与Canuto 等[32]结果吻合，略低于Ma= 0.086 6 时的结果。图14 为两组马赫数下升力系数时间发展历程对比，随着马赫数增加，稀薄程度增加，升力系数呈减小趋势。Canuto 等[32]采用宏观方法和无滑移边界分析了压缩性对静止圆柱绕流的流动特性影响，因方法的局限性无法考虑稀薄气体效应影响，数值结果在低雷诺数时存在一定偏差。但与本文静止圆柱和振动圆柱趋势一致，即Re=60 时马赫数对St影响较小，但最大升力系数降低。

图14 两组马赫数下圆柱强迫振动升力系数时间发展历程Fig. 14 Time evolutions of lift coefficients for flows around a forced-oscillating circular cylinder at two different Ma

2.6 Ma = 5.0 圆柱横向强迫振动稀薄流算例

2.2～2.5 节的算例验证了算法的低速跨流域运动

图15 Ma = 5.0 圆柱强迫振动物理网格和速度网格Fig. 15 The meshes in the (a) physical- and (b) velocityspace[21]for the flow around a forced-oscillating circular cylinder at Ma = 5.0

其中，A*为圆柱无量纲振幅，U*= 100 为折减速度。图17 为A*= 0.05 时的升力系数时间发展历程。由于振幅较小，振动频率较低，圆柱所受升力系数较小。此外，与静止圆柱绕流相同，马赫数固定时，来流越稀薄，升力越大。图18 为Kn= 0.1、A*分 别取0.05 和0.1 时的升力系数时间发展历程对比。由于来流为稀薄流，圆柱表面无分离涡和涡脱落现象，流动呈现为附着流的未失稳状态，此时升力系数与振幅表现为线性增长关系。考虑到不同马赫数、克努森数、振幅A*和折减速度U*组合下的圆柱强迫振动流动分析超出了本文的讨论范围，课题组将会在接下来的论文中进行系统讨论。

图17 两组努森数下圆柱强迫振动升力系数时间发展历程Fig. 17 Time evolutions of lift coefficients for flows around a forced-oscillating circular cylinder at two different Kn numbers

图 18 Kn = 0.1 时两组无量纲振幅A＊的升力系数时间发展历程Fig. 18 Time evolutions of lift coefficients for flows around a forced-oscillating circular cylinder at two different A＊ values

3 结 论

本文在原始DUGKS 基础上，耦合任意拉格朗日-欧拉方法和计算结构动力学，实现了可考虑稀薄气体效应的低速兼顾高速流动的流固耦合问题数值模拟框架。针对低速连续流域的圆柱强迫振动和自由振动两个典型算例，采用LBM 中高效的D2Q9 离散速度模型，该框架可获得与宏观方法一致的数值结果。鉴于DUGKS 的全流域、全速域计算能力，更改离散速度模型后，即可具备模拟跨流域运动边界模拟能力。在算法实施方面，本文初步验证了显式ALEDUGKS 模拟流固耦合问题能力，未来可发展隐式格式，进一步提高其计算效率。在应用方面，考虑到稀薄气体效应对新型跨流域飞行器低雷诺数可压缩绕流场的显著影响[35]，以及国内外较少开展低雷诺数可压缩圆柱绕流分析，本文后续将深入研究圆柱静、动态气动特性数值模拟。