钝体绕流，特别是柱体绕流问题在工程实际中经常遇到，如风吹过高层建筑物、海水流过石油钻井平台、海流流经海洋立管等。其共同的特点是，在一定雷诺数范围内，粘性流体流过柱体结构物时，会发生边界层的分离，进而对结构物产生周期性作用。因此，准确计算出流体对柱体结构的作用荷载及其在荷载作用下的动力响应，具有重要的科学价值和工程意义。

许多学者采用试验方法和数值模拟方法开展了钝体绕流问题的研究。Tritton^[1]在实验室内通过观察石英试管的弯曲度得到Re在0.5~100范围内圆柱绕流的阻力，与理论值吻合良好。Kawaguti等^[2]对Re在1~100圆柱绕流问题进行了数值模拟，并考虑了压力和摩擦力对阻力的贡献，结果表明压力阻力项在阻力中占有主导地位。Rajani等^[3]利用二阶精度的中心差分格式对对流方程进行离散，对圆柱绕流进行了二维和三维数值模拟，将两者的涡脱频率和升力系数进行对比，发现大于某一临界雷诺数时，二维结果偏大，三维效应不可忽略。方柱绕流也是钝体绕流的一种典型形式。Okajima^[4]研究了不同类型的矩形柱体的涡旋脱落频率，并指出Strouhal数是雷诺数的函数。Sohankar等^[5]引入二阶精度的Van Leer迎风格式对对流项进行求解，指出Re>50时会有稳定的卡门涡街产生。Breuer等^[6]将格子玻尔兹曼方法(lattice Boltzmann method，LBM)和有限体积法 (finite volume method，FVM)分别应用到二维方柱绕流数值模拟中，两种方法的计算结果显示了较好的一致性。Sen等^[7]利用有限元方法(finite element method，FEM)对Re<150时的方柱绕流进行二维数值模拟，发现最初的分离点是在基点处而不是在后角点处。潘小强等^[8]利用Fluent分别模拟了二维圆柱和方柱绕流的流场情况，结果表明同样雷诺数的圆柱和方柱绕流，圆柱绕流更易得到稳态非定常解。周云龙等^[9]采用FVM模拟了二维圆柱和方柱的流动分离状态，发现两者最大的不同点在于分离点的不同，圆柱绕流没有固定分离点而方柱绕流固定分离点在其前后角点处。不同柱体绕流过程中的水动力特性的差异，仍未明朗，这将是本文工作的重点。

本文采用自主研发的CIP方法模型^{[10, 11]}，通过浸入边界方法处理流-固耦合问题，对不同雷诺数(Re<300)的圆柱绕流和方柱绕流问题进行数值模拟，并通过与文献结果的比较验证模型的有效性。为了得到稳定的流场信息，本文采用了与传统差分格式不同的高阶格式CIP方法对N-S方程进行了离散，该方法通过对网格内结点函数值及其导数值联立，实现了三阶精度的差分求解对流方程，具有格式稳定和低耗散性的特点；同时，在直角网格系统内，通过引入浸入边界方法，有效处理了流-固耦合问题，可有效避免动态网格系统中的大量信息交换，保证了模型的计算效率，为进一步模拟运动物体做好准备。本文将从多个角度分析圆柱绕流与方柱绕流的区别。

1 数学模型的建立 1.1 控制方程

本流场模型以二维N-S方程为控制方程，其矢量形式为

$ \triangledown \cdot u=0 $

(1)

$ \frac{\partial u}{\partial t}+\left( u\cdot \triangledown \right)u=-\frac{1}{\rho }\triangledown p+\frac{\mu }{\rho }{{\triangledown }^{2}}u+F $

(2)

为了把CIP方法应用于对流方程，对方程(2)取空间导数，可得到

$ \frac{\partial \left( {{\partial }_{i}}u \right)}{\partial t}+\left( u\cdot \triangledown \right)\left( {{\partial }_{i}}u \right)=-{{\partial }_{i}}u\cdot \triangledown u-{{\partial }_{i}}\left( \frac{1}{\rho }\triangledown p+\frac{\mu }{\rho }{{\triangledown }^{2}}u-F \right) $

