水下压缩空气储能系统储气装置的CFD数值模拟

王金舜，王 虎，熊 伟，王志文

(大连海事大学 船舶机电研究所，辽宁 大连 116026)

引言

随着当前能源结构的变革，可再生能源正发挥着前所未有的重要作用。但可再生能源固有的“间歇性”、“随机性”和“低能量密度”的特点极大地阻碍了其高效可靠的利用[1-2]。这些由可再生能源产生的波动的电能和波动的电力负荷之间存在着严重不平衡。压缩储能技术是解决这种不平衡问题最有效的方法之一[3]，水下压缩空气储能技术是一种适合在沿海地区、海岛、海洋平台等区域建设的，规模化的可再生能源存储技术。与陆上压缩空气储能系统不同之处在于，水下压缩空气储能系统中的储气装置需要放置在水下环境中，可以保证即使发生失效事件也不会对整个系统或者周围设施与居民造成严重的危害[4-5]，这虽然在实现等压存储的基础上提高了系统的安全性，但复杂多变的海洋环境也对水下储气装置的设计提出了不小的挑战。因此，研究储气装置周围的流场分别在普通与极端恶劣水下环境时的变化特征是十分有必要的。

CFD数值模拟方法被广泛应用于海洋工程结构物的设计过程当中，如水下立管、海上平台桩腿等。在水下压缩空气储能技术领域，数值模拟同样是分析水下储气装置周围流场的有效研究方法之一。

VASEL-BE-HAGH等[6-9]使用了标准k-ω和LES Dyna-SM两种湍流模型分析了在雷诺数为2.3×105的流动环境中气球状储气装置周围的流动结构及装置的受力特性，阐明了装置后旋涡可能的发生形式及脱落过程，并得到了该装置的升阻力系数及涡脱频率。通过对比两种湍流模型的分析结果表明，使用标准k-ω模型对力系数的计算结果较大涡模拟结果偏大。WANG等[10-12]使用标准k-ω模型研究了在9种不同雷诺诺数下的气球状储气装置的受力特性，并将数值模拟和实验的力特性与流动结构结果进行了比较和讨论。此外，他还使用k-ωSST湍流模型对3个不同自由端形状的等效特征长度的钝体进行了数值模拟，并对比分析以研究自由端的形状效应，结果表明，k-ωSST湍流模型可以正确地预测尾流中气球状钝体和流动结构的时均力特性，但不能有效捕捉到瞬态流场特性。继续使用大涡模拟对气球状储气装置流体动力学特性进行了数值模拟，结果表明大涡模拟既能够准确预测流场时均特性，也能够有效捕捉到流场瞬态特性。此外，由于自由端和锥形效应，还在钝体后产生了独特的涡旋结构。

本研究采用大涡模拟方法，模拟在不同海流速度条件下的刚性水下储气装置的绕流动态过程。首先通过分析流场域内速度的分布特性，以了解储气装置周围不同位置处的流场信息以及该信息关于初始流速的敏感程度。然后以相同时间步长为间隔，分析一段时间内不同时刻的旋涡的发生及脱落的动态过程。最后与装置的受力变化曲线相互验证，分析旋涡脱落对装置受力的影响并得到装置所受流体力的主频。

1 数学模型

1.1 控制方程

本研究采用大涡模拟(LES)的湍流模型对水下储气装置的绕流流场进行数值模拟。大涡模拟通过滤波函数将旋涡在空间过滤，对尺度较大的旋涡进行直接模拟，再利用亚格尺度模型模拟小尺寸旋涡对大尺度旋涡的影响[13]。不可压缩流体的N-S方程为：

(1)

式中，ρ—— 流体密度

p—— 压力

ν—— 运动黏度

ui,uj—— 速度分量

xi,xj—— 位移分量

t—— 时间

经滤波函数处理N-S方程，可以得到大涡模拟的控制方程为：

(2)

(3)

经过滤所带来的附加未知项τij为SGS应力，用来描述小尺度涡对大涡的影响。通常使用由SMAGORINSKY[14]提出的SM模型来计算SGS应力，如式(4)所示：

(4)

