阻尼比对串列双圆柱二自由度涡激振动及传热影响的数值研究

Numerical Study on the Effect of Damping Ratio on Vortex-Induced Vibration and Heat Transfer of Two Degrees of Freedom Double Cylinders in Series

摘要: 采用数值方法研究了流体流过刚性串列双圆柱产生涡激振动(VIV)时的流动和传热。研究对象是二维的，雷诺数Re = 150，双圆柱的间距比为8，且都具有二平动自由度(2DOF)。采用有限体积法求解Navier-Stokes方程获得流场，同时用四阶龙格–库塔算法耦合求解圆柱的运动。主要研究了在不同折合速度Ur下，不同阻尼比(ζ = 0, 0.01, 0.1)时对上下游二自由度(2DOF)圆柱涡激振动时传热的影响。数值计算结果表明：当2 < Ur ≤ 6时，随着折合速度Ur的增大，在不同阻尼比时，上游圆柱的时间平均努塞尔数均先增大后减小。当Ur = 5时，阻尼比从ζ = 0增大为ζ = 0.1时，下游圆柱的时间平均努塞尔数减小8.74%。当折合速度Ur > 5后，下游圆柱的运动轨迹受上游圆柱尾迹区形成的波动的影响，下游圆柱的运动轨迹呈现不规则变化。

Abstract: The flow and heat transfer of vortex-induced vibration (VIV) generated by fluid flowing through a rigid tandem double cylinder are studied numerically. The Reynolds number is 150 and the clear-ance ratio of the two cylinders is 8, and both of them have two degree of freedom (2DOF). The finite volume method is used to solve Navier-Stokes equations to obtain the flow field, and the fourth order Runge-Kutta algorithm is used to solve the cylinder motion. The influence of different damping ratios (ζ = 0, 0.01, 0.1) on the heat transfer of a 2DOF cylinder under vortex-induced vibration at different reduced velocities Ur is studied. The numerical results show that when 2 < Ur ≤ 6, with the increase of reduced velocity Ur, the time mean Nusselt number of upstream cylinder increases first and then decreases with different damping ratios. When Ur = 5, the damping ratio increases from ζ = 0 to ζ = 0.1, and the time mean Nusselt number of downstream cylinder decreases by 8.74%. When the reduced velocity Ur > 5, the motion trajectory of the downstream cylinder is affected by the fluctuation of the wake region of the upstream cylinder, and the motion trajectory of the down-stream cylinder shows irregular changes.

1. 引言

将圆柱置于一定的自由速度的来流中，圆柱两侧会发生交替涡脱落，圆柱因此会受到流向(x方向)和横向(y方向)的脉动压力。当圆柱为弹性支撑时，脱落的涡会引发圆柱振动，圆柱的运动反过来会改变尾流结构，这种流体与圆柱结构之间相互作用被称为涡激振动(Vortex-Induced vibration, VIV)，有的文献中也称之为流致振动。当涡脱落频率接近圆柱结构的固有频率时，会出现共振或锁定现象，导致圆柱振荡显著 [1]。然而，在许多工程领域(如热交换器、海洋立管等)会发生圆柱结构的涡激振动。近几年来，这一现象吸引了大量研究者对其进行了广泛的研究 [2] [3] [4] [5] [6]。

关于圆柱涡激振动的研究可以分为流动特性和传热特性这两个方面。在研究双圆柱外涡激振动的流动特性方面。Zhao等 [7] 数值研究了两个不同直径的刚性圆柱在雷诺数为250时的二自由度涡激振动。在一定的直径比和质量比的条件下，研究了小圆柱位置角α对涡激振动时旋涡频率锁定区的影响。当小圆柱的位置角为α = 0˚，22.5˚，90˚，112.5˚时，与单圆柱相比，折合速度U_r的锁定区的范围明显扩大。Nepali等 [8] 对低雷诺数下串列两方柱二自由度涡激振动进行研究。在低Re时观察到典型的8字形轨迹，当Re较高时，由于尾迹干扰的影响，下游方柱的运动轨迹会变得不规则。

