0 引言
复杂山地流场受地形和大气条件影响，存在高湍流、分离及大气环流等现象，因此准确刻画复杂山地流场特性是一个巨大的挑战。由于具有丰富风资源的平坦地形区域较稀缺、海上风电机组建设和输电成本高等原因，复杂山地风电场的建设已成为一大趋势。由于高塔筒风力机具有度电成本低、发电量高的优势，目前复杂山地风力机也向大型化发展。随着大型风电机组的建立，复杂山地风电场风资源特性也亟待研究。
目前风工程中对风资源的刻画主要依赖各种大气边界层解析模型，包括平均风速、湍流强度和湍流能量等参数的工程应用模型。对于平均风速通常以风速廓线描述，常用的风速廓线表达式有对数律和指数律模型。为考虑大气稳定度对风速廓线的影响，学者们利用莫宁−奥布霍夫（Monin-Obukhov）相似性理论对这种模型提出了可行性修正方法[1-2]。为使近地层和艾克曼（Ekman）层的风速廓线具有良好的连续性形态，GRYNING等[3]通过引入不同高度边界层的长度标度参数，重新定义了对数律模型。在风速廓线模型的研究中，LU等[4]、LIU等[5]考虑地转偏向力对风向的影响，结果表明科氏力对平均风速和湍流强度在三个正交方向上的垂直廓线的影响不可忽略，且其对风速主流向的改变比其他两个正交方向更为明显。NARASIMHAN等[6]基于莫宁−奥布霍夫理论和艾克曼层大气边界层垂直结构特征，发展了一种适用于常规中性和稳定大气条件下的风速廓线解析模型，该模型考虑了大气边界层厚度和地转阻力。
对于湍流能量，最早由KOLMOGOROV[11]提出了Kolmogorov能量谱，此后又衍生出多种形式的脉动风谱，既有能量谱，也有功率谱。较为经典的脉动风功率谱密度模型有von Kármán谱模型[12]、Davenport谱模型[13]、Kaimal谱模型[14]等。为提高脉动风谱工程应用的精度，LARSÉN等[15]利用丹麦两个站点（一个在沿海，一个在海上）的大量平均气象数据和高频声波风速计数据，研究了边界层风全尺度频谱，结果表明谱隙存在且可以模拟。CHENG等[16]通过量纲分析和对以往理论的大量修正，提出了由小弗劳德数下的浮力子范围、过渡区和各向同性惯性子区范围三部分组成的湍流动能和温度谱形状。GYATSO等[17]基于青藏高原外场的大量实测风速，建立了高海拔流向脉动风速谱模型，准确地反映了高海拔地区大气低频湍流特性。
尽管上述大气边界层解析模型在应用过程中考虑了大气稳定度、地转偏向力及艾克曼层特性等，但其在复杂山地流场中的适应性尚不明确。为精细化评估复杂地形风电场流动特性，结合测风塔实测数据和兼顾中尺度效应和微尺度流动细节的高精度仿真数据对经典大气边界层解析模型在复杂地形流场中适应性开展研究非常重要。
本文以国际经典案例Askervien山为应用对象，将实测数据和基于天气研究与预报模式（weather research and forecasting model, WRF）耦合计算流体动力学（computational fluid dynamics, CFD）方法得到的高精度数值仿真数据作对比，研究经典的风速廓线模型、湍流强度廓线模型、脉动风谱模型等在复杂山地流场的适应能力。其中，风速廓线模型选择对数律、指数律和Gryning三种模型，湍流强度廓线选择对数律和指数律模型，脉动风谱模型选择常用的von Kármán谱、Davenport谱和Kaimal谱。研究结果可进一步用于复杂地形流场大气边界层解析模型建模及精细化风资源评估。
1 Askervein山案例
Askervein山位于英国苏格兰南斯特岛西海岸，海拔为116 m，且几乎呈椭圆形，长轴大致为西北−东南方向，短轴大致为西南−东北方向（长轴约2 000 m，短轴约1 000 m）。南斯特岛长年刮西南风，即风向大致沿Askervein山的短轴方向，同时，短轴方向的上游地形比较平坦，与海岸一致，山及其周围地表被短草和低矮灌木覆盖，因此地表粗糙度较小。Askervein山概貌如图1所示。
1.1 实测数据
Askervein山实验是国际能源署（International Energy Agency, IEA）组织的针对复杂山地进行的经典基准测试案例，由TAYLOR和TEUNISSEN等于1982年9月至1983年10月期间进行测量完成，数据于1983年和1985年发布[18-19]。在观测实验中，如图2所示，沿A线、AA线和B线分布了约40个测风塔，这三条线的交点位于山顶（hill top, HT）和中心点（center point, CP）。图2中，RS（reference site）为上游参考站，用于给定来流；ANE40（A线东北侧，海拔40 m）为背风坡测风塔，取其10 m离地风速。因该实验测风数据完整度很高，是用来验证大气边界层解析模型适用性最合适的资料之一。
1.2 高精度数值仿真数据
尽管Askervein山实验得到了很多测风数据，但是这些数据都是单点测量，且测风高度较低，无法准确反映复杂山地局地流场和高层风特性。为此，本文同时采用中−微尺度WRF-CFD耦合方法针对Askervein山进行了数值仿真计算，结合高精度的数值仿真数据，对经典大气边界层解析模型做全面的适应性对比分析与验证。
1.2.1 中尺度气象模式WRF设置
首先采用气象模式WRF进行中尺度计算。其中气候输入数据是欧洲中心天气预报中心（European Centre for Medium-Range Weather Forecasts, ECMWF）的第五代数据（ERA5），其空间分辨率为0.25°，时间分辨率为1 h。地形类别和地形高度静态数据取自美国地质勘探局（United States Geological Survey, USGS）的数据，其空间分辨率为30''。中尺度数值模拟时间段为1983年10月02日00时至1983年10月04日00时，采用四层嵌套，即图3中的D01、D02、D03、D04，投影方式设为兰勃特投影，模拟区域的中心点位于（57.183°N，7.37°W）。如图3，空间分辨率分别为27、9、3和1 km，时间分辨率分别为60、20、6.67和2.22 s，垂直高度为56层，靠近地面附近进行加密，其中100 m以下有10层，100 ~ 500 m范围有10层，500 ~ 1 000 m范围有5层，1 000 m以上有31层。模拟结果输出时间间隔为10 min。模拟过程中采用的物理参数化方案包括微物理过程（Thompson）、长波和短波辐射（RRTMG）、近地层（Monin-Obukhov）、陆面过程（Unified Noah）、行星边界层（MYJ）及积云对流参数（Kain-Fritsch）等[20-22]。
1.2.2 微尺度计算流体力学模式CFD设置
对于Askervein山案例，CFD计算域高度取2 000 m，水平域尺寸取6 km × 6 km，垂直分层生成网格。在覆盖山体的近地范围细化网格，近地层网格水平分辨率为5 m，垂直分辨率为2.5 m。背景网格水平分辨率为20 m，垂直分辨率为10 m。计算域网格如图4所示。
2 经典大气边界层解析模型
2.1 风速廓线模型
2.1.1 对数律模型
经典的对数律风速廓线模型是由BLACKADAR等[25]于1968年提出，该模型是在中性大气条件下基于长度尺度参数提出的一阶闭合模型。表达式为：
$U\left( z \right)=\frac{{{u}_{*}}}{\kappa }\ln \left( \frac{z}{{{z}_{0}}} \right)$
式中：u*为地表摩擦速度，可通过给定高度的参考风速得到；z为离地高度；κ为 von Kármán常数，0.41；z0为地表粗糙度。
考虑大气稳定度后，该模型可进一步发展为[2]：
$U\left( z \right)=\frac{{{u}_{*}}}{\kappa }\ln \left[ \frac{z}{{{z}_{0}}}-{{\psi }_{\text{M}}}\left( \frac{z}{L} \right) \right]$
式中：ΨM为大气稳定度修正项，ΨM表达式如下：
${{\psi }_{\text{M}}}=\left\{ \begin{matrix} -4.7\left( {z}/{L}\; \right) & L>0 & \\ 0 & \left| L \right|\to \infty & \\ f\left( \xi \right) & L<0 & \\ \end{matrix} \right.$
$\begin{align} & f\left( \xi \right)=2\ln \left[ \frac{1}{2}\left( 1+\xi \right)\left( 1+\xi \right) \right]+\ln \left[ \frac{1}{2}\left( 1+{{\xi }^{2}} \right) \right]- \\ & \ \ \ \ \ \ \ \ \ \ \ 2arctan\left( \xi \right)+\frac{\text{ }\!\!\pi\!\!\text{ }}{2} \\ \end{align}$
式中：$\xi ={{\left[ 1-15\left( {z}/{L}\; \right) \right]}^{1/4}}$；L为Monin-Obukhov长度，表达式为：
$L={-\left( u_{*}^{3}\theta \right)}/{\left[ \kappa g\left( \overline{{w}'{\theta }'} \right) \right]}\;$
式中：θ为当地位温；g为重力加速度；$\overline{{w}'{\theta }'}$为当地热通量。
2.1.2 指数律模型
指数律风速廓线是除对数律模型外应用最广泛的风速廓线模型，表达式为：
$\frac{U\left( {{z}_{2}} \right)}{U\left( {{z}_{1}} \right)}={{\left( \frac{{{z}_{2}}}{{{z}_{1}}} \right)}^{\alpha }}$
式中：z1和z2分别为风速廓线计算的两层高度；$\alpha $为风速廓线指数。
考虑大气稳定度和地表粗糙度之后，风速廓线指数$\alpha $表达式如下[2]：
$\alpha =\frac{{{\phi }_{\text{M}}}\left( {{\bar{z}}}/{L}\; \right)}{\left[ \ln \left( {{\bar{z}}}/{{{z}_{0}}}\; \right)-{{\psi }_{\text{M}}}\left( {{\bar{z}}}/{L}\; \right) \right]}$
式中：$\bar{z}$为几何平均高度，定义为$\bar{z}={{\left( {{z}_{1}}\cdot {{z}_{2}} \right)}^{1/2}}$；ΨM定义与对数律一致；${{\phi }_{\ M}}$定义如下：
${{\phi }_{\text{M}}}=\left\{ \begin{matrix} 1+4.7\left( \frac{{\bar{z}}}{L} \right) & L>0 & \\ 1 & \left| L \right|\to \infty & \\ {{\left[ 1-15\left( \frac{{\bar{z}}}{L} \right) \right]}^{-1/4}} & L<0 & \\\end{matrix} \right.$
2.1.3 Gryning模型
GRYNING等[3]指出，基于近地层理论和莫宁−奥布霍夫尺度的风廓线理论仅在50 ~ 80 m高度内有效。随着风力机大型化的发展，100 ~ 300 m范围内的风速特性亟待研究。为兼顾近地层和艾克曼层，GRYNING等根据高度范围，引入近地层长度分量LSL、中间层长度分量LML和上层长度分量LUL，最终确定风速廓线模型如下：
$U\left( z \right)=\left\{ \begin{matrix} {{f}_{1}}\left( \frac{z}{L} \right) & L>0 & {} \\ {{f}_{2}}\left( \frac{z}{L} \right) & \left| L \right|\to \infty & {} \\ {{f}_{3}}\left( \frac{z}{L} \right) & L<0 & {} \\ \end{matrix} \right.$
${{f}_{1}}\left( \frac{z}{L} \right)=\frac{{{u}_{*}}}{\kappa }\left[ \begin{align} & \ln \left( \frac{z}{{{z}_{0}}} \right)+\frac{5z}{L}\left( 1-\frac{z}{2{{z}_{\text{i}}}} \right)+ \\ & \frac{z}{{{L}_{\text{ML}}}}-\frac{z}{{{z}_{\text{i}}}}\left( \frac{z}{2{{L}_{\text{ML}}}} \right) \\ \end{align} \right]$
${{f}_{2}}\left( \frac{z}{L} \right)=\frac{{{u}_{*}}}{\kappa }\left[ \ln \left( \frac{z}{{{z}_{0}}} \right)+\frac{z}{{{L}_{\text{ML}}}}-\frac{z}{{{z}_{\text{i}}}}\left( \frac{z}{2{{L}_{\text{ML}}}} \right) \right]$
${{f}_{3}}\left( \frac{z}{L} \right)=\frac{{{u}_{*}}}{\kappa }\left\{ \begin{align} & \ln \left( \frac{z}{{{z}_{0}}} \right)-\varphi \left( \frac{z}{L} \right)-\frac{z}{{{z}_{\text{i}}}}\left( \frac{z}{2{{L}_{\text{ML}}}} \right)+ \\ & \frac{z}{{{z}_{\text{i}}}}\left[ 1+\frac{{{\left( 1-12\frac{z}{L} \right)}^{2/3}}-1}{8\frac{z}{L}} \right]+\frac{z}{{{L}_{\text{ML}}}} \\ \end{align} \right\}$
式中：zi为边界层高度，定义为大气边界层流动不受地面应力影响的高度，此时流动转变为地转风；$\varphi \left( \frac{z}{L} \right)$定义为：
$\begin{align} & \varphi \left( \frac{z}{L} \right)=\frac{3}{2}\ln \left( \frac{1+a+{{a}^{2}}}{3} \right)- \\ & \ \ \ \ \ \ \ \ \ \ \ \sqrt{3}\arctan \left( \frac{1+2a}{\sqrt{3}} \right)+\frac{\pi }{\sqrt{3}} \\ \end{align}$
式中：$a={{\left[ 1-12\left( {z}/{L}\; \right) \right]}^{1/3}}$。
LML的定义为：
$\frac{{{u}_{*}}}{{{f}_{\text{c}}}{{L}_{\text{ML}}}}=\left\{ \begin{matrix} -2\ln \left( \frac{{{u}_{*}}}{{{f}_{\text{c}}}{{z}_{0}}} \right)+55 & \\ \left[ -2\ln \left( \frac{{{u}_{*}}}{{{f}_{\text{c}}}{{z}_{0}}} \right)+55 \right]{{e}^{-\frac{{{\left( {{{u}_{*}}}/{{{f}_{\text{c}}}L}\; \right)}^{2}}}{400}}} & \\ \end{matrix} \right.$
式中：fc为科氏力系数，fc = 2ωsinθlat；ω为地球自转平均角速度，ω = 7.292 × 10−5 rad/s；θlat为当地纬度。
2.2 湍流强度廓线模型
2.2.1 对数律模型
湍流强度可反映脉动风速的相对强度，WIERINGA[26]根据对数律风速廓线模型给出了中性大气边界条件下的湍流强度廓线模型：
${{I}_{u}}\left( z \right)=\frac{{{\sigma }_{u}}}{U\left( z \right)}=\frac{{{\sigma }_{u}}}{{{u}_{*}}}\kappa {{\left[ \ln \left( \frac{z}{{{z}_{0}}} \right) \right]}^{-1}}$
式中：${{\sigma }_{u}}$为风速标准差。其他各参数含义参考对数律风廓线模型。
考虑大气稳定度之后，该模型可以写为：
${{I}_{u}}\left( z \right)=\frac{{{\sigma }_{u}}}{U\left( z \right)}=\frac{{{\sigma }_{u}}}{{{u}_{*}}}\kappa {{\left[ \ln \left( \frac{z}{{{z}_{0}}} \right)-{{\psi }_{\text{M}}}\left( \frac{z}{L} \right) \right]}^{-1}}$
式中：ΨM取值参考式（3）。
2.2.2 指数律模型
湍流强度廓线的指数律模型公式如下：
${{I}_{u}}\left( z \right)={{I}_{u}}\left( {{z}_{\text{ref}}} \right){{\left( \frac{z}{{{z}_{\text{ref}}}} \right)}^{-\alpha -0.05}}$
式中：zref 为参考高度；Iu (zref)为参考高度处湍流强度。
2.3 脉动风功率谱模型
脉动风功率谱密度是描述湍流风速的功率谱密度函数，是描述风速脉动特性的一种方法。在风工程中已有很多脉动风功率谱密度的经验表达式，其中应用较广的是与纵向相关的von Kármán谱、Davenport谱、Kaimal谱等。
2.3.1 von Kármán谱
$\frac{n{{S}_{u}}\left( z,n \right)}{\sigma _{u}^{2}}=\frac{4x}{{{\left( 1+70.8{{x}^{2}} \right)}^{5/6}}}$
式中：n为频率；x为无量纲频率或莫宁坐标，$x={n\cdot L_{u}^{x}}/{U\left( z \right)}\;$；$L_{u}^{x}$为纵向湍流积分尺度，取$L_{u}^{x}=100\sqrt{z\text{/30}}$。
2.3.2 Davenport谱
加拿大科学家DAVENPORT[13]根据世界不同地点、不同高度实测得到的70多次风速超过9 m/s的强风下的纵向湍流功率谱实测值取平均值后，建立了Davenport谱模型。Davenport谱是一种与高度无关的谱，经验表达式为：
$\frac{n{{S}_{u}}\left( z,n \right)}{\sigma _{u}^{2}}=\frac{2}{3}\frac{{{x}^{2}}}{{{\left( 1+{{x}^{2}} \right)}^{4/3}}}$
式中：x = 1200·n/U(10)。
2.3.3 Kaimal谱
Davenport谱模型在早期的风速模拟过程中应用广泛，但由于没有考虑高度的影响，且粗糙度系数的取值缺少明确的解析表达式，因此会在风速模拟结果中引入较大的误差，降低了该模型的可用性。为改进Davenport谱不足，KAIMAL等[14]基于美国堪萨斯州平坦均匀地形上的测量结果，建立了Kaimal谱：
$\frac{n{{S}_{u}}\left( z,n \right)}{\sigma _{u}^{2}}=\frac{100}{3}\frac{x}{{{\left( 1+50x \right)}^{5/3}}}$
式中：x = n·z/U(z)。
3 大气边界层解析模型在复杂山地流场应用中的适应性分析
3.1 平均风廓线
表1 Askervein山RS点风速误差指标统计Table 1 Statistics on wind speed error indexes at the RS point of the Askervein hill |
| 高度 | WRF-CFD | ||
|---|---|---|---|
| EMA/(m/s) | ERMS/(m/s) | PCR | |
| 50 m以下 | 0.980 8 | 1.136 9 | 0.974 4 |
| 50 m以上 | 0.443 3 | 0.580 3 | 0.987 6 |
表2 Askervein山HT点风速误差指标统计Table 2 Statistics on wind speed error indexes at the HT point of the Askervein hill |
| 高度/m | WRF-CFD | |
|---|---|---|
| EAB/(m/s) | ERB/% | |
| 8 | −0.82 | −5.01 |
| 15 | −0.36 | −2.19 |
| 24 | 0.57 | 3.52 |
| 34 | 1.61 | 10.22 |
${{E}_{\text{MA}}}=\frac{\text{1}}{N}\sum\limits_{i=1}^{N}{\left( \left| {{v}_{\text{s},i}}-{{v}_{\text{e,}}}_{i} \right| \right)}$
${{E}_{\text{AB}}}={{v}_{\text{s},i}}-{{v}_{\text{e},i}}$
${{E}_{\text{RB}}}=\frac{{{v}_{\text{s},i}}-{{v}_{\text{e},i}}}{{{v}_{\text{e},i}}}\times 100%$
${{E}_{\text{RMS}}}=\sqrt{\frac{\text{1}}{N}\sum\limits_{i=1}^{N}{{{\left( {{v}_{\text{s,}}}_{i}-{{v}_{\text{e},}}_{i} \right)}^{\text{2}}}}}$
${{P}_{\text{CR}}}=\frac{\sum\nolimits_{i\text{=1}}^{N}{\left( {{v}_{\text{s,}}}_{i}-{{v}_{\text{s,}}}_{\text{avg}} \right)\cdot \sum\nolimits_{i\text{=1}}^{N}{\left( {{v}_{\text{e,}}}_{i}-{{v}_{\text{e,}}}_{\text{avg}} \right)}}}{\sqrt{\sum\nolimits_{i\text{=1}}^{N}{{{\left( {{v}_{\text{s,}}}_{i}-{{v}_{\text{s,}}}_{\text{avg}} \right)}^{\text{2}}}\cdot \sum\nolimits_{i\text{=1}}^{N}{{{\left( {{v}_{\text{e,}}}_{i}-{{v}_{\text{e,}}}_{\text{avg}} \right)}^{\text{2}}}}}}}$
式中：vs,i为第i个点的风速模拟值；vs,avg为所有点的风速模拟平均值；ve,i为第i个点的风速实验值；ve,avg为所有点的风速实验平均值；N为观测点个数。
虽然Gryning模型在RS点比对数律和指数律模型表现较优，但是在高空与实测数据及高精度仿真数据仍有一些偏差。因此，在艾克曼层还需采用更复杂的风速廓线模型，在高层需要考虑非恒定的动量通量和摩擦速度等[28]。
结合山前平坦区域（RS点）和山顶（HT点）平均风速廓线分析，发现基于混合长理论的Gryning模型可以较准确地反映平坦地形风速廓线特征，但对于山顶风加速效应，Gryning模型无法描述；基于近地层理论和莫宁−奥布霍夫理论的对数律和指数律模型可以用于描述平坦地形50 m以下的近地层风速廓线，对于艾克曼层和风加速效应均无法反映。
3.2 湍流强度廓线
图6展示了Askervein山山前平坦区域RS点处的平均湍流强度廓线。其中，对数律湍流强度廓线模型采用式（15）拟合，σu取0.862 4 m/s，其余参数与对数律风速廓线的取值一致；指数律湍流强度廓线模型采用公式（17）拟合，式中zref 取40 m，对应的Iref 为0.077，风速廓线指数$\alpha $按式（7）计算。
由图6可以看到，在50 m以下范围采用的中/微尺度WRF-CFD耦合方法得到的数值结果与测量数据基本一致。在低空层（约50 m高度以下），对数律和指数律湍流强度廓线模型拟合效果相当，且与实测数据及高精度的WRF-CFD耦合仿真数据基本一致。但是随着高度增加（50 m高度以上），对数律和指数律湍流强度廓线模型与WRF-CFD耦合仿真结果的偏差增大：WRF-CFD耦合中，湍流强度随着高度增加具有显著衰减趋势，而在对数律和指数律湍流强度廓线模型中这一衰减不明显。
3.3 脉动风功率谱密度
虽然已有的谱方法可以刻画部分频率段（惯性子区）的能量，但无法准确刻画复杂山地高频段和低频段流场信息，存在高估高频区能量和低估低频区耗散的问题。因此，为建立适合复杂地形中脉动风特性的全频率功率谱模型，还需分别考虑地形对高频段和低频段的影响，需要考虑其他幂次衰减指数。
4 结论
以国际经典案例Askerven山为对象研究了经典大气边界层解析模型在复杂山地流场应用适应性问题。基于实验观测数据和WRF-CFD高精度仿真数据研究了经典的风速廓线、湍流强度廓线和脉动风谱等模型在复杂山地山前平坦区域、山顶及背风坡等位置的风资源评估效果，得出以下主要结论：
（1） 基于WRF-CFD方法得到的高精度仿真数据与实验观测结果有较好的一致性，可以作为观测数据的补充用于评估经典大气边界层解析模型在复杂地形流场中应用适应性问题。
（2） 基于近地层理论和莫宁−奥布霍夫理论的对数律和指数律模型可以用于描述山前平坦区域50 m高度以下的近地层风速廓线，但对于艾克曼层和山顶风加速效应均都无法反映；基于混合长理论的Gryning模型可以较准确地反映平坦地区低层和高层的风速廓线特征，但对于山顶风加速效应，Gryning模型也无法描述。
（3）在50 m高度以下，对数律和指数律湍流强度廓线模型可以较准确地反映平坦区域近地层湍流强度垂直分布；但在50 m高度以上，对数律和指数律模型低估了湍流强度在垂直方向的衰减速率。
（4）von Kármán谱、Davenport谱和Kaimal谱可以较准确地评估山顶和背风坡脉动风速惯性子区的能量特征，但存在高频部分高估能量、低频部分低估脉动衰减的问题；相比von Kármán谱和Kaimal谱，Davenport谱在高频区域表现略优。
在复杂地形流场分析中，经典的大气边界层解析模型存在一定的局限性。在未来，针对复杂地形风电场大气边界层解析模型研究，可以有以下三个方向：一是在Gryning模拟中考虑风加速效应和非恒定的湍流通量；二是考虑Kaimal谱全频率段的模型；三是研究山体特征对于大气边界层解析模型的影响。