在WRF中怎么算风能密度
前言
2022 年，全国风能资源为正常略偏小年景。10 米高度年平均风速较近 10 年（2012 ～ 2021 年，下同）平均值偏小 0.82%，较2021 年偏小 0.96%。70 米高度年平均风速约 5.4m/s，年平均风功率密度约 193.1W/m2；100 米高度年平均风速约 5.7m/s，年平均风功率密度约 227.4W/m2。其中，湖北、江西、湖南、重庆较近 10 年平均值偏大，贵州、山西、宁夏、江苏、山东、河北、天津、内蒙古、西藏、河南、云南偏小，其他地区与近 10 年平均值接近。————《2022年中国风能太阳能资源年景公报》
现目前环保概念已日益重要，风电资源的大力发展是趋势。
近日在微信群见有人问风能问题，略感兴趣，尝试一二。
风能密度是单位迎风面积可获得的风的功率，与风速的三次方和空气密度成正比关系。
其公式为
from wrf import uvmet, to_np, getvar, interplevel, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords,destagger
import numpy as np
from netCDF4 import Dataset
import metpy.calc as mpcalc
import xarray as xr
from metpy.units import units
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from matplotlib.colors import from_levels_and_colors
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
from cartopy.feature import NaturalEarthFeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
import cmaps
from glob import glob
import metpy.calc as mpcalc
import metpy.constants as constants
# 打开wrfout文件
ncfile = Dataset("/home/mw/input/wrfout3385/wrfout_d02_2022-07-14_0800.nc")
# 提取温度和气压
t = getvar(ncfile, "T")
p = getvar(ncfile, "pressure")
qvor = getvar(ncfile,"QVAPOR")
z = getvar(ncfile,"z")
# 计算密度
rho = mpcalc.density(p, t, qvor)
# 提取U和V
u = getvar(ncfile, "U")
v = getvar(ncfile, "V")
# 计算风速和风向
ws = getvar(ncfile, "wspd_wdir", units="m s-1")[0]
# 计算风能密度
we = 0.5 * rho * ws**3
print("风能密度:", we)Warning: ecCodes 2.21.0 or higher is recommended. You are running version 2.14.1
风能密度: <xarray.DataArray (bottom_top: 49, south_north: 90, west_east: 90)>
<Quantity([[[ 484.84293 529.96545 504.5152 ... 4552.3765 3606.1038
3129.516 ]
[ 463.362 499.4131 503.2667 ... 4232.495 3103.547
2655.305 ]
[ 459.04135 529.82355 527.0069 ... 3670.1274 2776.1077
2290.109 ]
...
[ 514.2071 188.49896 56.96391 ... 739.167 731.23145
730.46674 ]
[ 720.0671 287.3647 61.945694 ... 1212.0271 1218.5209
1264.0342 ]
[1095.568 559.9499 162.13527 ... 1755.4768 1827.7842
2036.5676 ]]
[[ 755.2138 872.288 838.95056 ... 8328.032 6803.8213
5975.269 ]
[ 751.6843 870.4069 859.779 ... 7875.6426 5938.154
5125.1387 ]
[ 759.6554 908.22144 893.49146 ... 6947.0215 5293.724
4410.594 ]
...
[ 24.392471 24.961102 26.36342 ... 40.52036 42.305367
47.723076]
[ 23.207579 23.494274 24.84367 ... 41.84898 46.0941
51.06096 ]
[ 21.764877 22.430603 24.13408 ... 43.112915 48.1985
52.11022 ]]
[[ 457.6555 463.40555 466.3164 ... 486.48077 493.40723
503.26 ]
[ 457.30045 462.20395 462.56293 ... 498.935 505.3736
513.84357 ]
[ 447.29855 453.62625 457.2152 ... 509.33307 513.9056
522.51184 ]
...
[ 22.765352 23.764181 24.53554 ... 41.917824 50.805874
55.152233]
[ 22.878866 23.800423 24.347864 ... 40.85533 47.187817
52.54639 ]
[ 22.951187 23.529852 24.114712 ... 42.591198 48.30677
53.363426]]], 'kilogram / meter ** 3')>
Coordinates:
XLONG (south_north, west_east) float32 100.9 101.0 ... 111.0 111.1
XLAT (south_north, west_east) float32 27.24 27.24 ... 34.43 34.42
XTIME float32 480.0
Time datetime64[ns] 2022-07-14T08:00:00
wspd_wdir <U4 'wspd'
Dimensions without coordinates: bottom_top, south_north, west_east常见的风机都蛮高的，因为算风能的都盯着七十米到一百米这块
下面取高度一百米的风能密度并可视
we100 = interplevel(we,z,100)
we100.plot()
看来直接插值插到地形里去了，全是nan
那只好问问神奇海螺如何计算离地高度一百米。
神奇海螺说有两种办法，那我们就都试试
方法一
ph = getvar(ncfile, "PH")
phb = getvar(ncfile, "PHB")
hgt = getvar(ncfile, "HGT")
P=ph+phb
P = destagger(P,0,meta=True)
gmp=P/9.81-hgt
we1 = interplevel(we, gmp ,100)
we1绘图效果如下
方法二
gmp1 = z - hgt
we2 = interplevel(we, gmp1 ,100)
we2绘图效果如下
小结
绘图可见两种方法相差非常小，作差之后也是如此。选取喜欢的使用吧。
diff= we1-we2
diff.plot()
完整代码与文件可回复”风能“查看
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2023-12-03，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录