在研究圆柱涡激振动的传热特性方面，Izadpanah等 [1] 研究了单自由度涡激振动对圆柱传热的影响。折合速度U_r = 4和阻尼比ζ = 0时传热最强。Yang等 [9] 对等温圆柱的二自由度流体诱导振动与传热进行了数值模拟。在研究的折合速度下，可以观察到2S (单)涡型和典型的“8”形轨迹。在VIV中，圆柱的振幅和表面平均努赛尔数随时间呈周期性变化。由于圆柱的振动，最大局部努塞尔数所对应的位置偏离了前驻点。折合速度U_r的变化会影响涡的形成和脱落，从而影响圆柱的振动响应和传热。

综上所述，尽管已有大量文献对圆柱涡激振动和相应的流体流动进行了研究，但是只有少数学者研究了圆柱涡激振动对传热的影响。在工程实践中，相较于固定圆柱(换热圆管)而言，圆柱的运动(振动)势必改变圆柱与流体之间的传热特性。另外，圆柱的振动又受到其结构阻尼系数的影响。所以，本文研究在不同阻尼比(ζ = 0, 0.01, 0.1)下串列双圆柱二自由度涡激振动时的对流传热特性。所考虑的流体为空气(Pr = 0.7)，考察了在不同阻尼比和不同折合速度时，流体的流动、圆柱的振动及流体与圆柱之间的对流传热特性。

2 数值模拟方法

2.1. 物理模型和控制方程

图1中为本文考虑的二维流固耦合运动和传热计算模型。如图1所示，圆柱的直径为D、质量为m，圆柱为有结构阻尼的弹性安装，弹性系数为k，阻尼比为ζ，双圆柱的参数均保持一致。温度为T_∞的低温流体由左侧入口流入计算区域后，流经表面维持恒定高温T_h的圆柱，并与之换热。根据Zhao等 [5] 和Sun等 [10] 的研究，计算区域的尺寸取为50D × 25D，此时上、下边界和右侧出口对圆柱运动及其附近流体流动影响可以忽略。入口边界到圆柱中心的距离为12.5D，两圆柱之间距离为8D出口边界到圆柱中心的距离为29.5D，上下边界到圆柱中心的距离各为12.5D。图1中的两个圆柱都具有x和y的两个方向的平动自由度。

Figure 1. Computational model

图1. 计算模型

假设流体为常物性、不可压缩的牛顿流体。忽略粘性耗散后，二维非定常无量纲控制方程(Navier-Stokes方程)为 [3]：

$\frac{\partial U}{\partial \tau }+U\frac{\partial U}{\partial X}+V\frac{\partial U}{\partial Y}=-\frac{\partial P}{\partial X}+\frac{U_r}{\mathrm{Re}}\left(\frac{{\partial }^{2}U}{\partial X^2}+\frac{{\partial }^{2}U}{\partial Y^2}\right)$ (1)

$\frac{\partial V}{\partial \tau }+U\frac{\partial V}{\partial X}+V\frac{\partial V}{\partial Y}=-\frac{\partial P}{\partial X}+\frac{U_r}{\mathrm{Re}}\left(\frac{{\partial }^{2}V}{\partial X^2}+\frac{{\partial }^{2}V}{\partial Y^2}\right)$ (2)

$\frac{\partial \Theta }{\partial \tau }+U\frac{\partial \Theta }{\partial X}+V\frac{\partial \Theta }{\partial Y}=\frac{U_r}{\mathrm{Pr}\mathrm{Re}}\left(\frac{{\partial }^2\Theta }{\partial X^2}+\frac{{\partial }^2\Theta }{\partial Y^2}\right)$ (3)

