计算成本及计算能力的限制使得我们在做计算时常常要将计算对象假定为不可压缩无粘流体.这一物理假设虽然可以在一定程度上减小计算成本，但同时也会在结果中引入一定程度的误差.为解决上述问题，本文提出了一种适用于求解弱可压缩粘性流动问题的数值解法.该方法的数学推导基于边界数据浸入法(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]：

$ \frac{{\partial p}}{{\partial t}} = - (\overrightarrow u \cdot \nabla p + \rho c_T^2\nabla \cdot \overrightarrow u ). $

(1)

式中，p为流体压强，$ \overrightarrow u $是流体速度，c_T是在温度T下的声速，且上述所有物理量单位均采用国际单位制.当以弱可压缩的水作为研究对象的时候，存在如下关系：

$ \rho c_T^2 = (p + \gamma B). $

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

$ {p^{n + 1}} = {p^n} - \Delta t \cdot [\nabla {p^n} + ({p^n} + \gamma B\nabla \cdot {\overrightarrow u ^n})], $

(2)

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

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

$ {\overrightarrow u ^{n + 1}} = {\overrightarrow u ^n} - \Delta t \cdot {\overrightarrow u ^n}\nabla \cdot {\overrightarrow u ^n} + \frac{{\Delta t}}{{{\rho ^n}}}[\mu {\nabla ^2}{\overrightarrow u ^n} - \nabla {p^{n + 1}}]. $

(3)

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

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

$ {\varphi _\varepsilon } = {b_f}(\varphi , \overrightarrow x , t)\mu _0^{_\varepsilon } + {b_s}(\varphi , \overrightarrow x , t)(1 - \mu _0^{_\varepsilon }). $

(4)

式中: ${b_f}(\varphi , \overrightarrow x , t) $及${b_s}(\varphi , \overrightarrow x , t) $分别为流体子域及固体子域的任意一个物理变量φ的控制方程；μ₀^ε为一无量纲的光滑过渡函数，其取值范围为0≤μ₀^ε≤1.

相较于流体速度，固体速度$\overrightarrow V (\overrightarrow x , t) $可以直接求得.将式(3)中求得的流体速度${\overrightarrow u ^{n + 1}}(\overrightarrow x , t) $与固体速度$ \overrightarrow V (\overrightarrow x , t)$代入式(4)后，可以得到适用于求解整个计算域中速度$\overrightarrow U (\overrightarrow x , t) $的方程如下：

$ \overrightarrow U (\overrightarrow x , t) = {\overrightarrow u ^{n + 1}}(\overrightarrow x , t)\mu _0^{_\varepsilon } + {\overrightarrow V ^{n + 1}}(\overrightarrow x , t)(1 - \mu _0^{_\varepsilon }). $

(5)

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

$ \frac{{Dp}}{{Dt}} = 0, $

将该式改写后可得

$ \frac{{\partial p}}{{\partial t}} = - \overrightarrow V \frac{{\partial p}}{{\partial x}}. $

(6)

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

$ \frac{{\partial p}}{{\partial t}} = - \overrightarrow V \frac{{\partial p}}{{\partial x}} + (p + \gamma B)\nabla \cdot \overrightarrow V . $

(7)

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

$ \frac{{\partial p}}{{\partial t}} = - \overrightarrow u \frac{{\partial p}}{{\partial x}} + (p + \gamma B)\nabla \cdot \overrightarrow u . $

(8)

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

$ p' = {p^n} - \Delta t \cdot [{\overrightarrow U ^m} \cdot \nabla {p^n} + ({p^n} + \gamma B)\nabla \cdot {\overrightarrow u ^m}]. $

(9)

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

2 算法验证

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

2.1 静止圆柱绕流

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

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

数值计算结果整理在表 1中，并且给出了使用其他数值算法所得结果以供对比.通过表 1中的对比，可以看出：总体而言，使用新算法得到的时均阻力系数C_D、升力系数振荡的幅值C′_L以及斯特劳哈尔数($ St = fd/{u_\infty }$)在量级上与使用其他算法所得到的结果相同.具体而言，使用我们的算法求得的时均阻力系数与其他算法所得结果的吻合程度较好；计算结果得到的升力系数振荡的幅值相比于其他算法得到的数值较小，但相差值最大不超过4.5%；相比于前两者，计算得到的斯特劳哈尔数与使用其他算法得到的数值相差较大，比其他值高约6.4%~10.91%左右，而这意味着我们算得的升力系数的振荡周期略小于其他数值算法算得的升力系数振荡周期.斯特劳哈尔数表征了涡脱落的频率，我们计算得到的结果较大，表明涡脱落的速度较快.产生这一差异化现象的原因是参考文献中采用数值计算方法均为不可压缩流体，其计算域的压力出口边界条件与本文采用的弱可压缩流体的压力边界条件不同.

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

2.2 振动圆柱绕流

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

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

$ y(t) = Asin(2\pi {f_0}t). $

表 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个系数均有不同程度的增幅，其中，升力系数均方根的增幅远大于时均阻力系数与阻力系数均方根，达到了十倍以上.

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

2.3 横向振动圆柱

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

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

$ y(t) = A{\rm sin}(2\pi {f_0}t). $

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

3 结论

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

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