全风向角下小间距双方柱绕流的大涡模拟

摘要

采用大涡模拟方法研究了间距比（P/B）为1.5（P为柱心间距，B为方柱边长）、风向角α 为 0°~90°等条件下双方柱在均匀来流作用下的气动力、流态划分、表面风压和流场特性。研究发现：小间距双方柱流动干扰效应显著，下游柱平均气动力随风向角的变化规律与上游柱和单方柱差异较大，且可能受更大的升力绝对值。将小间距比双方柱绕流分为前角分离流态（α = 0°~10°）、分离泡流态（α = 20°~30°）、附着流流态（α = 40°~60°）及间隙侧分离泡流态（α = 70°~90°）四种模式。附着流流态和间隙侧分离泡流态的间隙区出现较强的负压。尾流负压区的强度随风向角先增强后减弱，在α = 50°附近达到峰值。前角分离流态时双柱具有一个整体的尾流回流区，而在其他流态下上下游方柱均有独立的回流区。

关键词

双方柱; 小间距比; 大涡模拟; 气动干扰; 流场机理

超高层建筑群气动干扰效应显著，其表面风荷载与单体建筑差异较大^［1-2］。影响超高层建筑群风荷载的因素众多^［1］，如来流风特性、建筑外形、建筑群位置关系、建筑间距、风向角等，情况非常复杂，其气动干扰的流场机理尚不明确。为了更好地理解超高层建筑群的气动干扰效应及其流场机理，可将大气边界层下超高层建筑群周围流动干扰问题简化为均匀来流下双方柱绕流问题^［1-3］，研究双方柱的间距和风向角等参数对方柱气动力、表面风压和流场结构的影响规律。

针对小间距比双方柱绕流，以往文献主要围绕串列、并列和错列三种布置形式。Alam等^［3］的试验研究和杜晓庆等^［4］的数值研究发现，串列双方柱在间距比P/B = 1.5~1.9时（P为柱心间距，B为方柱边长）下游方柱会出现负阻力。Agrawal等^［5］研究发现小间距比并列双方柱会出现偏向流现象，即两个柱体的尾流区域具有显著的宽窄差异；而随着间距比增加，偏向流逐渐消失。陈素琴等^［6］进一步指出偏向流情况下尾流区域较窄的柱体其升阻力和涡脱落频率较大。针对错列双方柱绕流，Sakamoto和Haniu^［7］的试验研究发现当两方柱间距较小且下游方柱处于上游方柱尾流中时，下游方柱的脉动阻力较小。Aboueian和Sohankar ^［8］则对低雷诺数（Re 为 150）下，45°固定错列角的变间距双方柱进行了流态划分。陈素琴等^［9］的数值研究发现P/B = 2双方柱在错列角大于60°时会出现偏向流现象。Alam等^［10］指出P/B=1.5，错列角为80°的双方柱尾流会发生撞击‒反弹现象，使得下游方柱形成窄尾流。总的来说，小间距比双方柱绕流的流场干扰效应较为显著。

风向角对小间距比双方柱绕流的影响仍有待研究。已有较多文献研究了风向角对单方柱绕流的影响^［11-14］。变风向角下双方柱绕流的研究目前仍较少，Du等^［11］基于风洞试验，系统研究了在雷诺数Re = 8×10⁴下，风向角α = 0°~90°，间距比P/B = 1.25~5条件下双方柱的气动性能，并将双方柱绕流问题划分为小间距（P/B < 1.5）、中等间距（1.5 ≤ P/B ≤ 3）和大间距（3 < P/B ≤ 5）三种情况；发现小间距比时双柱气动性能随风向角的变化较为复杂，体现出较强的流动干扰效应，但未能从流场结构角度分析其气动干扰的机理。

