基于LES和DES的定日镜结构风致响应分析

卢春玲， 陈建通， 陈锦焜， 王 强

(1. 桂林理工大学 土木与建筑工程学院， 桂林 541004； 2. 广西建筑新能源与节能重点实验室， 桂林 541004；3. 桂林理工大学 有色金属矿产勘查与资源高效利用协同创新中心， 桂林 541004)

塔式太阳能发电系统需要在开阔场地内大范围安装定日镜，定日镜是其主要设备。当定日镜离中央集热塔较远时，微小的风致振动也会降低定日镜的聚光效率，造成效率损失。此外，在开阔的野外场地上，结构高度低，易受近地面复杂风场影响。并且定日镜受风面积大、刚度低，在风力动力荷载下容易引起定日镜振动、损坏，甚至倒塌。所以，对定日镜的风致响应研究具有重要意义。

针对定日镜的抗风问题，国内外学者进行了大量的研究。Vasquez-Arango等[1]结合实测与有限元分析，对定日镜的模态进行了验证，为定日镜有限元模拟提供了依据。Terrés-Nícoli等[2]通过风洞试验，对定日镜结构在风荷载作用下的应力分布进行了分析研究。Wolmarans等[3]采用单向流固耦合方法，研究了中型定日镜周围涡流脱落过程和瞬态风荷载。Gong等[4]在风洞试验和数值模拟的基础上研究了定日镜的脉动风压特征和风致动态响应。王延忠等[5]采用流固耦合数值模拟方法，研究了定日镜在不同仰角下的风致响应。冯煜等[6]采用AR模型模拟风场时程，对三维风场下定日镜风致动态响应进行了分析。王莺歌[7]将理论分析、风洞试验及数值模拟相结合，对定日镜结构的抗风问题进行了综合性的研究。相较于数值模拟，现场实测和风洞试验研究获取的流场信息相对较少，且数据后处理复杂、试验费用昂贵。对定日镜的数值模拟研究，雷诺平均法(Reynolds-averaged Navier-Stokes, RANS)较为常用，但RANS模型对脉动过程进行时均化处理，这对风致响应分析非常不利。而LES和DES避免了RANS抹平细节的缺陷，可以较准确的模拟脉动风压及风荷载。但其对网格精度及入口湍流的准确度要求较高，利用LES和DES对定日镜结构进行结构风荷载和风致响应的研究还较少。

本文选用LES和DES，结合离散再合成的随机湍流生成法对定日镜在0度风向角下0°、30°、60°镜面仰角(镜面法线与水平夹角)的多种典型工况[8]下进行数值模拟，得到定日镜正反面风荷载时程数据。建立有限元模型，将表面风荷载时程施加在定日镜有限元模型上，对其进行风致响应分析。将计算结果与风洞试验的相应数据进行对比，以验证LES和DES的准确性。

1 风洞试验

定日镜风洞测力试验在湖南大学风洞实验室的HD-2单回流边界层风洞中的低速试验段内进行。该风洞低速试验段尺寸为宽5.5 m，高4.4 m，长15 m。采用有机玻璃制成定日镜刚性模型，缩比为1/10，高1.2 m，宽0.9 m，模拟出镜面、转动轴、支撑臂组件及立柱等定日镜主要组成构件，使模型与实体结构在几何外形上保持相似，实体结构及模型如图1。

(a) 定日镜实体结构

风洞试验根据定日镜实体结构所处地区的地貌特征，采用格栅、尖梯、挡板、粗糙元装置模拟出我国建筑结构荷载规范(GB 50009—2001)规定的B类风场，采用IFA300热膜风速仪对大气边界层模拟风场进行调试和测定。风洞试验参考点选在 1.0 m高度处，对应于实际中 10 m 高度，采用皮托管测控流场的参考风速，参考风速为8 m /s。定日镜模型的基座处固定有六分力天平，从而获得模型六个方向上的测力数据。在定日镜正面和背面各布置144个测点，测点布置如图2，通过64通道DTCnet压力扫描阀、ESP-64HD型传感器及INV306大容量自动数据采集和信号处理系统测量各风压时程。采样频率330 Hz；每个测点记录20 000个数据长度，总采样时间为60.6 s。