其中： $X=x/D$， $Y=y/D$ 为无量纲坐标，时间 $\tau =f_{n}t$， $f_n=\frac{1}{2}\text{π}\sqrt{k/m}$ 为圆柱固有频率， $U_r=U_{\infty }/f_{n}D$ 为折合速度，速度 $U=\frac{u}{f_{n}D}$， $V=\frac{v}{f_{n}D}$，压强 $P=\frac{p}{\rho f_n^2{D}^2}$，温度 $\Theta =\frac{T-T_{\infty }}{T_h-T_{\infty }}$。雷诺数 $\mathrm{Re}=U_{\infty }D/\nu$，普朗特数 $\mathrm{Pr}=\nu /a$， $\nu$ 为运动粘度， $a$ 为热扩散系数。

圆柱运动的无量纲方程为：

$\frac{{\text{d}}^2{X}^{*}}{\text{d}{\tau }^2}+4\text{π}\zeta \frac{\text{d}{X}^{*}}{\text{d}\tau }+4{\text{π}}^2{X}^{*}=\frac{2}{\text{π}}\frac{U_r^2{C}_D}{m^{*}}$ (4)

$\frac{{\text{d}}^2{Y}^{*}}{\text{d}{\tau }^2}+4\text{π}\zeta \frac{\text{d}{Y}^{*}}{\text{d}\tau }+4{\text{π}}^2{Y}^{*}=\frac{2}{\text{π}}\frac{U_r^2{C}_L}{m^{*}}$ (5)

其中， $X^{*}=x^{*}/D$ 和 $Y^{*}=y^{*}/D$ 分别为圆柱在x和y方向上的位移， $m^{*}=m/m_d$ 为质量比，m为圆柱的质量， $m_d=\frac{1}{4}{\rho }_f\text{π}{D}^2$ 为圆柱排开的流体的质量。 $C_D=f_D\left(t\right)/\frac{1}{2}\rho U_{\infty }^{2}D$，其中 $C_D$ 为阻力系数， $f_D\left(t\right)$ 为圆柱所受阻力(x方向)，升力系数 $C_L=f_L\left(t\right)/\frac{1}{2}\rho U_{\infty }^{2}D$，其中 $f_L\left(t\right)$ 为圆柱所受升力(y方向)。

本文以空气为流动工质，入口速度恒定，出口为恒压力出口，Pr = 0.71，雷诺数Re为150。圆柱表面为无滑移固体壁面，质量比m^* = 3。上下两边界为对称边界条件，边界上V = 0，其余各变量梯度沿法向方向的分量为零。流体初始速度和压力在整个计算域内为零，温度为低温，圆柱初始时刻在平衡位置处于静止。以上所述边界条件和初始条件对应的数学表达为：

入口： $U=U_r$， $V=0$， $\Theta =0$

上边界和下边界： $\frac{\partial U}{\partial Y}=0$， $V=0$， $\frac{\partial \Theta }{\partial Y}=0$

圆柱壁面： $U=0$， $V=0$， $\Theta =1$

出口： $\frac{\partial V}{\partial X}=0$， $P=0$， $\frac{\partial \Theta }{\partial X}=0$

初始时刻： $U=0$， $V=0$， $P=0$， $\Theta =0$ ； $X^{*}=0$， $Y^{*}=0$

2.2. 数值方法

本文中的流动采用了SIMPLE算法求解Navier-Stokes方程。对流项的离散采用QUICK式，瞬态项采用了二阶隐式离散格式 [11]。固体运动的控制方程采用四阶龙格–库塔法求解。采用了动网格技术，通过固体表面无滑移及流体和固体之间作用力相等实现流体流动和固体运动的耦合计算。无量纲时间步长为0.005。在一个时间步长内，Nu_s相对偏差在10⁻⁵以下时，认为迭代在该时间步内收敛。

局部努塞尔数Nu的定义为：

$\text{Nu}=\frac{h\cdot D}{\lambda }$ (6)

