边界数据浸入法在弱可压缩流动中的应用

赵体豪, 赵 欣

(1.北京理工大学 机械与车辆学院,北京 100081；2.北京理工大学 宇航学院,北京 100081)

计算成本及计算能力的限制使得我们在做计算时常常要将计算对象假定为不可压缩无粘流体. 这一物理假设虽然可以在一定程度上减小计算成本，但同时也会在结果中引入一定程度的误差. 为解决上述问题，本文提出了一种适用于求解弱可压缩粘性流动问题的数值解法. 该方法的数学推导基于边界数据浸入法(Boundary Data Immersion Method, BDIM)，并且该方法充分考虑了运动固体对流场计算结果的影响.

对于计算域内含有动边界的问题，浸入边界法(Immersed Boundary Method, IBM)是一种准确有效的求解方法. 该法由Peskin[1]于1972年首次提出，用于对人类心脏中血液流动的模拟. 其主要思想是将流场的复杂边界等效为一个体积力并施加于N-S动量方程中. 此后的几十年中，该法得到了不断发展且被广泛使用. 为改善体积不守恒问题，Peskin[2]引入新的投影算子替代原有方程中的哈密顿算子；自适应加密算法与IBM的结合[3-5]使得计算效率得到了提高;针对流体与运动刚体之间相互作用的研究，这些年来逐渐发展出了离散直接力法[6-8]、虚拟边界法[9-10]、反馈力法[11-13]等不同方法. 此外，如何提高边界浸入法的计算精度也引起了不少研究学者的关注. Yusof[14]将B样条(B-spline)与浸入边界法结合以获得高阶精度，且消除了时间步长的限制；Tseng[10]通过在边界上靠近固体一侧的区域进行虚拟网格重建的方法来提高计算精度；Yang[15]通过对离散delta函数的光滑处理，有效抑制了体积力中的非物理振荡；宫兆新等[16]通过虚拟解法验证了正则化δ函数的改变对于计算精度的提高没有任何帮助；李秋实等[17]通过改变力源的构造以及边界速度插值的方式获得二阶精度结果；及春宁等[18-20]将力源项的求解内嵌到压力泊松方程的迭代过程，得到了一种高精度的嵌入式迭代浸入边界法.

然而与原始浸入边界法相同，上述研究对象大多为不可压缩无粘流体，并未成功地将边界浸入法推广到可压缩流动中来. 对于不可压缩流动来说，使用边界浸入法只需要给出速度边界条件(即无滑移速度边界)即可；但是对于可压缩流动来说，因流体的一些变量，如密度、温度等都会在运动中发生较大的变化，所以需要除了速度边界条件以外的其他边界条件才可以顺利完成流场的求解. 为了将浸入边界法应用到可压缩流动中来，研究人员进行了大量工作[21-24]. Palma等[21-22]将IBM与求解预置可压缩N-S方程的解法相结合用于求解流场中复杂几何结构的问题；Qiu等[23]依次使用无滑移速度边界条件对流体的密度、速度、温度以及压力依次进行修正，推导出了一种边界条件浸入法(Boundary condition-enforced)；Schlanderer[25]提出一种适用于求解可压缩流动的浸入边界法——边界数据浸入法(Boundary Data Immersion Method). 该方法使用广义积分核公式将计算域内不同子域的控制方程进行重组，得到具有统一形式的方程，并且使得重组后的方程不仅在各子域内有效，且在边界处也可以实现光滑过渡. 此外，现在对于弱可压缩粘性流动问题的求解方法主要是predictor-corrector方法[26-27]或者人工可压缩的方法[28]. 在这里我们将借鉴Caltagirone[29]对流体弱可压缩粘性的处理方式，并结合Schlanderer[25]的边界数据浸入法推导出一种可以求解水下弱可压缩粘性流动问题的简单算法.

1 算法推导

对于单相弱可压缩粘性的流体，从连续性方程与流体状态方程出发，经过推导可以得到如下关于压力项的方程[29]：

(1)

式中:γ为比热容比系数，这里取6.1；B为水的刚度系数，其具有与压力相同的量纲， 取值345 MPa. 将上式代入式(1)，并对其进行离散化处理，有

(2)

式中，上标n表示物理量的当前时刻，n+1表示下一时刻.

对动量方程进行离散可以得到速度项的离散方程[25]如下：

(3)

式中:ρn为n时刻的流体密度；μ为流体的动力黏度系数，且上述物理量单位均采用国际单位制.

在得到流体的压力与速度的离散方程之后，可以对流场进行求解. 但在实际的流场计算中，常常会遇到整个计算域中存在多个不同物理性质的子域，因此需要得到求解整个流场的计算方法. 为解决上述问题，Schlanderer等[25]提出了基于广义积分核公式的边界数据浸入法. 这种方法通过积分核公式，将不同子域物理变量的控制方程和界面条件结合起来，实现了对整个计算域的模拟. 使用这种方法得到的具有统一形式的控制方程不仅保留了原系统中的物理变量特性，还做到了不同子域间的光滑过渡，使得我们可以对包含运动物体的可压缩粘性流动进行计算. 本文在这里直接使用可将各子域控制方程连结起来的控制方程[25]：

