基于时间步长自适应技术的海上浮基风电平台绕流数值模拟分析

王开加，程建生，段金辉，王景全，王 亮

（解放军理工大学野战工程学院，南京 210007）

0 引 言

随着能源问题的日趋严峻，开发海上风能并向远海、深海发展是必然趋势，漂浮式结构是深海基础的主要方案。而解决圆柱绕流问题是分析深远海漂浮式风电基础平台在波浪和海流的作用下的动力特性的关键。对于多圆柱的研究一直是流体力学领域中一个重要研究课题。由于绕流场的位置和形状的复杂性，数值求解成为十分重要的手段[1]。

目前，对圆柱绕流采用了不同的方法。对不同雷诺数情况下的圆柱绕流问题包括二维数值模拟作了研究。而在二维数值模拟研究中，易获得较好的阻力系数DC 和斯特鲁哈数St，但升力系数LC的数值计算结果与真实值相差较大，这是由于流场网格（特别是边界层附近区域的流场网格）不够细密及时间步长的选取不合理造成的。在流固耦合的数值模拟中，网格设计及时间步长是影响数值模拟精度的重要因素。本文在应用自适应时间步长理论及小雷诺数（Re=100）情况下，采用有限体积法对圆柱绕流进行了数值模拟，并对流场特性作了分析。

1 数学模型

1.1 概念化设计

以REpower公司研制开发的5MW的水平轴风力发电机为设计依据，自主设计了一种风电浮式基础平台见图1，结构主体采用半潜式浮式平台，为了降低整个系统的重心，底部采用悬挂重物的办法。概念设计的风电浮基系统由3个部分组成：风电机组和塔筒组成的风机系统，甲板、立柱和双浮体组成的半潜式平台和包括6根连接杆件在内的压载物部分。确定了初步设计尺寸参数，并建立了几何模型，见图2。

图1 海上浮基风电平台概念化设计

图2 海上浮基风电平台模型

1.2 参数的定义

在计算流体力学中，一般用无量纲的升力系数LC和阻力系数DC 来衡量某物体的受力情况。升力系数和阻力系数定义如下：

式中：DF和LF——分别为作用于单位长度圆柱上的阻力和升力；U——入流速度。

斯特鲁哈数St是描述圆柱绕流的一个非常重要的无量纲量，表征旋涡脱离情况。其定义为：

式中：T——圆柱单侧旋涡脱落周期；D——圆柱直径；U——未受干扰的自由来流速度。

1.3 控制方程

对于不可压缩黏性流体，在笛卡尔坐标系下，其运动的数值解受控于N-S方程，平面流动的连续性方程和动量方程为：

式中：u、v——分别为x，y方向的速度分量，p——压力，υ——流体的运动黏性系数，ρ——流体密度。

2 数值计算

2.1 边界条件

流动方向与垂直流动方向上的均匀速度分别定义为u和v。

进口处：给定速度进口，u=U,v=0；

出口处：在计算域的出口处，设置压力出口边界条件；

上下边界：采用对称边界条件；

圆柱表面：给定无滑移、无穿透边界条件，u=U,v=0。

2.2 离散格式及求解

计算雷诺数（Re=100）为小雷诺数，属层流范围，故采用Laminar模型。采用基于压力基的分离式求解器进行求解，计算中采用具有二阶隐式时间格式的非定常流动进行计算。压力项与速度项的耦合项计算采用PISO算法实现，压力项离散采用具有二阶精度的格式离散，动量方程采用二阶迎风格式离散。计算中压力、密度、体积力和动量项的欠松弛因子分别为0.3、1、1和0.7。

2.3 几何模型

结合本文研究目标，取圆柱直径 2mD= ，计算区域为5030DD×的矩形区域，如图3所示。上游尺寸20D，下游尺寸30D，D为物体垂直于来流方向平面上的特征尺寸，对圆柱一般取直径。

图3 计算模型