其中， $h=-\frac{\lambda }{D}\frac{\partial \Theta }{\partial n}$ 为对流换热系数， $\lambda$ 为导热系数，n为圆柱表面外法线方向。

表面平均努塞尔数(Nu_s)定义为：

${\mathrm{Nu}}_s=\frac{1}{2\text{π}}{\displaystyle \underset{0}{\overset{2\text{π}}{\int }}\mathrm{Nu}\text{d}\theta }$ (7)

时间平均努塞尔数(Nu_A)定义为：

${\mathrm{Nu}}_A=\frac{1}{{\tau }_c}{\displaystyle \underset{0}{\overset{{\tau }_c}{\int }}{\mathrm{Nu}}_s\text{d}\tau }$ (8)

其中 ${\tau }_c$ 是一个稳定周期。

圆柱最大横向振幅 $A_y/D=\left(y_{\max }^{\ast }-y_{\min }^{\ast }\right)/2D$，最大流向振幅 $A_x/D=\left(x_{\max }^{\ast }-x_{\min }^{\ast }\right)/2D$。

2.3. 网格划分与独立性验证

如图2所示，本文采用了非均分的多块结构化四边形网格。在圆柱壁面附近设置了较密集的网格，使网格尺寸平稳地过渡，以分辨其此处较剧烈的速度和温度变化。选用了四套网格密度不同的设置，在Re = 150，m^* = 3，ζ = 0，U_r = 7时的工况下进行网格独立性测试。按照网格总数由小到大(83,600、114,475、150,932、189,528)依次编号为G1~G4。

(a) 在计算域内的网格分布 (b) 初始位置的放大图

Figure 2. Cylindrical grid

图2. 圆柱网格

网格独立性的计算结果列于表1，括号内是以网格G4为基准的相对偏差。从表中可以看出，从G1到G4，随着网格数的增加，最大振幅A_y/D和局部努塞尔数Nu的变化趋缓。当网格数大于G3的设置后，偏差小于2%，数值计算结果已经满足网格独立性的要求。综合考虑计算精度和效率，本文所有工况选取G3的网格设置进行计算，生成的计算网格图2中图2(a)及其放大图2(b)。

Table 1. Grid independence study when Re = 150, m* = 3, ζ = 0, Ur = 7

表1. Re = 150，m^* = 3，ζ = 0，U_r = 7时网格独立性研究

2.4. 数值方法的验证

为了验证本文所采用数值方法的正确性，对一自由度单圆柱涡激振动问题进行了数值模拟，并与采用了不同数值方法的文献 [3] [12] [13] 的结果进行了对比。所考虑的雷诺数Re = 150，质量m^* = 2.55，无量纲阻尼比ζ = 0，折合速度U_r在3.0~8.0范围内变化。

本文和其他研究者计算所得的圆柱最大振幅随U_r的变化曲线绘于如图3。与其他文献相一致地，在所研究的参数范围内，振幅的最大值均未超过圆柱的直径。在U_r = 4.0~7.0时，出现了“锁定”现象，圆柱振幅最大值为A_y/D = 0.39~0.57。在该区域以外，圆柱的最大振幅约为Ay/D = 0.07。图中可以看出，本文的计算结果与已有文献很好地吻合。此外，本文和文献使用了不同的数值方法，证明了本文所采用的数值方法和计算结果的可靠性和正确性。

Figure 3. Maximum amplitude response curve

图3. 最大振幅响应曲线

3. 结果分析

3.1. 涡激振动和流动特性