(a) 迎风面测点布置

(b) 背风面测点布置图2 风洞试验模型风压测点布置Fig.2 Distributions of pressure measuring points in the model of wind tunnel test

2 数值模拟方法

2.1 计算模型及网格划分

定日镜结构的计算模型为全尺寸，定日镜整体高度为(H)11.93 m，镜面宽度(L)9 m，转动轴和立柱为直径1 m的圆柱，镜面与支架之间留有空隙。计算域z方向高为5H,Y方向宽度为8L,镜面法向(X)方向长度为20L,流域入口距模型镜面5L，其结构阻塞度为2.6%，计算域尺寸如图3。

图3 计算域尺寸Fig.3 Computational domain size

流场网格设置采用混合网格。图4为定日镜表面网格划分方案，图5为流场网格划分方案。定日镜结构处于X方向长3L，Y方向宽2L，Z方向高1.5H的网格加密区中。加密区采用非结构化网格，网格尺寸为0.01～0.5 m，在定日镜结构表面的网格最小，并以5%的增长率逐渐增大。加密区嵌套于非加密区之中，非加密区采用结构化网格，总网格数约为510万。这样既能保证计算效率又能保证定日镜周围的流场计算精度。

图4 定日镜表面网格划分Fig.4 Mesh for heliostat surface

图5 流场整体网格划分Fig.5 Global mesh generation for flow field

2.2 湍流模型

本文采用LES和DES对定日镜结构进行数值模拟计算。对于LES，FLUENT提供了四种亚格子模型，分别是Smagorinsky-Lilly模型、动态Smagorinsky-Lilly模型、壁面自适应局部涡流黏性模型和动态动能亚网格模型。动态动能亚网格模型通过亚网格尺度湍流动能的传输来描述亚网格尺度湍流。弥补了Smagorinsky-Lilly模型在高雷诺数下，能量耗散难符合预期的缺点。因此，采用了动态亚网格动能模型。其应力尺度为

(1)

其中，ksgs通过求解其输运方程(2)得到

(2)

DES中，采用分离涡的一种优化模型——延迟分离涡(delayed detached eddy simulation, DDES)，分离涡模型是Spalart等[9]提出的一种混合数值模型，该方法结合了RANS和LES，通过对网格尺寸相关的长度比例l的判断，自动在两种方法之间转换。如式(3)、(4)所示

l=min(d,CDESΔDES)

(3)

ΔDES=max(δx，δy，δz)

(4)

式中:CDES=0.65，δx、δy、δz表示三个方向的网格比例，d表示网格到壁面的距离。当l=d时，启用RANS计算；当网格从壁面开始增长时，当l=CDESΔDES时，转换为LES进行计算。某些情况下，网格纵横比在近壁面不够大，即使网格厚度很小，LES也会被启用。为了保证计算精度和效率，DDES对d进行修正如式(5)、(6)

d≡d-fdmax(0,d-CDESΔ)

(5)

fd=1-tanh[(8rd)3]

(6)

(7)

式中：rd是湍流尺度与网格到壁面的距离之比，其中vt是动力黏度，v是分子黏度，Uij是速度梯度，k是Von Karman常数。当网格靠近壁面时fd=0，启用RANS计算。当网格远离壁面时rd=1，即fd=1，湍流模型转化为LES。通过对d的修正，DDES延迟了LES的发生，提高模拟的准确度和计算效率。DES亚网格模型沿用动态动能亚网格模型，RANS模型设定为Realizablek-ε两方程模型。

2.3 求解设置

在数值求解中，采用不可压缩流。在LES和DES计算前，导入瞬态化处理的RANS模型计算结果作为初始流场，以加快计算收敛速度。数值计算采用SIMPLE算法，时间离散采用二阶隐式格式，空间离散采用二阶中心格式，计算收敛残差标准取为1×10-4。进行流体计算时，采用Fluent提供的monitor模块对定日镜结构上的每个镜面进行风压监测并对其积分，得到风荷载时程数据。监测面布置如图6。

图6 数值模拟模型监测面布置Fig.6 Pressure monitoring surface on the model for the numerical simulation