(4)

(5)

对于刚体内部，其压力方程具有如下形式：

将该式改写后可得

(6)

为满足刚体与流体边界处上压力连续变化的条件，刚体边界上压力方程修正后有如下形式：

(7)

其中式(7)中右边第二项恒为0. 流体的压力方程有以下形式:

(8)

将式(7)和式(8)进行离散化处理，并且结合式(4)整理可得适用于整个计算域中压强p′的方程为

(9)

至此，我们在弱可压缩粘性流动与BDIM的基础上，得到了可用于求解包含运动固体的可压缩粘性流动的完整算法.

2 算法验证

为对算法的有效性及准确性进行充分说明，我们在本章计算了3个二维圆柱的经典算例，并将计算结果与前人在相同工况下的计算结果进行了对比.

2.1 静止圆柱绕流

在本节，我们首先计算了静止圆柱绕流问题，验证了算法在计算包含静止固体的可压缩粘性流动问题时的有效性及准确性.

基于入流速度u与圆柱直径d，设置流动的雷诺数Re=ud/υ=100；设置计算域的大小为30d×8d，并在整个计算域上均布尺寸大小为1/30d的计算网格；时间步长为1×10-4. 计算采用入口速度条件，出口压力条件，其余边界采用无滑移边界条件.

表1Re=100时静止圆柱绕流的阻力、升力系数与斯特劳哈尔数

Tab.1 Drag and lift coefficients and Strouhal number for flow around a stationary cylinder atRe=100

方法C-DCL'St本文方法1.3600.3340.183文献[30]1.5010.3490.172 文献 [30]*1.4530.3390.169文献[31]1.3500.3390.165文献[32]1.3300.332—文献[33]1.3420.2970.166

注：*为增大计算域后所得结果.

图1中给出了当流动稳定时，圆柱表面的时均压力系数CP=(P-P)/u2的分布情况. 这里，提取数据的压力测点均为圆柱外侧最近的网格节点. 图中θ=0°与θ=180°分别对应圆柱上距离来流最近与最远的两点. 在图1中，我们将分别使用网格尺寸大小为1/20d、1/30d及1/40d得到的计算结果与实验结果[34-35]进行比对. 近些年技术水平的提高使得学者们将关注点放在高雷诺数的工况上，较少有人对低雷诺数的工况进行实验，故本文引用的参考文献时间较早. 该对比结果表明：使用本算法得到的计算结果与实验数据[34-35]整体趋势吻合较好. 将两个实验结果对比可以发现，较大的雷诺数会使得压力系数CP在50-180°小于雷诺数较小时的压力系数. 该现象可以解释为什么使用数值计算得到的结果在θ角度较大时会略小于实验数据. 此外，对比发现使用3种不同网格尺寸大小得到的计算结果相差不大，由此可得出计算结果与网格尺寸大小无关性的结论. 为使结果更加直观，在图2中提供了流动稳定后的涡量图. 至此，验证了算法在计算包含静止物体的可压缩粘性流动问题时的有效性及准确性.

2.2 振动圆柱绕流

上节使用算法计算了包含静止固体的可压缩粘性流动的算例，本节更进一步，验证当物体发生运动时算法的有效性及准确性. 本节考虑另一个二维算例——横向振动圆柱绕流的算例.

图1Re=100时，流动稳定后静止圆柱时均表面压力系数CP分布图

Fig.1 Distribution of time-averaged surface pressure coefficient on stationary cylinder when flow is stable atRe=100

图2 Re=100时静止圆柱扰流稳定时涡量图

Fig.2 Vortex diagram for flow past a stationary cylinder when the flow is stable atRe=100

基于入流速度u与圆柱直径d，设置流动的雷诺数Re=ud/υ=185. 设置计算域的大小为30d×10d，并在整个计算域上均布尺寸为1/30d的计算网格；时间步长为2×10-4. 计算采用入口速度条件，出口压力条件，其余边界采用无滑移边界条件. 圆柱在垂直于来流方向上以正弦规律进行振荡，其振荡幅度A=0.2d，振动频率为fo=0.156. 圆柱位移随时间的变化可用如下数学函数关系进行表达：

y(t)=Asin(2πfot).

表2对比了采用不同数值方法计算得到时均阻力系数以及阻力系数和升力系数的均方根. 通过对比，可以发现：使用我们的算法计算得到的3个数值在量级上与其他算法所得结果相同，具有较好的一致性；具体而言，计算得到的时均阻力系数值比其他算法得到的数值小1.1%～7.8%，阻力系数的均方根同文献[30,37]中结果几乎相同，而升力系数的均方根则略大于文献[30,36-37]中的计算结果. 产生这一差异化现象的原因是参考文献中采用数值计算方法均为不可压缩流体，其计算域的压力出口边界条件与本文采用的弱可压缩流体的压力边界条件不同.