图4中图4(a)、图4(b)分别为上下游圆柱在不同阻尼比ζ和不同U_r下的横向(y)和流向(x)振幅值。从图4(a)中可以看出，上游圆柱在不同阻尼比ζ下的横向振幅值随U_r的增大先增大后减小，并且上游圆柱在三种不同阻尼比ζ下均在U_r = 3时取得横向振幅最大值。在本文研究的U_r范围内，对于某一特定的U_r下，随着阻尼比的增大，上游圆柱的振幅值逐渐减小。下游圆柱在不同阻尼比ζ下的横向振幅值也随U_r的增大先增大后减小，但是，阻尼比ζ = 0和阻尼比ζ = 0.01的曲线均在U_r = 5时取得横向振幅最大值，而阻尼比ζ = 0.1的曲线在U_r = 4时取得横向振幅最大值，同时，在图4(b)中，阻尼比ζ = 0.1的曲线在U_r = 4时，下游圆柱取得流向振幅最大值。当阻尼比增大时，圆柱在横向和流向方向上，较大的U_r下会形成较大的阻力，使得振幅减小。下游圆柱处在上游圆柱尾迹涡流中，其流动特性也会受上游圆柱的影响。图4(b)中，上游圆柱在不同阻尼比ζ下的流向振幅值随U_r的增大先增大后减小，并且上游圆柱在三种不同阻尼比ζ下均在U_r = 3时取得流向振幅最大值。由于下游圆柱处在上游圆柱尾迹涡流中，其流动特性也会受上游圆柱的影响，所以，下游圆柱的流向振幅在不同阻尼比ζ下呈现出不同的变化形式。

(a) 横向(y)振幅(b) 流向(x)振幅

Figure 4. Transverse (y) and streamwise (x) amplitudes under different U_r

图4. 不同U_r下的横向(y)和流向(x)振幅

为了显示涡激振动对圆柱周围流场影响，图5给出了U_r = 1、U_r = 3、U_r = 5、U_r = 6、U_r = 8、U_r = 9在不同阻尼比下的瞬时涡量等值线。

如图5所示，在三种不同阻尼比和不同折合速度U_r下，上下游圆柱在每个振动周期内从圆柱体上脱落的两个单独的涡，两个涡的强度几乎相等，但旋转方向相反，这样涡脱落的模式称为2S。如图中所示，圆柱在正方向和负方向的运动幅度是相等的，圆柱运动的平均值为零。随着折合速度增大，涡旋形成长度增加。从图中还可以看出，上游圆柱在每个振动周期内从圆柱体上脱落的两个单独的涡直接又作用在下游圆柱上，这样使得下游圆柱脱落的涡和上游圆柱的脱落涡在尾部形成了很大区域的涡流，进而影响了下游圆柱的传热过程。因为每个涡旋都会带走一定的热量，所以对流传热直接受到涡脱落的影响。当上下游圆柱都具有二自由度(2DOF)振动时，上游圆柱涡脱落带走的热量又返回到了下游圆柱，导致下游圆柱的传热过程受影响。交替的涡脱落导致圆柱上出现周期性流体力，从而引起圆柱周期性振动，即VIV。同时，由于涡脱落和圆柱的周期性运动导致了Nu_S的周期性变化，这也使得在不同阻尼比和不同折合速度U_r下，上下游圆柱的时间平均努塞尔数有差别。

Figure 5. The instantaneous vorticity contour at (a) U_r = 1, (b) U_r = 3, (c) U_r = 5, (d) U_r = 6, (e) U_r = 8, (f) U_r = 9

图5. 瞬时涡量等值线：(a) U_r = 1，(b) U_r = 3，(c) U_r = 5，(d) U_r = 6，(e) U_r = 8，(f) U_r = 9

圆柱的振动响应和传热特性随折合速度的变化而变化。图6给出了上下游圆柱在阻尼比ζ = 0时，不同折合速度U_r下的运动轨迹图。从图中可以看出，在本文研究的折合速度范围内，除了折合速度U_r = 6时，上游圆柱的运动轨迹呈现不规则变化，在其余折合速度下，上游圆柱的运动轨迹呈现“8”形轨迹，而此轨迹与Rabiee等 [14] 和Kang等 [15] 所得结论类似。“8”形轨迹也说明圆柱在x方向的振动频率是y方向的2倍，这是因为，圆柱后的尾涡是上、下交替脱落的，涡脱落时对圆柱在流向上的振动激励总是沿着x轴负方向，而同时发生在横向上的振动激励则是沿着y方向正负交替的。而当折合速度U_r = 6时，此时上游圆柱的流向振幅值达到最大，横向振幅值也较大，使得圆柱振动比较剧烈，而运动轨迹则呈现不规则变化。