本文以雷诺数Re = 8×10⁴（Re = ρUB/μ。式中，U、ρ、μ分别为来流速度、流体的密度和黏性系数）、间距比P/B = 1.5、风向角α = 0°~90°的双方柱绕流为对象（图1），采用大涡模拟方法研究了风向角对双方柱的气动性能和绕流场特性的影响规律，重点从流场角度阐明了不同风向角下双方柱流动干扰的机理。

图1  计算模型和边界条件示意图

Fig. 1  Sketch of computational model and boundary conditions

1 数值模型与验证

1.1　数值方法

大涡模拟法（large-eddy simulation，LES）是一种在网格尺度上运用空间平均的方法除去高频成分的物理量计算，此方法可实现湍流流场的再现。大尺度涡通过滤波后的不可压缩Navier-Stokes方程直接求解：

$$\frac{\partial {\overline{u}}_i}{\partial x_i}=\mathrm{0}$$

（1）

$$\frac{\partial {\overline{u}}_i}{\partial t}+\frac{\partial {\overline{u}}_i{\overline{u}}_j}{\partial x_j}=-\frac{\mathrm{1}}{\rho }\frac{\partial \overline{p}}{\partial x_i}+\nu \frac{{\partial }^{\mathrm{2}}{\overline{u}}_i}{\partial x_j\partial x_j}-\frac{\partial {\tau }_{ij}}{\partial x_j}$$

（2）

式中：${\overline{u}}_i$为滤波后速度分量，x_i为空间位置分量，t为时间；$\overline{p}$为滤波后压力，v为流体运动黏度。τ_ij为亚格子尺度应力，它体现了小尺度涡的运动对所求解的运动方程的影响。综合考虑计算效率和精度，本文采用了壁面自适应局部涡粘（wall-adapted local eddy-viscosity，WALE）亚格子尺度模型模拟τ_ij，WALE模型常数为0.325。计算中采用SIMPLEC算法求解压力速度耦合方程，采用中心差分格式进行空间离散，采用二阶隐式的时间离散方案。

1.2　计算模型与边界条件

图1为全风向下双方柱计算模型和边界条件的示意图，其中上游方柱C1与下游方柱C2的相对位置不随风向角α变化。C_L1和C_D1为C1的升力和阻力。C_L2和C_D2为C2的升力和阻力。柱体截面的角点位置分别用a、b、c、d标出。本文取雷诺数Re = 8×10⁴，采用的柱心间距固定为P = 1.5B，风向角α = 0°~90°（每隔10°一个工况）。流场采用O型计算域，其直径为60B；双柱居中布置于计算域中，随着风向角变化阻塞率的变化范围为1.67 % ~ 3.33 %，计算域的展向长度为L = 4B。计算中采用了速度入口边界条件和自由出口边界；在计算域展向采用了周期性边界，方柱表面则采用了无滑移壁面边界。图2展示了结构化计算网格，每个方柱沿周向布置300个网格并在角部适当加密；近壁面最小网格厚度为0.001B，满足$y^{+}=\mathrm{\Delta }{y}_p\sqrt[]{{\tau }_{\mathrm{w}}}/v\sqrt[]{\rho }\approx \mathrm{1}$（Δy_p是壁面节点到壁面的距离，τ_w是壁面剪应力）；展向等间距划分40层网格。双方柱的计算模型网格总单元数为337万。计算中采用固定的量纲一时间步长Δt^* = ΔtU/B=0.005，其中Δt为有量纲时间步长；在流场计算稳定后进行采样并计算统计值，采样时长不少于60个漩涡脱落周期。

图2  网格示意图

Fig. 2  Details of computational mesh

1.3　模型验证