δij—— Kroneker符号

CS—— 常数，一般取0.1

Δ—— 网格滤波尺度

1.2 装置的几何模型与计算网格

本研究中所使用的储气装置模型如图1a、图1b所示，模型经简化处理为直径D=10 m、总高度H=15 m的球顶圆柱体结构，可在内部储存大约1000 m3的压缩空气以满足设备储存容量要求，如图1c所示,装置外径10 m，同时也是计算不同工况下的雷诺数Re时所使用的特征长度。

图1 储气装置几何模型与计算模型

在中国南海海域内，根据当地水文条件及洋流活动资料显示，该海域内海底水流速度全年在0.5～2.0 kn 范围内，冬季极端条件下海流速度可达2.8 kn以上。此外，水下温度全年在8 ℃以下[15]。因此根据条件，本研究内设置4种工况以对应不同流速的海洋环境，各个工况由代号Un表示(代号中n对应当前工况下的海流速度)。不同工况的设置信息及雷诺数计算由表1给出。

表1 工况设置信息表

储气装置计算域模型如图2所示，采用笛卡尔坐标系为参考坐标，X，Y，Z分别代表顺流方向、高度方向和横流方向，坐标原点在速度入口面底部中心处，3个方向上的计算与尺寸分别为X∈(0，20D)，Y∈(0，3D)，Z∈(-5D，5D)。模型底面中心坐标为(5D，0，0)，到出口边界的距离为15D以降低出口边界对尾迹流场的影响，同时保证能够监测到更大范围的绕流尾迹。经计算，阻塞比为4.64%。根据不同工况设置对应流速的速度入口，流域底面及装置外表面均为固定壁面，其余上、侧表面为滑移壁面以降低壁面边界对内部流场的影响，壁面移动速度与相应工况下的海流速度一致，流域出口边界条件为自由流出口，流量权重为1。

图2 计算域模型

图3所示为装置计算域网格划分，在划分网格时，将装置周围及后方流域做加密处理，并在装置壁面以及流域底面划分边界层网格，用以保证装置壁面受力的准确监测。图3a～图3c分别为网格的+X，+Y，+Z视图，图3d、图3e为装置边界层网格。

图3 计算域网格设置

2 结果与分析

2.1 流动结构

通过监测流场的时均速度沿各个方向上的分布情况，可以较为准确地划分流场内不同流动特征的特殊流场区域，从而完整地描述整体流场，再根据所得信息，对特殊流场区域进行重点分析。

经过无量纲化处理后，图4a、图4b分别表示在装置半高处，流向方向时均速度分量u与竖直方向时均速度分量v沿流向的分布情况。

图4表明，装置在不同流速的海洋环境中所形成的尾迹流场特征基本一致，但在形成位置以及速度量值上有所区别。图4a所示的流向方向速度分量u在到达装置迎风面之前迅速下降至0，在绕过装置后，在距离背风面较近的尾流区内出现负速度，表示在这个流场范围内，绕过装置的流体的流动方向与来流速度相反，朝着装置背风面流动，该区域称之为回流区，随后，向流场远处发展的流体速度逐渐增加至正值并向来流速度接近直至达到稳定，并且通过4种工况的速度曲线对比，可以看出在0.5 m/s流速下，海流绕过装置后形成的回流区较其他3种工况范围更大，并且达到稳定的位置有所滞后。

图4 不同流速下时均速度沿流速方向分布示意图

不同工况下的u特征信息在表2中给出。从表中数据可以得出，4种工况下的最大回流速度umax均在0.15 Un附近，回流区的长度随速度的增加而减小，不论是在何种工况下，在装置后方均存在回流现象且逐渐稳定。0.5 m/s流速下达到稳定的位置在更远的10D处。

表2 不同工况下的u特征信息表