此外，本文引入不同的圆柱振动特性进行计算，探究其变化对计算结果的影响. 文中采用改变单一变量的方法，计算了振动频率为0.156时，圆柱振动幅度A分别为0.2d、0.3d、0.4d、0.5d的4种情况以及振动幅度A为0.2d时，圆柱振动频率分别为0.156、0.3、0.45的3种情况. 对比结果见表3及表4. 从表3可以看出，振动幅度的增大，使得3种系数均呈现不断增大的趋势；而从表4可以看出，随着圆柱振动频率的增大，3个系数均有不同程度的增幅，其中，升力系数均方根的增幅远大于时均阻力系数与阻力系数均方根，达到了十倍以上.

表2Re=185时振动圆柱绕流的时均阻力系数、阻力及升力系数的均方根

Tab.2 Time-averaged drag coefficient and root mean square of drag and lift coefficient for flow around an oscillating cylinder atRe=185

方法C-D(CD)rms(CL)rms本文方法1.1820.0620.231文献[30,37]1.2820.0620.223文献[36]1.25—0.18

表3 圆柱振动频率为0.156时，不同振幅对应的时均阻力系数、阻力及升力系数的均方根

Tab.3 Time-averaged drag coefficient and root mean square of drag and lift coefficient associated with different oscillating amplitudes when the oscillating frequency is 0.156

圆柱振幅/dC-D(CD)rms(CL)rms0.21.1820.0620.2310.31.2650.1230.2750.41.3470.1870.2880.51.4260.2570.306

表4 圆柱振动幅度为0.2d时，不同振动频率对应的时均阻力系数、阻力及升力系数的均方根

Tab.4 Time-averageddrag coefficient and root mean square of drag and lift coefficient associated with different oscillating frequencies when the oscillating amplitude is 0.2d

圆柱振频C-D(CD)rms(CL)rms0.1561.1820.0620.2310.301.3790.1451.3090.451.5060.2432.494

图3和图4给出了阻力系数及升力系数随圆柱位移变化的图像. 从图中可以看出，无论是阻力还是升力，其数值大小均为对称分布，其中，前者是以静止位置(yc=0)轴对称，而后者则以(0,0)点呈中心对称. 此外，在图5中也给出了流动稳定时的涡量图.

图3 Re=185时振动圆柱扰流的阻力系数随位移变化图

Fig.3 Variation of drag coefficient with displacement for flow past an oscillating cylinder atRe=185

图4 Re=185时振动圆柱扰流的升力系数随位移变化图

Fig.4 Variation of lift coefficient with displacement for flow past an oscillating cylinder atRe=185

2.3 横向振动圆柱

上述2个算例对算法的有效性及准确性做了验证，本节通过验证横向振动圆柱，对算法在声学层面的应用进行验证.

计算域的大小为30d×30d，整个计算域的网格尺寸为0.05d，时间步长为0.01. 圆柱在横向上的振动规律同2.2节，即圆柱位移与时间的函数关系式为

y(t)=Asin(2πfot).

图5 Re=185时振荡圆柱流动稳定时涡量图

Fig.5 Vortex diagram for flow past an oscillating cylinder when flow is stable atRe=185

图6给出了圆柱在振动稳定时，计算域内压力云图变化的整个周期. 通过图6可以看出，圆柱在发生振动时所产生的压力脉动在声学上等同于一个偶极子，这与理论分析相一致. 此外，从图7中可以看出，对方向谱而言，计算结果同偶极子振动的理论解有很好的一致性. 从而验证了算法在计算包含运动物体的可压缩粘性流动时声学层面的有效性.

(a)振动周期阶段一 (b)振动周期阶段二图6 圆柱振动压力云图Fig.6 Pressure contour of oscillating cylinder

3 结 论

相比于传统的浸入边界法，本文所采用的边界数据浸入法不仅可对包含运动物体的流域进行求解，更可与弱可压缩粘性流结合起来以考虑流体的可压缩性及粘性对计算结果的影响. 本文从可压缩粘性流的连续性方程及水的状态方程出发，使用边界数据浸入法将运动物体的固体子域与流体子域的变量方程进行重组，严格推导了适用于整个流场的压力与速度方程. 该方程不仅有效保留了在各子域内的所有物理性质，而且成功地在各子域之间交界面上实现了光滑过渡.

此外，通过3个二维经典算例，在力学与声学层面上验证了算法在计算包含静止物体及运动物体的流场时的有效性. 通过将计算所得结果与其他算法所得结果进行对比，对算法的准确性进行了验证说明. 相较于目前已有的可压缩粘性流动问题的求解方法，本算法推导简单，无需过多假设，且计算结果与实验数据具有良好的吻合程度.