模拟计算在高性能工作站上进行，工作站搭载18核36线程的Intel Xeon Gold 6140处理器，运行内存192 GB。LES与DES计算时间步长取0.05 s，每个时间步长最大迭代次数为20次，每个工况进行了10 000个时间步长的计算。LES模拟每个时间步长计算耗时约40 s，每个工况耗时约110 h；DES模拟每个时间步长计算耗时约35 s，每个工况耗时约95 h。

2.4 边界条件

本文LES和DES模拟入口湍流生成方法采用DSRFG法[10]，相对于以往的方法，DSRFG法具有以下优点：该方法基于严格的理论推导，能产生满足任意形式功率谱及各向异性的湍流脉动风速场，能够严格满足流体连续性条件，从而保证了大涡模拟计算的稳定性，易于并行化处理等。采用DSRFG方法生成的入口瞬时风速分布云图，如图7所示。

图7 DSRFG法产生的入口瞬时速度场Fig.7 Instantaneous velocity field at the inflow boundary generated by DSRFG method

计算域的入口处边界条件设定为速度入口，按照定日镜实际工作环境，大气边界层平均风速度剖面的模拟采用如式(8)的指数率形式

(8)

(a) 平均风速剖面图

(b) 湍流强度剖面图图8 数值模拟与风洞试验入流条件对比图Fig.8 Inflow boundary condition of numerical simulation and wind tunnel test

湍流强度参考日本规范[11]中的第II类地貌，取值如式(9)

(9)

式中：I0为由规范确定的湍流强度；zg为梯度风高度，m；B类风场下，I0=0.23，zg=350 m。数值模拟与风洞试验湍流强度剖面图如图8(b)。

2.5 有限元模型

通过ANSYS软件选用不同的材料单元模拟了定日镜主要组成构件，建立了定日镜有限元模型。采用shell63单元模拟镜面和转动轴、立柱；采用beam188单元模拟出不同型号空心钢管构成的后部支撑钢架；利用耦合的方法，模拟出支撑钢架与镜面和镜面与镜面之间的驳接抓连接。通过命令流，可快速实现有限元模型不同工况的转换。图9给出了镜面仰角为0°的有限元模型示意图。实测频率由业主方提供，实测试验在24块镜面中点、立柱中点和顶端各布置一个加速度传感器，采用环境激励，通过模态分析软件将采集数据进行处理分析，得到定日镜实体结构的模态参数。对有限元模型进行模态分析，将有限元模型前三阶自振频率与实测频率进行对比，结果如表1。将Fluent流体计算得到的定日镜表面风荷载时程加载到有限元模型上进行瞬态动力分析，得到定日镜各测点的位移响应，位移测点的布置如图10。

图9 有限元模型Fig.9 Finite element model of heliostat

图10 位移测点布置Fig.10 Displacement monitoring points on the FEM

表1 不同镜面仰角下有限元模型与实体结构的自振频率Tab.1 Natural frequency of heliostat finite element model and Entity structure different elevation angle Hz

3 结果分析

3.1 流场分析

风洞试验得到的流场信息十分有限，而数值模拟方法可以得到丰富的可视化流场信息，可以更好探究流场作用机理。

图11为两种模拟方法的各工况下在Y=0、Z=0.15H和Z=0.85H平面上的平均风速流线图。从图中可以看出，两种模拟方法都表现出相似的流场特性，这反映了定日镜结构在不同镜面仰角下流场的一般规律。气流在定日镜边缘发生流动分离，涡旋作用于定日镜背面形成负压区。随着仰角增大，定日镜后方的两侧尾流区域变长且集中。这在定日镜群中，可能会对周边定日镜产生不利影响。

(a) LES-00

(b) DES-00

(c) LES-30

(d) DES-30

(e) LES-60

(f) DES-60图11 各工况下平均风速流线图Fig.11 Mean veloicty contours around heliostat with different elevation angle