在折合速度U_r = 1和U_r = 2，时，此时上游圆柱振动幅度较小，使得横向振幅值和流向振幅值也较小，尾迹区形成的波动对下游圆柱的影响也较小，所以，下游圆柱在较小的折合速度U_r下，圆柱的运动轨迹呈现“8”形轨迹。当折合速度U_r = 3时，上游圆柱振动相对比较剧烈，圆柱的横向振幅值达到最大，使得尾迹区形成的较大的波动，而此时下游圆柱的横向振幅值较小，所以下游圆柱的运动轨迹受上游圆柱振动情况的影响比较显著，下游圆柱的运动轨迹呈现不规则变化。当折合速度U_r = 5时，上游圆柱的振动幅度相对减弱，而此时下游圆柱自身的振动幅度较大，横向振幅值达到最大，所以下游圆柱的运动轨迹受上游圆柱振动情况的影响较小，下游圆柱的运动轨迹呈现“8”形轨迹。当折合速度U_r > 5后，下游圆柱的运动轨迹受上游圆柱尾迹区形成的波动的影响，下游圆柱的运动轨迹呈现不规则变化。圆柱在流向方向(x方向)运动时，随着折合速度U_r增加，圆柱的平衡位置明显地向下游(x + 向)移动。对于圆柱在横向方向(y方向)的振动，圆柱的平衡位置在y = 0处。x方向位移和y方向位移的相位差会影响圆柱运动轨迹的方向和形状。

(a) (b)(c) (d)

Figure 6. The trajectory of upstream and downstream cylinders under different U_r when damping ratio ζ = 0: (a), (c) upstream, (b), (d) downstream

图6. 阻尼比ζ = 0时不同U_r下上游和下游圆柱的轨迹图：(a)，(c)为上游，(b)，(d)为下游

3.2. 传热特性

如图7所示，图为上下游圆柱在不同阻尼比下，时间平均努塞尔数随U_r的变化情况。从图中可以看出，阻尼比对传热有显著影响。当U_r ≤ 2时，在不同阻尼比时，上下游圆柱的时间平均努塞尔数相差都很小，表明了在此范围内折合速度下，不同阻尼比对圆柱的传热影响较小。当2 < U_r ≤ 6时，随着折合速度U_r的增大，上游圆柱在不同阻尼比时，时间平均努塞尔数的差值先增大后减小，说明在此范围内折合速度下，不同阻尼比对上游圆柱的传热影响先增强后逐渐减弱，当在U_r = 3时，阻尼比为ζ = 0时比阻尼比为ζ = 0.1时的时间平均努塞尔数大3.33%，说明在此折合速度时，不同阻尼比对上游圆柱的传热影响最显著。当6 < U_r ≤ 9时，上游圆柱在不同阻尼比时的时间平均努塞尔数的差值都很小，说明在此范围内折合速度下，不同阻尼比对圆柱的传热影响较小。当2 < U_r ≤ 5时，随着折合速度U_r的增大，下游圆柱在不同阻尼比时，时间平均努塞尔数的差值越明显，而当在U_r = 5时，阻尼比为ζ = 0时比阻尼比为ζ = 0.1时的时间平均努塞尔数大8.74%，说明在此折合速度时，不同阻尼比对下游圆柱的传热影响最显著。当5 < U_r ≤ 9时，随着折合速度U_r的增大，下游圆柱在不同阻尼比时的时间平均努塞尔数的差值逐渐减小，说明在此范围内折合速度下，不同阻尼比对圆柱的传热影响也逐渐减弱。