(3)

模型在直角笛卡尔系统下建立，采用多相流理论处理固-液的相互作用，满足下面的控制方程：

$ \frac{\partial {{\varphi }_{m}}}{\partial t}+u\cdot \triangledown {{\varphi }_{m}}=0 $

(4)

式中：m=1,2,φ₁为液体相，φ₂为固体相，在一个网格内满足φ₁+φ₂=1。固体相的处理基于浸入边界方法得到，将在下面的部分给出解释。网格内的流体特性可用下式来表示：

$ \lambda =\sum\limits_{m=1}^{2}{{{\varphi }_{m}}{{\lambda }_{m}}} $

(5)

式中：λ为密度ρ或者粘性系数μ。

1.2 数值方法 1.2.1 CIP方法

CIP 方法是Takewaki等^[12]于20世纪80年代中期提出并发展起来，用于求解双曲型偏微分方程的一种有效的数值计算方法。其基本原理是基于空间网格点的变量值及其空间导数值，利用三次多项式进行插值近似，反演出网格单元内部变量的真实信息，得到三阶精度的显式格式。为了更好地解释CIP方法，考虑简单的常系数一维对流方程：

$ \frac{{\partial f}}{{\partial t}} + u\frac{{\partial f}}{{\partial x}} = 0 $

(6)

式中：u为大于零的常数。此类偏微分方程可用不同差分方法求解。下面将给出CIP方法的工作原理。图 1(a)~(c)给出了一阶迎风差分方法的工作原理，对于一阶差分格式来说，它利用直线的方式联立相邻两个节点的信息来工作，而忽略了网格内部的信息，导致较大的数值耗散。为了真实再现网格内部的信息，有必要求助高阶差分方法，而对于常规高阶差分的建立需要更多的网格点的信息。CIP采用一种独特的方式，在一个网格内实现了高阶差分格式，通过利用空间网格点的变量值及其空间导数值，来描述该网格内的信息，可真实再现网格内的信息。

通过对方程(6)求关于x的偏导，得到如下的空间导数方程：

$ \frac{{\partial g}}{{\partial t}} + u\frac{{\partial g}}{{\partial x}} = - \frac{{\partial u}}{{\partial x}} $

(7)

式中：$g = \frac{{\partial f}}{{\partial x}}$。因对流速度u为常数，方程(7)右边项为0。这样，方程(7)与方程(6)有相同的形式。对于u>0的情况，在迎风向网格单元[x_i-1,x_i]内n时刻的剖面函数fⁿ可以近似为

$ {F_i}^n = {a_i}{\left( {x - {x_i}} \right)^3} + {b_i}{\left( {x - {x_i}} \right)^2} + {c_i}\left( {x - {x_i}} \right) + {d_i} $

(8)

在n+1时刻的单元格剖面函数fⁿ⁺¹可以通过将n时刻的剖面函数fⁿ平移-u△t得到，函数f和g的时间演变可以通过下面的拉格朗日变换得到

$ {{f}_{i}}^{n+1}={{F}_{i}}^{n}\left( {{x}_{i}}-u\vartriangle t \right),{{g}_{i}}^{n+1}=\text{d}{{F}_{i}}^{n}\left( {{x}_{i}}-u\vartriangle t \right)/dx $

(9)

因此，从CIP对流格式采用了拉格朗日常量解法(Lagrangian invariant solution)的角度来看，CIP对流格式又称为半拉格朗日方法，如图 2所示。式(8)中的4个未知系数可由已知量f_iⁿ、f_i-1ⁿ、g_iⁿ和g_i-1ⁿ通过下列关系式来确定：