图12为两种模拟方法的各工况下Y=3处涡量图。0°和30°时，气流冲击定日镜后，气流从流场驻点处分别向定日镜顶端和底端流去，部分气流在镜面顶端发生分离，涡旋向上发展，再附着现象不明显；定日镜底部也相似的现象。随着仰角增大，60°时，气流在底部分离，发生强烈的涡旋脱落，形成众多强度较大的小涡后再附着在镜面的中上部。与DES相比，LES在结构近壁面形成的涡旋尺度更加细碎。

(a) LES-00

(b) DES-00

(c) LES-30

(d) DES-30

(e) LES-60

(f) DES-60图12 各工况下Y=3 m平面涡量图Fig.12 Vorticity magnitude contours on X-Z plane at Y=3 m with different elevation angle

3.2 风压相关系数

在对建筑结构风致响应的考虑中，脉动风压的空间相关性起这重要的作用。两点之间的脉动风压相关系数计算式为

(10)

式中：Cov(X,Y)表示两点脉动风压协方差；Var[X]和Var[Y]分别表示两点脉动风压方差。

以定日镜中间B列测点为代表进行比较，选取了8个测点来研究脉动风压随高度变化的规律，分析定日镜表面脉动风压竖向相关性。不同工况的脉动风压相关系数如图13。

图13 B列各测点竖向相关性Fig.13 Vertical correlation of column B

迎风面脉动风压的竖向相关系数随着测点间距离的增大而减小。背风面的相关性要强于迎风面，背风面测点的相关性均在0.9以上。随着镜面仰角的增大，迎风面相关系数逐渐增大，背风面相关系数逐渐减小。对比发现，DES结果的相关性均强于LES的结果。

3.3 风荷载

风荷载时程具有随机性，对风荷载时程进行谱分析可以较直观的看出风荷载的特性。图14为定日镜0°工况镜面总体、顶层和底层的归一化风荷载功率谱的对比图。定日镜属于低矮结构，其风荷载谱受特征湍流影响较大，即定日镜表面涡旋脱落的影响。从图中可以看出，在低频段时，风荷载功率谱与卡门谱有吻合度较好，而风荷载功率谱峰值出现0.7 Hz， 说明定日镜顶部与底部的涡旋脱离频率集中在0.7 Hz附近的低频段。风荷载功率谱与卡门在1 Hz附近出现分离，并且衰减较快。这在Yan等[12]和Daniels等[13]的研究也有相似的情况，这是由于数值模拟方法有限的网格分辨率造成的[14]。

图14 顺风向荷载时程及功率谱Fig.14 Time history and power spectrum of along-wind load

3.4 等效静力风荷载

为验证数值风洞模拟的准确性，对模拟结果进行时域分析和频域分析计算出等效静力风荷载，并与风洞试验测力结果对比，以验证数值模拟结果的准确性。等效静风荷载的计算方法有多种，本文分别采用惯性风荷载法、惯性力-荷载响应相关法[15]对定日镜的等效静力风荷载进行时域分析和频域分析。

惯性风荷载方法中等效静力风荷载可表示为

(11)

惯性力-荷载响应相关法采用定日镜前三阶振型对应的自振频率进行计算，通过编写Matlab程序，对等效静力风荷载进行求解计算。

表2为各工况的等效风荷载计算结果。时域法和频域法结果比较接近，但和风洞结果相比都较小。这是由于数值模拟测点设置有限，只考虑了作用于定日镜镜面上的风荷载，忽略了支撑臂、转动轴和立柱的风荷载导致的。频域分析结果小于时域分析，这是由于频域法只考虑了三阶振型。LES和DES的结果差距较小，表面这两种湍流模型对等效风荷载的预测有较好的一致性。

表2 定日镜顺风向等效风荷载数值模拟与风洞试验对比Tab.2 Comparison between numerical simulation and wind tunnel test of equivalent along-wind load of heliostat N

3.5 位移响应分析

各工况不同结点的位移时程趋势基本相同，图15为LES0°镜面仰角下的各测点的位移时程。从图中可以看出，定日镜底部的镜面位移最小，顶部镜面的位移最大。

图15 各测点位移时程Fig.15 Displacement time history of each monitoring point