图4 网格划分

3 网格和计算时间步长

3.1 网格划分

3.2 自适应时间步长

在流固耦合的数值模拟中，时间步长是影响数值模拟精度的重要因素。如果时间步长过大，则圆柱和流场的耦合计算不够紧密，会导致计算结果精度不高，甚至会出现错误的结果。如果时间步长太小，则无谓地增加计算量。目前有关文献资料中对时间步长的选择方案有两种：一种是取圆柱自振周期的0.2%[2]，另外一种是时间步长取远处来流速度经过圆柱直径长度所耗费时间(对流时间单位)的0.2%[3]。

在Fluent中，当求解器采用PISO时，即使时间步长取得较大，计算仍然会收敛，但结果的准确性无法得到保证。关于时间步长的选取，引入CFL稳定性限制条件，即

式中：cV——网格单元的体积（二维时取为面积）；d——维数（二维计算取为2，三维计算取为3）；——相应网格单元内流速。式(6)表示计算的时间步长tΔ近似等于流速从一个网格点传播到另一个网格点所需的时间。由于采用的是非均匀网格以及每个网格点处当地的流速不一样，所以每个网格点处的tΔ是不一样的。因此，在给定的时刻t和给定的网格点i上，式(6)应写成：

为了保证流场在推进求解中不发生“扭曲”现象，选择在所有网格点（总数设为M）上计算(i=1,2,…,M)，从中选择最小的一个，再附一调节系数α，即取：

式中：α——自适应时间步长适应系数。

采用上述的自适应时间步长方案，通过用户自定义功能UDF的二次开发接口对Fluent软件二次开发来实现自适应时间步长方案在圆柱周围流场的数值模拟中应用。

4 计算结果及分析

4.1 模拟工况

采用同一条件进行绕流计算：入流速度为1m/s，圆柱直径D=2m，Re=100。对自适应时间步长适应系数α为0.5、1.0、1.5和固定时间步长为0.025、0.050、0.100共6种工况进行了数值模拟。数值模拟结果见图5～7，同其他数值模拟结果和实验结果的比较见表1、2。

4.2 结果分析

4.2.1 时间步长的影响

针对上述6种工况进行数值模拟，给出了采用自适应时间步长时α对计算结果的影响，如表1所示，通过与前人研究结果比较知，该计算结果是可靠的并提高了计算精度。从表1中可以看到，随着α的增大，St数、升、阻力系数都会逐渐减小。当α=1.0时，St数、升、阻力系数相对于α=1.5时的数据变化不大，但实际运算中计算量却要大50%，因此全面考虑计算时间、效率以及经验，取α=1.5时自适应时间步长算法为宜。采用固定时间步长时，St数、升、阻力系数随时间步长dt增大而减小的趋势比较明显。如当dt= 0 .100，计算所得的 St数和阻力系数相比于 dt= 0 .025时减小得并不太多，但升力系数的振幅却下降了约三分之一。

图5是Re=100，1.5α=情况下的时间步长时程曲线。从图中可以明显地看出，1.5α=情况下采用自适应时间步长算法的时间步长均大于0.025，并且计算结果优于后者。

表1 升阻力系数及斯特鲁哈数对比表

图5 Re=100，α=1.5情况下的时间步长时程曲线

图6给出了采用自适应时间步长(1.5α=)时圆柱升力系数和阻力系数的时程曲线。由于圆柱尾流形成了卡门涡街，当周期性的涡脱处于圆柱上方时圆柱的升力最大，而后涡脱会逐渐向下游运动，经过圆柱中心时升力系数为0，当涡脱运动到圆柱下游时，升力系数达到最小负值，如此往复，故升力系数的均值为0。阻力系数在涡脱达到稳定后的均值为1.342，这与文献[4]中Re为100的计算结果非常吻合。