为验证计算方法和计算参数的合理性和准确性，首先以单方柱为对象，研究了周向网格数、展向长度和量纲一时间步长等参数对平均阻力系数C_D、脉动阻力系数C_D'、脉动升力系数C_L'和斯托罗哈数St的影响，并将计算结果与文献值^［1，12，15-17］进行了比较，结果如表1所示。主要气动参数的定义如下：平均阻力系数C_D = 2F_D/ρU²BL，平均升力系数C_L = 2F_L/ρU²BL，平均风压系数C_p = 2p/ρU²，斯托罗哈数St = f_vB/U，其中F_D，F_L和p分别为方柱所受阻力、升力和表面压力的时间平均值，f_v为漩涡脱落频率；而脉动阻力系数C_D'、脉动升力系数C_L' 和脉动风压系数C_p'则分别对应阻力、升力和表面压力的标准差。可见，周向网格精度、时间步精度和展向长度对C_D、C_L'和St的影响相对不显著，各个算例的计算结果均处于文献结果的分布范围内。另一方面，随着周向网格精度、时间步精度和展向长度的提高，C_D'结果同文献的误差明显缩小。总的来说，周向300个网格，时间步长Δt^* = 0.005，展向长度为L = 4B的工况2可以得到同文献结果较为接近的气动力参数，且计算资源的消耗相对较少。

表1  单方柱绕流计算验证

Tab. 1  Verification and validation of flow past a single square cylinder

数据 来源	工况	Re /104	周向 网格	Δt*	L	CD	CD'	CL'	St
本文 （LES）	周向网格数量比较
	工况1	8	200	0.005	4B	2.08	0.25	1.38	0.119
	工况2	8	300	0.005	4B	2.00	0.22	1.36	0.118
	工况3	8	360	0.005	4B	1.99	0.19	1.36	0.120
	量纲一时间步长比较
	工况4	8	300	0.002 5	4B	2.02	0.18	1.32	0.120
	工况2	8	300	0.005	4B	2.00	0.22	1.36	0.118
	工况5	8	300	0.01	4B	2.04	0.21	1.32	0.118
	计算域展向长度比较
	工况6	8	300	0.005	1B	1.98	0.35	1.32	0.118
	工况7	8	300	0.005	2B	2.01	0.30	1.34	0.119
	工况2	8	300	0.005	4B	2.00	0.22	1.36	0.118
试验	杜晓庆等［1］	8			13.5	2.25	0.22	1.35	0.132
	Vickery［15］	10			14	2.05	0.17	1.30	0.120
	Lee［12］	18			9.2	2.04	0.22	1.19	0.120
LES	Rodi［16］	2.2	－	－	4	2.02~2.77	0.14~0.27	1.15~1.79	0.09~0.15
	Sohankar等［17］	2.2	－	－	4	2.03~2.32	0.16~0.20	1.23~1.54	0.126~0.132

图3a、图3b分别给出了Case 2单方柱表面平均风压系数C_p和脉动风压系数C_p'分布（由方柱展向中间截面周向测点的风压时程曲线得到），并与文献^［1，12，18-19］结果进行比较。由图3可知，工况2的平均风压系数与和脉动风压系数均与文献的结果吻合良好，分布趋势一致。综合考虑计算精度和计算效率，本文采用工况 2的计算参数进一步分析双方柱的绕流问题。

图3  单方柱的表面风压系数分布

Fig. 3  Distribution of pressure coefficients on a single square cylinder

2 结果与分析

2.1　气动力系数

图4为双柱平均阻力系数C_D和平均升力系数C_L随风向角的变化曲线，同时给出了杜晓庆等^［1］的试验结果作为对比。整体来看，双方柱的平均阻力系数与文献值的变化趋势一致，吻合良好。上游方柱（C1）的变化趋势与单方柱的基本相似，而下游方柱（C2）的趋势与单方柱相比明显不同，这是由于方柱C1的尾流干扰所致。在图4a中，上游方柱C1在α = 10°和70°时平均阻力出现极小值；方柱C2的平均阻力随风向角的增大而增大，且仅在α = 0°时出现负阻力，其中在α = 0°~10°时的变化更为剧烈，这是由于方柱C2受到方柱C1的干扰更为强烈。在图4b中，方柱C1的平均升力在α = 10°和70°分别出现极大值和极小值；方柱C2的平均升力总体呈减小的趋势，其中在α = 0°~10°时的变化较为剧烈，在α = 70°时达到了极小值。需要注意的时，随着风向角的变化，下游方柱的平均阻力和升力一般小于上游方柱和单方柱，但α = 20°~90°时其升力绝对值可能大于后者。