位移时程数据具有随机性，对定日镜各工况选取对应测点的顺风向位移响应谱进行分析，可以更好地揭示其结构动力反映机理。各工况对应结点归一化位移响应谱如图16，位移响应谱共振峰值能量如表3。可以看出，第一阶振动模态对定日镜的风致响应起主要作用。对于定日镜上部RT测点对应一阶模态的共振峰值较为明显，随着仰角增大，定日镜上部RT测点的一阶共振峰值表现出更高的能量。0°和30°时，定日镜下部存在两个振动峰值，第一个振动峰值对应一阶模态，第二个振动峰值对应三阶模态。随着仰角增大，定日镜底部第二个峰值逐渐减弱。从图11可以看出，由于仰角较小时，气流流向定日镜表面后，较多的气流从定日镜底端绕流，在底部与地面之间形成丰富的旋涡，发生剧烈的旋涡脱落现象，引起多阶共振响应，动力响应较为强烈。而随着仰角增大，气流在定日镜底端分离，脱落漩在镜面中上部再附着，上部的共振峰值表现出更高的能量。两种模拟方法对功率谱能量峰值频率的预测有较好的一致性，而LES比DES在定日镜上部和下部的共振峰值预测出更高的能量。如图17所示，这是由于LES比DES在定日镜近壁面形成更多尺度细碎的小涡，并再附着在定日镜背风面，使LES预测出更高的脉动能量造成的。

图16 定日镜各工况归一化位移功率谱Fig.16 Normalized displacement power spectrum of heliostat with different elevation angle

(a) LES-00

(b) DES-00图17 Y=3 m的X-Z剖面上LES和DES涡量分布图Fig.17 Vorticity magnitude contours on X-Z plane at Y=3 m by LES and DES

表3 各工况归一化位移功率谱共振峰值Tab.3 Peak value of normalized displacement power spectrum with different elevation angle

3.6 风振系数

进行结构风致响应分析可以得到设计人员关注的风振系数，其反映的是脉动风对结构的影响。各测点位移风振系数可由下式求出

(12)

定日镜各部位风振系数随镜面仰角而变化，对定日镜各工况选取对应测点的风振系数进行分析，得到定日镜不同测点各工况风振系数如图18。镜面仰角为0°时，定日镜底端有丰富的涡旋脱落现象，脉动响应影响较大，此时下部的风振系数最大，约为3.1。随着镜面仰角增大到60°，在镜面底端分离的涡旋再附在定日镜中上部，中部的风振系数增大到2.0，上部的风振系数增大到3.4。

图18 定日镜不同测点各工况风振系数Fig.18 Wind vibration coefficient of heliostat with different elevation angle

4 结 论

本文采用LES和DES，对0°风向角下0°、30°、60°镜面仰角的多种定日镜典型工况进行数值风洞模拟计算，利用ANSYS软件建立有限模型，得到定日镜结构的动力特性，将模拟结果与风洞结果进行对比分析，主要结论如下：

(1) 运用LES和DES可以较准确地预测定日镜顺风向等效风荷载。LES的结果更接近风洞试验结果，但两种模拟方法差距不大，其结果差距在2%以内。

(2) 两种方法模拟的流场表现出相同的规律。随着仰角的增大，定日镜底部的涡旋出现脱离再附着现象，且尾迹区发展得更狭长且集中。这在定日镜群中，可能会对后方定日镜产生不利影响。

(3) 随着仰角增大，流场对定日镜影响发生变化，定日镜下部的共振峰值能量逐渐减小，上部的共振峰值能量逐渐增大。两种方法的位移响应谱都表现出相同的趋势，但LES在结构近壁面的涡旋尺度更加细碎，对定日镜的共振峰值预测出更高的能量。

(4) 定日镜风振系数各部位最不利工况随镜面仰角变化。镜面仰角为0°时，下部的风振系数最大为3.1。中上部的最不利工况发生在镜面仰角为60°时，风振系数分别为2.0、3.4。

(5) 对于竖向相关系数，两种方法的模拟结果都有相同的变化趋势，DES结果的相关性均强于LES的结果。

(6) 综合本文的模拟结果，结合DSRFG与LES、DES的模拟方法，可以有效获得定日镜各测点的风荷载时程。结合有限元分析，可为定日镜的结构抗风设计提供初步的参考。