$ \left\{ \begin{matrix} {{a}_{i}}=\frac{{{g}_{i}}^{n}+{{g}_{i}}^{n}}{\vartriangle {{x}^{2}}}-\frac{2\left( {{f}_{i}}^{n}-{{f}_{i-1}}^{n} \right)}{\vartriangle {{x}^{3}}},{{c}_{i}}={{g}_{i}}^{n} \\ {{b}_{i}}=\frac{2{{g}_{i}}^{n}+{{g}_{i}}^{n}}{\vartriangle x}-\frac{3\left( {{f}_{i}}^{n}-{{f}_{i-1}}^{n} \right)}{\vartriangle {{x}^{2}}},{{d}_{i}}={{f}_{i}}^{n} \\ \end{matrix} \right. $

(10)

CIP方法采用一种独特的方式，在一个网格内实现了高阶差分格式，使得本文的数值模型可以应用于复杂流动问题的模拟，这也是本文所采用的差分方法与其他方法的不同之处。

1.2.2 时间积分

采用分步算法对动量方程进行时间积分。首先忽略扩散项和压力项，只考虑对流项的中间速度；其次求解扩散项；然后求解压力方程，计算下一时间步的压力；最后考虑压力梯度项，计算速度的最后值。设△t为时间步长，在t=n△t到t=(n+1)△t时刻的计算时间内，具体的时间积分过程如下：

1) 对流项(Ⅰ)：