图4  双方柱的平均气动力系数随风向角的变化

Fig. 4  Variation of mean aerodynamic coefficients of two square cylinders with incidence angle

2.2　流场结构

图5~图8给出了P/B = 1.5，α = 0°~90°串列双方柱周围的平均流线图。依据上游方柱表面剪切层的分离与再附等流场特征，参考单方柱绕流^［20］，并结合双方柱的气动性能，可将不同风向角下小间距比双方柱绕流划分为4种流态，即前角分离流态（α = 0°， 10°）、分离泡流态（α = 20°， 30°）、附着流流态（α = 40°， 50°， 60°）以及间隙侧分离泡流态（α = 70°， 80°， 90°）。

图5  前角分离流态平均流线图

Fig. 5  Mean streamlines in leading-edge separation mode

图6  分离泡流态平均流线图

Fig. 6  Mean streamlines in separation-bubble mode

图7  附着流流态平均流线图

Fig. 7  Mean streamlines in attached-flow mode

图8  间隙侧分离泡流态平均流线图

Fig. 8  Mean streamlines in gap-side separation-bubble mode

（1） α = 0°， 10°为前角分离流态（图5）。方柱C1前角点a，d处发生流动分离，且分离剪切层未在C1发生再附；形成的回流区延伸到方柱C2的尾流，并形成一个整体的尾流回流区，尾流的宽度较小，双方柱的绕流类似单个钝体。其中α = 0°时C1流动分离直接影响C2的c-d和a-b面，而对间隙影响较小。而α = 10°时，C1柱d点的分离剪切层在C2柱d点附近发生了再附并直接影响间隙区域，但C1柱未形成独立的尾流回流区。

（2） α = 20°， 30°为分离泡流态（图6）。流体自C1柱前角点d分离后在C1的c-d面发生再附，形成分离泡；C1柱c-d面的附着流在c点二次分离后进入间隙侧，并最终与C1柱a点分离的流体共同形成次要的尾流回流区；此时柱C1、C2分别形成了次要和主要的尾流回流区，尾流的宽度显著增加。

（3） α = 40°， 50°， 60°为附着流流态（图7），方柱C1的c-d、d-a面没有出现分离泡，流体紧贴其表面流动并在c、a点发生分离；柱C1、C2也分别形成了次要和主要的尾流回流区。

（4） α = 70°， 80°， 90°为间隙侧分离泡流态（图8），方柱C1流动分离的位置改变为c、d点，其中c点的分离剪切层在间隙侧b-c面发生再附，并形成分离泡；b-c面附着流在b点二次分离后，与C1柱下侧的尾流共同组成次要尾流回流区；柱C1、C2分别形成了次要和主要的尾流回流区，其中次要尾流回流区的长度随着风向角增加而增加。其中α = 70°时，C1柱d点分离的流体在d-a面形成分离泡。根据陈素琴等^［6］的研究，小间距比并列双方柱会由于间歇流的偏转而形成宽窄不同的尾流，本文在α = 90°时的C2柱尾流宽度大于C1柱，出现了较为显著的偏向流现象。

2.3　表面压力和流场特性