Figure 7. Change of time average Nusselt number with U_r

图7. 时间平均努塞尔数随U_r的变化情况

为了显示涡激振动对圆柱周围热场的影响，图8给出了U_r = 1、U_r = 3、U_r = 5、U_r = 6、U_r = 8、U_r = 9在不同阻尼比下的温度等值线。

如图8所示，从图中的温度分布上可以看出流体与VIV圆柱之间的换热过程。冷流体冲刷热上游圆柱的上游侧，经过加热后，一部分热流体流向了下游圆柱，另外一部分热流体了流向尾迹区。由于上下游圆柱的振动并不是保持同步的，而驻点的位置也会发生变化。同时，流动剪切层的分离点在上下游圆柱表面移动，这与来流速度和圆柱的VIV响应有关。剪切层的驻点和分离点都改变了热边界层的发展，影响了上下游圆柱与流体之间的传热。流动剪切层的分离、涡脱落、结构振动和热边界层的演化相互作用，实现了流动、结构振动和传热的耦合。结合图5所示的尾迹涡结构可知，尾部停滞点附近的流动特性是影响上下游圆柱与流体传热的关键因素。

Figure 8. The temperature contour at U_r = 1, (b) U_r = 3, (c) U_r = 5, (d) U_r = 6, (e) U_r = 8, (f) U_r = 9

图8. 温度等值线：(a) U_r = 1，(b) U_r = 3，(c) U_r = 5，(d) U_r = 6，(e) U_r = 8，(f) U_r = 9

从图5和图8可以看出温度等值线的变化与涡量等值线的变化相似，这是由于涡对圆柱面附近热流体的输运作用。当流动剪切层被拉伸时，热边界层也发生相同的变化。来流速度对涡旋的形成和脱落有显著影响，从而影响热边界的厚度。

图9给出了阻尼比ζ = 0时，上下游圆柱在折合速度U_r = 3和U_r = 5时一个稳定周期T内的局部努塞尔数沿圆柱壁面分布。从图中图9(a)、图9(b)明显可以看出，当U_r = 3时，上游圆柱的局部努塞尔数大于下游圆柱的局部努塞尔数，这与图7中的时间平均努塞尔数分布图很好的对应。当U_r = 3时，从图9(a)中还可以看出，当最开始0/4T周期时，下游圆柱局部努塞尔数最大值出现在前滞止点靠上的位置，之后，随着时间的改变，局部努塞尔数最大值出现的位置发生了改变，当下游圆柱振动到3/4T周期时，下游圆柱局部努塞尔数最大值出现在前滞止点靠下的位置，从图9(b)中也可以得出类似结论。所以，上下游圆柱局部努塞尔数最大值出现的位置发生在前滞止点靠上或者靠下位置，一般呈对称分布。

当U_r = 5时，从图9(c)、图9(d)中可以看出，在某些时刻，上下游圆柱的局部努塞尔数最大值相差很小。但是，图9(c)中在一个稳定周期T内局部努塞尔数曲线所包围的面积值大于图9(d)中对应一个稳定周期T内局部努塞尔数曲线所包围的面积值。分析得出上游圆柱的时间平均努塞尔数大于下游圆柱的时间平均努塞尔数，此结论与图7中时间平均努塞尔数分布图很好的对应。图9(c)、图9(d)中上下游圆柱的局部努塞尔数最大值也是出现在前滞止点靠上或者靠下的位置，并且呈对称分布。

(a) 上游圆柱U_r = 3 (b) 下游圆柱U_r = 3 (c) 上游圆柱U_r = 5 (b) 下游圆柱U_r = 5

Figure 9. The distribution of the local Nusselt number of the upstream and downstream cylinders at ζ = 0

图9. ζ = 0时上下游圆柱局部努塞尔数的分布情况