图4c表示竖直方向的时均速度分量v沿流向的分布，从图中可以看出，与u的分布情况类似，装置在不同流速的海洋环境中该方向速度分量分布情况基本一致，但在特征出现位置以及速度量值上有所区别。海流在到达装置迎风面前方时有沿竖直方向向上运动的趋势，流过装置侧面后，在背风面附近v迅速降低至负值，这表明在此区域内(影响区)，大量的海水向装置底部运动(具体成因将在2.2节中说明)，随着与装置距离的增加，v逐渐从负值上升并趋于稳定。不同工况下的v特征信息如表3所示。从表3中可以看出，虽然4种工况下的v分布曲线形态相似，均在距离装置12倍直径的位置附近达到稳定，稳定值近似于0，这表明在此区域外的流场变化与装置之间的关系很小，但在0.5 m/s工况下，海水绕流后向下运动的影响范围较长且强度较弱。表4列出了在不同工况下、位于装置的不同高度处，向下运动的海流的最大合速度信息，从数据中可以看出，位于装置半高处附近，各个工况下的该速度达到最大值，表示在此区域内，下行海流强度最大。

表4 不同工况下、不同高度处的v特征信息表

根据图4b与表3的数据，在v影响区内，距离装置中心1D的位置，于流场中面沿装置竖直方向再次监测时均速度分量v，虽然在流场区域内的v峰值不在此处，但也可以以此判断向下运动的海流在竖直方向上的分布情况，如图5a所示。在影响区外，距离装置中心2D的位置，于装置半高处沿横流方向监测横流方向时均速度分量w的分布情况，如图5b所示。

表3 不同工况下的v特征信息表

由图5a可知，v的分布沿装置高度呈现抛物线的分布情况，在U1.0，U3.0，U5.0工况下，绕流后海水向下运动的速度峰值出现在装置中高处偏下的位置，并显示出与来流速度之间的良好的线性关系。U0.5工况下，该速度峰值出现的位置在中高处偏上，这是由于在低速流场下，绕过装置顶端的下行流体在来流方向向上位置偏后，在此截面上(X=1D，YOZ截面)流动强度较弱、速度分量v较小所造成的。

图5b显示出流方向时均速度分量w在横流方向上的分布关于流场中面是反向对称的，表明在装置后方有反向旋转的旋涡产生，由旋涡产生的切向速度约为来流速度的0.3倍，在数值上，w的大小受来流速度影响较小。

图5 不同流速下时均速度沿其他方向分布示意图

2.2 旋涡结构

根据对各个工况下的时均速度场分析，可以判断海水流过装置时的整体过程及尾迹流场的时均形态，以此推测出在装置后方，存在具有一定规则的旋涡结构，形成回流、海水下行、与横流速度反向对称现象等局部流场特征。

以0.5 m/s流速的工况为例，图6是在装置半高处、背风面附近的时均流线图。如图所示，装置后方的回流区本质上是一对反向旋涡作用的结果，这是由于流体绕过曲面时速度先升后减，压力先减后升，在曲面后方某一位置产生的边界层分离现象所导致的，流体质点在流场中面附近与来流反向运动，因而速度为负值，形成了回流区。

图6 Y=0.5 L，XOZ截面时均流线图

图7a和图7b分别表示在流场中面及X=1D，YOZ截面的时均流线图。如图所7示，由于自由端的结构是球形，在流场中面内以及与此截面平行的各个截面内自由端的投影为圆形，所以流体绕过自由端时同样会发生边界层分离现象，形成一个顺时针方向旋转的旋涡结构，在背风面附近速度向下。与此同时，在与YOZ截面平行的截面内，如图7b所示，同样存在一对反向对称的旋涡结构，在中面附近速度向下，因此，竖直方向的时均速度分量v的负值是图7a、图7b中旋涡结构共同作用的结果，形成了方向向下的海流。

图7 其他截面时均流线图

图8a和图8b分别表示在X=1D，X=2D的YOZ截面内，Q-criterion显示的涡量与速度矢量合成图的时间快照序列。如图8a所示，随着以1/4T为间隔的时间变化，在该截面虽然能够清楚地观察到旋涡结构，但由于旋涡结构复杂难以直观地发现其发展规律。但是旋涡结构出现的位置相对固定，分别在装置半高处流场中面两侧，以及装置底部流场中面两侧。