由图6阻力系数与升力系数时程曲线可以看出，升力系数发生1个周期变化的同时，阻力系数就会发生2个周期的变化，即阻力系数的频率约为升力系数频率的2倍。这是由于圆柱发生周期性涡脱时，从上下表面脱离的涡会引起阻力改变1次，而这2个涡共同影响升力变化1次。

图6 Re=100，α=1.5情况下的圆柱阻力系数与升力系数时程曲线

图7给出了圆柱边界层附近的矢量场，从中可以看到在边界层内，矢量沿圆柱外法线方向逐渐增大，准确地模拟出速度矢量在边界层内的变化规律。

图7 圆柱边界层附近矢量场

4.2.2 网格设计的影响

在静止圆柱绕流的数值模拟中，比较容易获得较好的阻力系数DC和斯特哈尔数tS，升力系数LC的数值计算结果有时同真实值相比则偏低。原因是因流场网格划分不够精细，特别是边界层网格划分不够精细造成的。在边界层网格的设计中，除了沿圆周有足够的网格数，而且第一层网格到圆柱壁面的距离也要足够小。表2给出了一个比较：

表2 粗糙网格和精细网格对数值模拟结果的影响

从表2的比较可以看出，采用精细网格比采用粗糙网格能获得更好的数值模拟结果，特别是升力系数的幅值有了提高，更加接近精确值。因此要想获得较好的数值模拟结果，流场计算区域要足够大，流场网格划分要足够精细，边界层附近区域的流场网格除了沿圆周要有足够多的网格数，而且第一层网格到圆柱壁面的距离要足够小。

5 结 语

深海风电具有风资源丰富、无噪音污染、风电设备利用率高、不受水深和海底地质条件限制等优点，是未来海上风电发展的主要方向。风电浮式基础，是一种专门为深海区域风电开发的新型海洋基础，但由于深海区域的各种复杂环境条件引起的作用力对基础结构的安全造成很大的影响，目前国内外对于深海风电浮基的研究尚处于起步阶段。采用时间步长自适应技术和精细的边界层网格处理技术能够获得较为精确的数值结果。

在雷诺数较低、圆柱周围的流动主要呈现二维流动的情况下，采用基于有限体积法的自适应时间步长算法，流场网格划分较好、特别是边界层网格划分较精细，完全可以用Fluent软件准确地模拟静止圆柱绕流问题。本文应用Fluent软件对单圆柱绕流流场进行数值模拟，所采用的自适应时间步长算法在现有的计算设备下可以获得较好的数值模拟结果，与固定时间步长算法相比提高了数值计算的精度。为后续风电浮式基础结构的优化设计和将来深海风电浮基安装、现场作业等具有指导性的意义。

[1] 王 健，李海涛. 计算流体力学方法在船舶领域的实用性研究[J]. 船舶与海洋工程，2012, (4): 6-11.

[2] 曹丰产，项海帆. 圆柱非定常绕流及涡致振动的数值计算[J]. 水动力学研究与进展（A辑），2001, 16(1): 111-118.

[3] Newman D J, Karniadakis G E. A direct numerical simulation study of flow past a freely vibrating cable[J]. Journal of Fluid Mechanics, 1997, 344: 95-136.

[4] Kim J, Kim D, Choi H. An immersed-boundary finite-volume method for simulations of flow in complex geometries [J]. J Comput Phys, 2001,171:132-150.

[5] Braza M, Chassaing P, Ha Mind H. Numerical study and physical analysis of the pressure and velocity fields in the near wake of a circular cylinder [J]. J Fluid Mech, 1986, 165: 79-130.

[6] Lu XY, Dalton C, Yhang JF. Application of large eddy simulation to an oscillating flow past a circular cylinder [J]. Journal of Fluids Engineering, Transactions of the ASME. 1997, 119(3): 519-525.

[7] 苏铭德，康钦军. 亚临界雷诺数下圆柱绕流的大涡模拟[J]. 力学学报，1999, 31(1): 100-105.