根据图5~图8的流态划分和各个风向角双柱绕流的特点，挑选了6个具有代表性的工况α = 0°， 10°， 20°， 50°， 70°， 90°作进一步的分析。图9给出了这6个风向角上下游方柱表面平均风压系数C_p的分布。（1）对于前角分离流态（α = 0°和10°），上游柱C1的迎风面d-a为正压，C1的其余面和C2的全部表面整体上都是负压，体现了C1尾流对C2的较强影响；其中α = 0°时双柱间隙回流在C2柱d-a面产生的负压强于尾流侧b-c面，使C2受到负阻力的作用，与Alam等^［3］和杜晓庆等^［4］的试验结果一致；而α = 10°时C2柱c-d面由于不处于C1的尾流中，其压力值明显上升，并在d点附近出现局部正压。（2）对于分离泡流态（α = 20°），一方面C2柱d点附近正压相比α = 10°继续增强；另一方面C1柱c-d面压力明显升高，而d-a面正压的不对称性显著增加。（3）对于附着流流态（α = 50°），C1柱的c-d面已经转变为正压区，原先的正压区d-a面的压力整体上变为负压，但靠近d点仍维持正压；C2柱的c-d面压力进一步增加，并整体上表现为正压。（4）对于间隙侧分离泡流态（α = 70°和90°），上、下游柱的表面压力分布相似，均在迎风侧（C1和C2的c-d面）表现为正压，而其余面为负压；其中α = 70°时方柱C1的b-c与d-a面皆出现较大负压且数值相近，其受到的负升力（见图4）主要是由c-d面与a-b面压差的y方向分量导致的。

图9  典型风向角下双方柱表面平均风压系数分布

Fig. 9  Distribution of mean pressure coefficient on two square cylinders at typical incidences

图10为6个典型风向角下双方柱周围平均风压系数C_p云图。对于较小风向角，即前角分离流态和分离泡流态（α=0°，10°和20°），整个双柱的绕流场基本处于较弱的负压中；C2尾流区负压强度随着风向角增加而增加。对于较大风向角，即附着流流态和间隙侧分离泡流态（α=50°，70°和90°），间隙区域的负压相比较小风向角时明显增强；C1柱后出现了明显的尾流负压区，其负压强度随着风向角增加而增强；C2柱尾流区的负压强度则随着风向角的增加而减弱；α=50°时，C2尾流的负压最强，C1柱的d-a面附近具有较强的负压；α=90°时双方柱的平均风压场呈不对称分布，同图8的平均流线结果一致，表现出偏向流的特点，主要由间歇流的偏转导致^［6］。

图10  典型风向角下双方柱周围平均风压系数云图

Fig. 10  Contours of mean pressure coefficients around two square cylinders at typical incidences

3 结论

本文在均匀来流、雷诺数Re = 8×10⁴、间距比P/B = 1.5、风向角α = 0°~90°的条件下，研究了小间距比双方柱的气动力、流场形态、表面压力和流场特性，分析了干扰效应的流动机理。具体结论如下：

（1） 下游方柱气动力受流动干扰的影响较大，随着风向角增加其平均阻力和升力分别为上升和下降趋势，其数值一般小于上游方柱和单方柱，但α = 20°~90°时其升力绝对值可能大于后者；而上游方柱平均升力和阻力在各个风向角下与单方柱相近。

（2） 根据上游方柱表面流动分离和再附着的情况，可将风向角α = 0°~90°下小间距比双方柱绕流分为前角分离流态（α = 0°~10°）、分离泡流态（α = 20°~30°）、附着流流态（α = 40°~60°）及间隙侧分离泡流态（α = 70°~90°） 4种模式。

（3） 随着风向角的增加，尾流负压区的强度先增强后减弱，在α = 50°附近达到峰值；附着流流态和间隙侧分离泡流态的间隙区负压显著强于前两个流态；在前角分离流态下，双柱具有一个整体的尾流回流区，而在其他流态下上下游方柱均有独立的回流区。

作者贡献声明

杜晓庆：研究选题，研究思路和技术指导，论文撰写。

许庆：数值计算，数据处理与论文撰写。

董浩天：研究选题，技术指导和论文构思，论文修改。

陈丽萍：数值计算，数据处理。