在图8b中，距离装置稍远处，在该截面能够清楚地观察到旋涡结构，如果以t=t0为特征流场，包括向左发展的下行海流以及在下行海流两侧存在一对反向旋涡，则当x=5/4T时刻，该特征流场再次出现。在2000～3600 s的时间内，特征流场出现的次数及出现的时间点如表5所示，将相邻2个时刻做差得到的周期分布图如图9所示。经过计算，该周期分布的平均值为95.33 s，即为流场旋涡脱落的主频。图8a所示结果没有规律的原因是距离装置较近，截面处在回流区及下行海流影响区内，该区域是旋涡形成及发展的区域，流场杂乱无章。

表5 特征流场出现时间表

图8 Q-criterion与速度矢量合成图的时间快照序列

图9 特征流场出现的周期分布图

2.3 受力特征

图10a～图10d分别展示了4种工况下的装置受力特征，阻力与升力通过式(5)、式(6)换算为相应的力系数。

图10 不同工况下的装置力系数变化曲线图

(5)

(6)

随着流动的发展，虽然在微观上仍有波动，但装置受到的阻力与升力逐渐趋于稳定。通过监测入口流量与出口流量的平衡与残差收敛保证计算结果的收敛性。经过统计可以得出各个工况下阻力系数与升力系数的均值，见表6。当雷诺数超过3×106时，在自模区范围内力系数不再随Re变化，因此力系数应当不随工况的变化而变化，对各个工况下的力系数进行对比的意义是以这样一种方式进行“重复性实验”，从而再次验证计算结果的准确性，结果表明，储气装置的阻力系数约为0.45，升力系数约为0.60。

表6 各个工况下阻力系数与升力系数的均值表

由于涡街的存在，装置在横流方向上所受的流体力随着时间而周期性变化，同样称之为升力，用lift-z表示，以U0.5工况的计算结果为例，图11是装置所受横流方向升力的升力系数时程图，升力系数在0值附近正负波动证明了装置侧向存在一对交替产生的旋涡。经傅里叶变换后，换算成斯特罗哈尔数-幅值的关系曲线如图12所示。

图11 U0.5工况下横流方向升力系数时程图

图12 U0.5工况下斯特罗哈尔数-幅值曲线图

斯特罗哈尔数是表征流动周期性的相似准则，对于周期性的非定常流动，斯特罗哈尔数可以表示为速度与特征频率之间的关系，从而确定旋涡脱落的频率。图12中图像的最大幅值对应的频率就是旋涡脱落的主频，斯特罗哈尔数计算结果为0.17，表明在此工况下，涡脱的周期约为100 s，与特征流场出现的周期基本一致(见小节2.2)。

3 结论

使用LES模型对装置处于不同流速的海洋环境中所形成的尾迹流场及受力特征进行了数值模拟，通过对流场内各个方向上的时均速度分量(u,v,w)、速度矢量、涡量，以及流体力等物理量的监测，对海流通过装置后的流场及装置的受力情况进行了研究，主要研究结论如下：

(1)不同速度的海流通过装置后形成的流场特征相似，特征速度量值与特殊流域的范围与海流速度有关，低流速海洋环境下的绕流影响区较广，各个特殊流域的位置均有滞后；

(2)绕流形成的旋涡具有明显的三维结构特征，经储气装置侧向流过的海流在装置背风面附近逐渐向底部扭曲形成旋涡结构，背风面附近区域的两侧旋涡关于流场中面对称分布；

(3)绕过装置球顶后形成的分离涡导致海水向下运动，在装置半高处附近与侧向向下扭曲的旋涡合流，所形成的下行海流是尾迹流场中竖直方向负速度流体出现的主要原因之一；

(4)装置所受阻力与升力随流动的发展趋于稳定，但这种稳定是一系列无规则波动的力叠加而成。通过统计得出装置的阻力与升力系数分别为0.45和0.6，以此可以计算在任意流速海洋环境下装置所受到的流体力；

(5)侧向形成的三维旋涡在距离中心2倍直径的位置附近开始发生脱落，以横流方向升力系数与流场速度矢量、涡量的变化规律为依据，对旋涡脱落周期定性定量分析，结果表明该旋涡脱落的周期在100 s左右，会对装置造成低频扰动。