图10给出阻尼比ζ = 0时，上下游圆柱在U_r = 3和U_r = 5时温度等值线的局部放大图。从图10(a)中可以看出，在U_r = 3时，上下游圆柱的边界层厚度明显不同，在前滞止点附近处，上游圆柱边界层厚度明显薄于下游圆柱边界层厚度，而热边界层较薄，热阻较小，使得局部努塞尔数较大，所以，上游圆柱壁面局部努塞尔数大于下游圆柱壁面局部努塞尔数，这与图9中给出的上下游圆柱在U_r = 3的局部努塞尔数沿圆柱壁面分布图相对应。从图10(b)中可以看出，当U_r = 5时，上下游圆柱在前滞止点附近处的热边界层厚度都较薄，而前滞止点处的热边界层较薄，热阻较小，使得局部努塞尔数较大。但是，在前滞止点之后，边界层逐渐变厚，热阻随之增大，局部努塞尔数逐渐减小，明显可以看出在下游圆柱的尾迹处的边界层逐渐变厚，并且要大于上游圆柱尾迹处的边界层的厚度，这就使得上游圆柱的时间平均努塞尔数大于下游圆柱的时间平均努塞尔数，此结论与图7中时间平均努塞尔数随U_r的变化情况所展现出的结论完全一致。由于上下游圆柱的振动，驻点的位置也发生了变化。剪切层的驻点和分离点都改变了热边界层的发展，影响了上下游圆柱与流体之间的传热。流动剪切层的分离、涡脱落、结构振动和热边界层的演化相互作用，实现了流动、结构振动和传热的耦合。

Figure 10. Local magnification of temperature isolines of upstream and downstream cylinders at ζ = 0: (a) U_r = 3, (b) U_r = 5

图10. ζ = 0时上下游圆柱温度等值线的局部放大图：(a) U_r = 3，(b) U_r = 5

4. 结论

本文采用数值模拟方法，在低雷诺数Re = 150时，串列双圆柱的间隙比为8，质量比m^* = 3，折合速度在1 ≤ U_r ≤ 9范围内，阻尼比ζ = 0, 0.01, 0.1下，对串列双圆柱二自由度涡激振动时的流动特性和对流传热进行了数值研究。分析了振荡系统阻尼比ζ和U_r对流动与传热以及圆柱振动的影响，主要研究结果如下：

1) 上游圆柱在不同阻尼比ζ下的横向振幅值和流向振幅值随U_r的增大先增大后减小，并且上游圆柱在三种不同阻尼比ζ下均在U_r = 3时取得横向振幅最大值和流向振幅值。下游圆柱在阻尼比ζ = 0和阻尼比ζ = 0.01的曲线均在U_r = 5时取得横向振幅最大值和流向振幅最大值，而阻尼比ζ = 0.1的曲线在U_r = 4时取得横向振幅最大值和流向振幅最大值。

2) 上游圆柱在阻尼比ζ = 0时，除了折合速度U_r = 6时，上游圆柱的运动轨迹呈现不规则变化，在其余折合速度下，上游圆柱的运动轨迹呈现“8”形轨迹。当折合速度U_r > 5后，下游圆柱的运动轨迹受上游圆柱尾迹区形成的波动的影响，下游圆柱的运动轨迹呈现不规则变化。

3) 对于上游圆柱二自由度振动时，在较小的折合速度2 < U_r ≤ 5范围内，不同阻尼比对圆柱传热的影响比较显著，表现在U_r = 3时最为明显。对于下游圆柱二自由度振动时，由于上游圆柱二自由度振动的影响，折合速度在3 ≤ U_r ≤ 8范围内，不同阻尼比对圆柱传热的影响比较显著，也扩大了折合速度U_r对圆柱传热影响的范围。并且在U_r = 5时，阻尼比为ζ = 0时的时间平均努塞尔数比阻尼比为ζ = 0.1时的时间平均努塞尔数大8.74%。