$ \left\{ \begin{matrix} \frac{\partial u}{\partial t}+\left( u\cdot \triangledown \right)u=0 \\ \frac{\partial \left( {{\partial }_{i}}u \right)}{\partial t}+\left( u\cdot \triangledown \right)\left( {{\partial }_{i}}u \right)=0 \\ \end{matrix} \right. $

(11)

通过CIP方法^[13]可得到方程的解：

$ \left\{ \begin{matrix} {{u}^{*}}=\text{X}\left( x-u\vartriangle t \right) \\ {{\left( {{\partial }_{i}}u \right)}^{*}}\left( x \right)=\frac{\partial {{\text{X}}^{*}}}{\partial {{x}_{i}}}\left( x-u\vartriangle t \right) \\ \end{matrix} \right. $

(12)

式中：‘*’为对流项计算结束后的中间时间标志。

2)非对流项(Ⅰ):

扩散项的计算：

$ \left\{ \begin{matrix} \frac{\partial u}{\partial t}=\frac{\mu }{\rho }{{\triangledown }^{2}}u+F \\ \frac{\partial \left( {{\partial }_{i}}u \right)}{\partial t}=-{{\partial }_{i}}u\cdot \triangledown u{{\partial }_{i}}\left( \frac{\mu }{\rho }{{\triangledown }^{2}}u+F \right) \\ \end{matrix} \right. $

(13)

扩散项的时间离散采用显示格式，中间速度可表示为下面的形式：

$ \left\{ \begin{matrix} \frac{{{u}^{**}}-{{u}^{*}}}{\vartriangle t}=\frac{\mu }{\rho }{{\triangledown }^{2}}u+F \\ \frac{{{\left( {{\partial }_{i}}u \right)}^{**}}-{{\left( {{\partial }_{i}}u \right)}^{*}}}{\vartriangle t}= \\ -{{\partial }_{i}}{{u}^{*}}\triangledown \cdot {{u}^{*}}+{{\partial }_{i}}\left( \frac{\mu }{\rho }{{\triangledown }^{2}}u+F \right) \\ \end{matrix} \right. $

(14)

其中,方程(14)的空间离散采用中心差分格式。

3)非对流项(Ⅱ):

压力和速度的匹配:

$ \left\{ \begin{matrix} \frac{\partial u}{\partial t}=-\frac{1}{\rho }\triangledown p \\ \frac{\partial \left( {{\partial }_{i}}u \right)}{\partial t}=-{{\partial }_{i}}\left( \frac{1}{\rho }\triangledown p \right) \\ \end{matrix} \right. $

(15)

通过对方程 (15)第一行公式取散度，并引入连续方程，可得到如下形式的泊松方程：

$ \triangledown \cdot \left( \frac{1}{\rho }\triangledown {{p}^{n+1}} \right)=\frac{1}{\vartriangle t}\triangledown \cdot {{u}^{**}} $

(16)

泊松方程的求解通过SOR迭代得到。

考虑动量方程的压力梯度项，计算速度的最终值，计算式为

$ \left\{ \begin{matrix} {{u}^{n+1}}={{u}^{**}}-\frac{\vartriangle t}{\rho }\triangledown {{p}^{n+1}} \\ {{\left( {{\partial }_{i}}u \right)}^{n+1}}={{\left( {{\partial }_{i}}u \right)}^{**}}-\vartriangle t{{\partial }_{i}}\left( \frac{1}{\rho }\triangledown {{p}^{n+1}} \right) \\ \end{matrix} \right. $

(17)

根据式(4)、(5)和(18)更新界面和网格内流体信息，然后返回到步骤1，这样就完成了整个计算过程，一直到设定的步骤结束。

1.2.3 固体边界处理

采用浸入边界法(immersed boundary method)^[14]处理固体对流场的影响:

$ {u^{n + 1}} = {\varphi _2}{U_b}^{n + 1} + \left( {1 - {\varphi _2}} \right){u_f}^{n + 1} $

(18)

式中：u_fⁿ⁺¹为在欧拉坐标体系中通过式(17)计算得到的“局部”流场速度；U_bⁿ⁺¹为通过拉格朗日方法得到的固体速度^[15]，本文中柱体静止，固体的速度设定为零；uⁿ⁺¹为二者耦合作用的“整体”速度，这样就得到了流-固耦合流场的求解。上述方法也可称为固体体积流-固耦合处理方法。基于浸入边界方法的流-固耦合技术，采用固定网格，在计算过程中无需更新网格，可保证模型的高效率，使得模型在处理大位移的流固耦合计算中具有强大优越性。本文采用的方法具有较易实施的优势，方便模型向三维拓展。

1.3 参数定义

升力F_l为

$ {{F}_{i}}={{F}_{i}}^{\left( p \right)}+{{F}_{i}}^{\left( v \right)}=\mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}\limits_A \frac{\partial p}{\partial {{x}_{i}}} \text{d}A+\mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}\limits_A \frac{\partial \left( 2\mu {{S}_{ij}} \right)}{\partial {{x}_{i}}} \text{d}A $

(19)

式中：Re为雷诺数,C_d为阻力系数,F_d为阻力,C_l为升力系数,S_t为Strouhal数,上角标p和v分别表示由于压力和摩擦力所产生的阻力或升力，S_ij=(əu_i/əx_j+əu_j/əx_i)/2。

1.4 数值模型

计算区域及网格划分如图 3和图 4所示。圆柱与方柱的特征长度均为D=1 m，上游边界距离柱体中心7.5D，下游边界距离柱体中心22.5D，两侧边界距离柱体中心均为7.5D。为了捕捉柱体周围的复杂流动，柱体近壁面采用细网格。

1.5 边界条件和初始条件设置

计算流场的边界条件以及初始条件设置如下：初始条件：初始设定一个从左向右的速度场，u=U₀，v=0；其中，u、v分别为x、y方向的速度；初始压力场都设定为零。边界条件：来流边界保持为均匀来流速度U₀；出口边界条件，开边界；柱体表面边界条件：无滑移边界，即u_n=0；两侧壁面：自由滑移边界。

2 模型验证及数据分析 2.1 网格收敛性分析

为验证模型的有效性及网格收敛性，采用三套网格对Re=100时圆柱绕流进行模拟。三种网格从网格1到网格3，最小网格分别为0.04、0.02和0.01 m，模拟结果如表 1所示。由于网格1的网格尺寸较大，C_d值与文献结果相差较为明显，而网格2与网格3误差较小。考虑到计算效率和精度，选择网格2对下文其他工况进行模拟。

2.2 圆柱绕流参数分析

下面对Re=1~300的圆柱绕流进行模拟，结果如表 2所示。从表中可看出，Re<50时，没有漩涡脱落，S_t数为零。

图 5给出了圆柱绕流主要水动力参数随Re的变化情况。在Re<40时，没有漩涡脱落为层流状态，如图 5(a)所示在开始段C_d骤减，主要是因为4<Re<5^[19]，在圆柱的后方开始形成对称的漩涡，而漩涡处在负压区，从而导致C_d骤减，待漩涡稳定附着在圆柱后方时，C_d也逐渐稳定。C_d^(p)(压力阻力项)和C_d^(v)(摩擦阻力项)的变化趋势与C_d类似；在Re=1时C_d^(v)占主导地位，随着雷诺数的增大，C_d^(p)所占比例逐渐增加，2条直线趋于平行，与文献结果吻合较好。在50<Re<300时，圆柱后方有卡门涡街产生，层流转化为非层流；C_d变化如图 5(b)所示，先下降后有所上升，在Re≈150时达到极小值，此处可理解为三维影响的起始点；从图 5(c)可看出在非层流状态下，C_d^(p)已完全占据主导地位。由于层流状态下没有漩涡脱落，所以不存在升力，对应的C_l为零，由图 5(d)所示，非层流状态下，C_l随着的雷诺数的增大而增加，间接体现出漩涡的能量随着雷诺数的增加而增加，与文献吻合较好。

图 6给出了50<Re<300时，Strouhal数和主波长λ随Re的变化曲线。S_t数随着Re的增大而增大，呈现单调特性，与Rajani等^[3]的结果吻合较好，与Williamson^[16]实验的结果在Re<100吻合较好；Re>100时由于考虑到实际情况的三维效应，数模的二维结果会高估S_t数。主波长λ为卡门涡街一侧前后2个涡之间的距离，如图 7所示，λ与S_t数的关系可近似表示为λ/D=1/S_t^[22]。由图可知随着Re的增大，λ呈单调下降趋势，卡门涡街漩涡之间的距离越来越小。

2.3 方柱绕流参数分析

下面将对Re<300情况下的方柱绕流进行模拟，结果如表 3所示。与圆柱绕流结果类似，本文模拟出的结果在Re<50时显示没有漩涡脱落，依旧为层流，S_t数值为零；而在Re=50时，不稳定状态首次显现，伴随有漩涡脱落，流动转化为非层流，说明在Re=40~50，漩涡开始脱落，与Sohankar等^[23]通过数值模拟测得的Re=52较为接近。

与圆柱绕流类似，在Re<40时，是层流状态。图 8(a)给出了C_d、C_d^(p)、C_d^(v)的变化，可发现整体趋势相同，最后近似发展为3条平行的直线，与Sen等^[7]在2<Re<40的结果吻合较好。图 8(b)给出了平均阻力系数随Re的变化，从趋势来看与文献结果基本一致，从数值来看，在Re<100时，与其他文献吻合较好，Re>100时，与毕继红等^[25]的结果相差较小，与其他文献相差较大，其中的影响因素较多，像离散方法不同、计算域的大小、边界条件设置的不同均可导致计算结果的误差。在Re≈140时，阻力系数达到极小值，此时分离点由后角点处转移到前角点处(本文将此时定义为“转折点”，方便后文叙述)，使得方柱侧面受力发生变化，导致C_d^(v)的合力方向发生变化，在Re<140时，C_d^(v)数值为正；而如表 4所示，在Re>140时，C_d^(v)数值为负，方向发生了变化，这也是方柱绕流与圆柱绕流最大的区别之一。由图 8(c)可看出方柱绕流中，压力阻力项仍然占有主导地位，且随着雷诺数的增大，摩擦阻力项相比于压力阻力项基本可以忽略。从图 8(d)中可以看出在Re≈140时，升力系数的增长率逐渐减小，这也是分离点变化导致的结果，过了转折点之后，升力系数的增长率变大，结果与现有文献吻合良好。从图 8(e)可看出S_t数的变化趋势与现有文献基本一致，在转折点处达到最大值，而数值方面的误差可能由多种原因导致，Breuer等^[6]指出阻塞比的增加会导致S_t数的增大。主波长λ的变化与S_t数相反，先减小再增大，在转折点处的主波长达到最小。

2.4 圆柱与方柱绕流参数对比

为了清晰说明圆柱绕流与方柱绕流的区别，将结果中的重要参数进行对比。

图 9给出了圆柱、方柱绕流主要水动力参数随Re的变化曲线。由图 9(a)可以看出，C_dmean随Re的变化趋势大体一致，在Re<10时，C_dmean随Re的变大急剧减小，然后趋于平稳，且方柱的阻力系数大于圆柱，这是两者形状不同所导致的。两者在Re≈140~150达到极小值，此时方柱的分离点由后角点转移到前角点，方柱摩擦阻力的方向发生变化，而圆柱没有此现象发生，其分离点只围绕一点往复运动，且摩擦阻力保持为正值，方向不变。图 9(b)给出了圆柱和方柱C_lmean随Re的变化曲线。在Re<150时，两者的大小基本一致，经过转折点之后，方柱的升力系数显著增大，大体呈线性增长，分离点在前角点时，方柱两侧面将受到附面涡的作用，使得升力增大，而圆柱不存在此现象，升力系数变化相对平稳。图 9(c)给出了圆柱和方柱S_t随Re的变化曲线，圆柱的S_t随着Re的增大呈单调上升趋势，而圆柱的S_t先增加，在转折点处到达极大值，后随着Re的增大而减小。方柱的S_t值一直小于圆柱的S_t值，这也说明方柱的漩涡脱落频率比圆柱的漩涡脱落频率小。图 9(d)给出了圆柱和方柱λ/D随Re的变化曲线。可看出，圆柱绕流的主波长随着Re的增大而单调下降，而方柱绕流的主波长先随着Re的增大先下降再上升，在转折点附近到达极小值，且方柱绕流的主波长一直大于圆柱绕流的主波长。

图 10(a)、(b)给出了圆柱、方柱阻力和升力系数振幅随Re的变化。可看出，升力系数的振幅比其对应的阻力系数的振幅大一个数量级，两者均随着Re的增大而增大，且在方柱绕流的转折点之前，圆柱和方柱的阻力升力系数振幅大体一致，转折点之后，方柱的阻力升力系数振幅显著增大，圆柱的变化依旧平缓，且方柱的振幅大于圆柱的振幅。图 10(c)给出了C_d^(p)/C_d^(v)随Re的变化曲线，可看出C_d^(p)相比于C_d^(v)一直处于主导地位，且随着Re的增大越明显。圆柱绕流C_d^(p)/C_d^(v)随Re呈线性增长，而方柱绕流C_d^(p)/C_d^(v)随Re呈非线性增长，随着Re的增加，变化率越来越大，且在方柱绕流中，C_d^(p)的主导地位相对圆柱绕流要明显。

3 结 论

本文利用自主研发的CIP-ZJU模型，对低雷诺数(Re<300)圆柱和方柱绕流问题开展了数值模拟。模型在直角坐标系统下建立，采用紧致插值曲线CIP方法作为流场的基本求解器，离散了Navier-Stokes方程,利用浸入边界方法处理固体边界，基于多相流的理论实现流-固耦合同步求解。分别分析了圆柱和方柱绕流的绕流水动力特性，并将两者的主要水动力参数进行比较，得出如下结论：

1)圆柱和方柱绕流的C_d在Re≈140~150达到极小值，此时方柱的分离点由后角点转移到前角点，方柱的摩擦阻力的方向发生变化，而圆柱没有此现象发生，其分离点只围绕一点往复运动，且摩擦阻力保持为正值。

2)圆柱绕流的C_d^(p)、C_d^(v)均为正值，而方柱绕流在转折点之后，由于分离点的变化，导致摩擦阻力的合力方向发生改变，C_d^(v)变为负值。

3)方柱绕流的S_t值一直小于圆柱绕流的S_t值，说明方柱的漩涡脱落频率比圆柱的漩涡脱落频率小。相反方柱绕流的主波长大于圆柱绕流的主波长。

4)圆柱绕流与方柱绕流升力系数的振幅比其对应的阻力系数的振幅大一个数量级，两者均随着Re的增大而增大，且在方柱绕流的转折点之前，圆柱和方柱的阻力升力系数振幅大体一致，转折点之后，方柱的阻力升力系数振幅显著增大，圆柱的变化依旧平缓，且方柱的振幅大于圆柱的振幅。

本文的计算结果与文献结果吻合良好，验证了此方法的可行性，为进一步探讨柱体形状对绕流阻力的影响，开展流致振动等研究打下良好基础。