一种建筑物可分辨尺度的城市三维风场快速分析方法与流程

1.本发明涉及风环境分析技术领域，尤其涉及一种建筑物可分辨尺度的城市三维风场快速分析方法。


背景技术：

2.随着城市化水平的提高以及城市建筑密度、人口密度、机动车保有量等不断加大，城市区域局地风灾环境以及城市冠层内风场与污染扩散等逐渐成为人们日益关注的问题。研究城市冠层内风环境对城市建设、城市风灾、城市空气质量改善均具有重要作用。
3.过去对城市风的研究，一般忽略城市建筑物的影响，然而世界经济的发展催生大量的高层建筑，给高层建筑带来了一些列的风灾事故，这使得人们不得不重视建筑周围风环境会受到建筑影响的问题。建筑周围的风环境与建筑单体的外形、尺寸、建筑间距及周围的地形地貌等因素都有密切关系，由于建筑物钝体的阻挡作用，造成建筑物周围风场在时间和空间上都具有非常复杂的非定常流性状。如何更好地观测和模拟城市风环境，一直是城市微气候研究的难题。
4.当前研究建筑周围风环境的方法分为实地测试、物理模拟和计算机数值模拟。三种方法各有优缺点，在现今的研究或应用当中均有所体现。总体而言，当前研究方法共同存在的问题是研究成本(包括计算成本和实验成本)较高，无法长时间模拟或快速模拟建筑物周围风场。在一些特殊情形下，比如突发事故，三种研究方法无法满足实际研究需求。并且城市建筑物中还有大量为斜坡型屋顶结构，在目前的城市风模拟中未充分考虑坡度影响。因此，针对上述问题，有必要提出进一步地解决方案。


技术实现要素：

5.本发明旨在提供一种建筑物可分辨尺度的城市三维风场快速分析方法，以克服现有技术中存在的不足。
6.为解决上述技术问题，本发明的技术方案是：
7.一种建筑物可分辨尺度的城市三维风场快速分析方法，包括：
8.获取目标区域气象数据，所述气象数据包括气象数据时间和入流风信息；
9.获取目标区域建筑物数据，并确认各个建筑物在模型网格内的尺寸、相对位置及每个建筑物的迎风区、背风区和顶部涡旋区的尺度；
10.根据气象数据和建筑物数据构建初始风场；
11.根据质量守恒方程调整所述初始风场，得到最终风场。
12.本发明的一个较佳实施例中，所述质量守恒方程中包括斜坡前方、后方、右侧及左侧边界条件系数，所述斜坡包括建筑物倾斜屋顶及地表斜坡。
13.本发明的一个较佳实施例中，在所述质量守恒方程中对于底面及建筑物采用第二类边界条件其他边界采用第一类自由边界条件λ＝0。
14.本发明的一个较佳实施例中，所述根据质量守恒方程调整所述初始风场，得到最
终风场，包括：
15.采用变分方法求解所述质量守恒方程，且α1＝α2＝0.5。
16.本发明的一个较佳实施例中，所述建筑物数据包括屋顶是否为坡度、建筑物高度、建筑物长度及建筑物最长边和正东方向夹角角度。
17.本发明的一个较佳实施例中，所述确认各个建筑物在模型网格内的尺寸，包括：
18.当入流风和建筑物的夹角θ不超过20
°
时，建筑物在模型网格内的尺寸w
eff
、l
eff
与其实际尺寸w、l相等；
19.当入流风和建筑物的夹角大于20
°
时，w
eff
＝wsinθ lcosθ，l
eff
＝wcosθ lsinθ；
20.其中，w为宽度，l为长度。
21.本发明的一个较佳实施例中，所述根据气象数据和建筑物数据构建初始风场，包括：
22.计算空旷区域和建筑物区域的入流风的垂直廓线；
23.对单体建筑物周围风场参数化；
24.对建筑街渠风场参数化。
25.本发明的一个较佳实施例中，对单体建筑物周围风场参数化，包括：
26.当入流风和建筑物的夹角θ大于20
°
时，采用迎风面的入流风速分解为平行迎风面分量和垂直迎风面分量进行插值。
27.本发明的一个较佳实施例中，对建筑街渠风场参数化，包括：
28.当s/h》2.5或1.4《s/h《2.4时，采用单体建筑物周围风场参数化方案；
29.当s/h《1.4时，采用sr参数化方案或mr参数化方案；
30.其中，h为建筑物高度，s为建筑间距离。
31.与现有技术相比，本发明的有益效果是：
32.本发明通过利用简单的气象输入以及基本的建筑物信息能够快速模拟出超大城市中超高层建筑物对周围风场环境的影响，引入的参数化方案能够简化计算过程、减少计算量，同时加入坡度因素，使模拟结果更加准确。
附图说明
33.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
34.图1为本发明流程图；
35.图2为建筑物顶部涡旋区示意图；
36.图3为另一方案中建筑物顶部涡旋区示意图；
37.图4为入流风倾斜入射夹角示意图；
38.图5为建筑物等效长宽示意图；
39.图6为本发明的模拟网格图；
40.图7为本发明实施例中模拟区域建筑物分布图；
41.图8为本发明实施例中2-140m高度水平风场；
42.图9a为不同平面cedval风洞观测，图9b为本方法模拟风矢量图；
43.图10为不同垂直平面全风速廓线比较图，三角为cedval风洞观测数据，实线为本方法模式结果；
44.图11为不同垂直平面u、w风速廓线比较图，三角为cedval 风洞观测数据u，实线为本方法模式结果u，圆圈为cedval风洞观测数据w，虚线为本方法模式结果w。
具体实施方式
45.为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式对本发明进行进一步的详细描述。需要说明的是，在不冲突的情况下，本技术的实施例及实施例中的特征可以相互组合。
46.在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述的其他方式来实施，因此，本发明的保护范围并不受下面公开的具体实施例的限制。
47.如图1所示，一种建筑物可分辨尺度的城市三维风场快速分析方法，包括：
48.s100获取目标区域气象数据，所述气象数据包括气象数据时间和入流风信息。
49.s200获取目标区域建筑物数据，并确认各个建筑物在模型网格内的尺寸、相对位置及每个建筑物的迎风区、背风区和顶部涡旋区的尺度。
50.其中，所述建筑物数据包括屋顶是否为坡度、建筑物高度、建筑物长度及建筑物最长边和正东方向夹角角度，以在模拟网格中构建目标区域建筑物分布。
51.s300根据气象数据和建筑物数据构建初始风场。
52.s400根据质量守恒方程调整所述初始风场，得到最终风场。
53.具体而言，s300根据气象数据和建筑物数据构建初始风场，包括：
54.s310计算空旷区域和建筑物区域的入流风的垂直廓线，如下： (1)空旷区域
[0055][0056]
(2)建筑物区域
[0057][0058][0059]
在(1)-(3)中，h
can
表示冠层高度，u
can
表示冠层顶h
can
风速大小，z0表示粗糙系数，
一般约为0.1～0.2h
can
，根据一般城市冠层的粗糙度特征，以及现有城市模式的模拟经验，本技术取 z0＝0.1h
can
，d表示位移高度，一般约为0.7h
can
，本技术取d＝0.67h
can
。α(z)为衰减系数，为了简化模型，本技术将α(z)设为常数，α的范围一般为1～3，本技术取α＝1。α(z)和建筑物密度λf成正比，α＝kλf，其中k＝9.6。
[0060]
s320对单体建筑物周围风场参数化，如下：
[0061]
单体建筑物周围风场可分以下三个关键区域：迎风区(包括迎风位移区和迎风涡旋区)，尾流区(包括空腔区)，顶部涡旋区。
[0062]
s321建筑物前迎风区
[0063]
当入流风垂直于建筑物迎风墙面时，会在建筑物前迎风区域内形成位移区和迎风涡旋区。本模式采用的建筑物迎风区插值方案如下
[0064][0065][0066][0067][0068]
在(4)-(7)中，h表示建筑物高度，l表示建筑物迎风方向的长度，w表示建筑物垂直于迎风方向的长度。(4)(5)公式左项表示l
fx
、l
fx1
与h的比值。从水平方向上看，l
fx
表示迎风区(upwind zone) 长度，l
fx1
表示涡旋区长度，迎风区和涡旋区的空间大小近似看作 1/4椭球，椭球公式为(6)(7)，l
fx
即迎风区的椭球长轴，l
fx1
即迎风区内涡旋区的椭球长轴，椭球z方向的半轴为0.6h，y方向上半轴为0.5w。椭球的原点为建筑物迎风面底边中点。假设迎风区内涡旋区为a区域，迎风区为b区域。
[0069]
在涡旋区内，即a区域，当位置坐标(x,y,z)满足水平和垂直方向风速的插值方案如下：
[0070][0071][0072]
其中xf、zf表示x正方向、z正方向的位置，xf＝x，zf＝z。
[0073]
位移区指b-a区域，当位置坐标(x,y,z)满足且不在a区域内，水平方向风速的插值方案如下：
[0074]
u0(z)＝c
dz
u(z)(10)
[0075]
其中u(z)为s310中插值风速，c
dz
是一个衰减系数，根据一般城市冠层的粗糙度特征，本技术取c
dz
＝0.4。
[0076]
s322建筑物尾流区
[0077]
当入流风垂直于建筑物迎风墙面时，会在建筑物背后形成尾流区。本研究提供了两种背风区插值方案以供选择，通过风洞实验资料与模拟结果的对比试验，挑选合适方案。
[0078]
第一种背风区插值方案
[0079][0080]
lw＝3lc(12)
[0081][0082][0083]
在(11)-(14)中，h表示建筑物高度，l表示建筑物迎风方向的长度，w表示建筑物垂直于迎风方向的长度。(11)公式左项表示lc与h的比值。从水平方向上看，lw表示尾流区(wakezone)长度，lc表示空腔区(cavityzone)长度，尾流区和空腔区的空间大小近似看作1/4椭球，椭球公式为(13)(14)，lw即背风区的椭球长轴，lc即背风区内空腔区的椭球长轴，椭球z方向的半轴为h，y方向上半轴为0.5w。椭球的原点为建筑物背风面底边中点。假设尾流区内空腔区为a区域，尾流区为b区域。
[0084]
在空腔区内，即a区域，当位置坐标(x,y,z)满足水平方向风速的插值方案如下：
[0085][0086][0087]
在b-a区域，当位置坐标(x,y,z)满足且不在a区域内，水平方向风速的插值方案如下。
[0088][0089]
在(15)(17)中，u(h)和u(z)表示建筑物顶和z高度的边界层风速。
[0090]
第二种背风区插值方案
[0091]
将尾流边界分为垂直区域和水平区域。
[0092]
垂直区域指原尾流区的正上方，位置坐标(x,y,z)处的风速插值方案如下。
[0093]
u(x,y,z)＝(u
ref_vert-u(z))exp(-ξ
1.5
) u(z)(18)
[0094][0095][0096][0097]
(18)中u(z)表示z高度边界层风速，ξ是一个无量纲参数。(21)中ht_ellipse(x,y)表示当前位置的椭球高度，本技术中ht_ellipse(x,y)＝z_vert，因此(21)可以转化为
[0098][0099]
水平区域指原尾流区的两侧，位置坐标(x,y,z)处的风速插值方案如下。
[0100]
u(x,y,z)＝(u
ref_hor-u(z))exp(-η2) u(z)(23)
[0101][0102][0103][0104]
其中，u(z)表示z高度边界层风速，η是一个无量纲参数。
[0105]
s323建筑物顶部涡旋区
[0106]
本研究提供了两种顶部涡旋区插值方案以供选择，通过风洞实验资料与模拟结果的对比试验，挑选合适方案。
[0107]
第一种顶部涡旋区插值方案
[0108]
将建筑物顶部的插值区域看作椭圆柱体，长度和高度由建筑物的几何性质确定，
尺度用r来表示。
[0109][0110]hcm
(heightofvortex)＝0.22*r(28)
[0111]
lc(lengthofvortex)＝0.9*r(29)
[0112]bs
＝smallerofupwindbuildingheightorwidth(30)
[0113]bl
＝largerofupwindbuildingheightorwidth(31)
[0114]
如图2所示，建筑物顶部涡旋区可分为两个区域，区域1也被成为逆流区，采用调整的负对数风速插值方案，区域2采用正对数风速插值方案。
[0115]
区域1
[0116]
u(z)＝-u(z)*c1(32)
[0117][0118]
区域2
[0119]
u(z)＝u(z)(34)
[0120]
(27)-(34)中h是建筑物高度，h
cm
是插值区的高度，u(z)表示z高度边界层风速。
[0121]
第二种顶部涡旋区插值方案
[0122]
如图3所示，逆流区的风廓线形状应为倒楔形。回流区内采用高度线性垂直变化的插值公式。
[0123][0124]z′
＝h h
cm-z(36)
[0125]
其中，h是建筑物高度，h
cm
是插值区(回流区)的高度，u(z)表示z高度的来流边界层风速。
[0126]
注意的是需要考虑夹角入流情形。如图4所示，当入流风和建筑物的夹角θ过大时(一般认为θ》20
°
)，即当θ》20
°
时，迎风面的入流风速分解为平行迎风面分量u0(z)
//
和垂直迎风面分量u0(z)
⊥
，u0(z)
⊥
代替入流风速根据s321、s322和s323有关参数化方案进行插值，u0(z)
//
保持不变。
[0127]
如图5所示，背部空腔区尾流区插值时，用建筑等效长度l
eff
和等效宽度w
eff
替代原有长宽。
[0128]weff
＝wsinθ lcosθ(37)
[0129]
l
eff
＝wcosθ lsinθ(38)
[0130]
综上，当入流风和建筑物的夹角过大时，建筑物前后的插值方案如表1、表2所示。
[0131]
表1夹角入流建筑物位移区插值
[0132][0133]
表2夹角入流建筑物空腔区、尾流区插值
[0134][0135]
s330对建筑街渠风场参数化
[0136]
当s/h》2.5或1.4《s/h《2.4时，采用单体建筑物周围风场参数化方案；
[0137]
当s/h《1.4时，采用sr参数化方案或mr参数化方案。
[0138]
具体而言，sr参数化方案(标准rockle插值方案，standardrocklestreetcanyonmodel)采用空间距离参数s
*
和比例形态参数w/h来确定街渠流场类型，s
*
的计算公式为
[0139][0140]
当s《s
*
，建筑街渠内形成一个顺时针涡旋，水平方向和垂直方向的速度插值公式为
[0141][0142][0143]
上式中s为两个建筑物之间距离，d为位置坐标到上风建筑物墙面的距离，u(h)表示上风向建筑物顶的边界层风速。
[0144]
mr参数化方案(改进的rockle插值方案，modifiedrocklestreetcanyon
model)重新计算空间距离参数s
*
。
[0145][0146]
并将街渠流场分为三个区域：(1)中心区域，风速插值方案和传统街渠插值方案相同；(2)垂直湍流扩算区域，需考虑街渠上方的动量输入和输出；(3)水平湍流扩散区域，需考虑区域(1)水平两侧的动量输入和输出。区域(2)(3)的空间形态是由入流方向建筑物延伸出来的三角柱体或楔形体。区域(1)(2)(3)的边长分别为
[0147]
δ
sc
(x
can
)＝w-2δ
lw
(x
can
)(43)
[0148]
δ
vw
(x
can
)＝0.2x
can
(44)
[0149]
δ
lw
(x
can
)＝0.2x
can
(45)
[0150]
区域(1)的风速参数化方案为
[0151][0152][0153]
区域(2)的风速参数化方案为
[0154][0155]
区域(3)的风速参数化方案为
[0156][0157]ylw
＝[|y
can
|-δ
sc
(x
can
)/2](50)
[0158]
上式中，x
can
表示位置坐标到上风建筑物墙面的距离，y
can
表示位置坐标到上风建筑物墙面中心平面的距离，p为经验参数，用于调整风洞实验结果与模拟结果的差异，此处取0.25，u
roof
表示建筑物顶风速，u
bl
(z)表示z高度下边界层风速，γ为经验衰减系数，可根据观测资料或者风洞实验结果决定，此处取0.3。
[0159]
s400根据质量守恒方程调整所述初始风场，得到最终风场。
[0160]
由于上文所用风场参数化方案不能反映建筑物墙壁对气流的作用，因此在插值区域(特别是建筑物周围)不满足质量守恒连续方程。质量守恒是指速度的散度为0，假设插值风场为u0(u0,v0,w0)，有
[0161][0162]
本研究采用变分方法求解满足(51)的最终风场u(u,v,w)。首先，以u(u,v,w)和u0(u0,v0,w0)构建的变分函数公式
[0163][0164]
(52)中，α1和α2定义为高斯精度模，用于控制变分公式对水平 u、v和垂直w风速分量的调整力度。α1和α2的值越大，说明变分公式对速度的调整幅度较小，即相对于初始风速的改变较小。
[0165]
引入拉格朗日乘数λ，联立(51)(52)可得
[0166][0167]
根据欧拉-拉格朗日方程组，u(u,v,w)和u0(u0,v0,w0)满足以下条件时，e(u,v,w,λ)取极小值。
[0168][0169][0170][0171]
将(54)代入(51)，得到关于λ的泊松方程。
[0172][0173]
α1和α2取常数，α1＝α2＝0.5。根据分辨率dx、dy、dz，以及计算区域长宽高的格点数nx、ny、nz，生成nx*ny*nz大小的三维网格。计算效率采用1/2网格，λ
i,j,k
表示(i,j,k)网格的值，速度则表示为网格边界值(图6)。对(55)按照图6网格离散化，对于非边界格点，有
[0174][0175]
等式两边简化后，有
[0176][0177]
令并定义边界系数为e，f，g，h， m，n，o，p，q。
(57)可简化为
[0178][0179]
整理后得到
[0180][0181]
其中
[0182][0183]
采用超松弛算法(sor)，根据(60)，有
[0184][0185]
ω为sor收敛速度参数，本技术ω取1.78。
[0186]
当时，根据(62)求得λ的最优解。ε是一个小量，可按照计算速度要求进行调整，为了提高计算效率，加快收敛速度，本技术ε取0.01。
[0187]
对于边界格点，离散化时需要考虑边界条件。本研究对于底面及建筑物采用第二类边界条件其他边界采用第一类自由边界条件λ＝0。
[0188]
比如，若正i方向存在建筑物，即建筑物墙壁位于i 1/2处，有
[0189][0190]
根据第二类边界条件，此时
[0191][0192]
代入(55)并按照图6网格离散化，有
[0193]
[0194]
简化后得
[0195][0196]
此时，o＝0.5，p＝1，q＝1，e＝0，f＝1，g＝1，h＝1，m＝1，n＝1。同理可得其他边界条件系数。
[0197]
边界条件系数在质量守恒方程中，代表当前格点与建筑物墙面边界的相对位置。本研究将相对位置分为14类，详见表3。
[0198]
表3边界条件系数表
[0199][0200]
本实施例中，具体对国内城市某地部分街区进行了模拟，建筑物信息来自地理信息系统(gis)软件，模拟后的建筑物分布图如图7。通过以上方法计算出对应高度下水平面的格点风场，k表示输出的高度为2/4/6/8/10/20/40/.../140m，通过matlab画图程序 read_wind.m，得到图8。
[0201]
本实施例中气象数据来自欧洲数值预报中心(ecwmf)，每日08 时、20时预报结果，每一次预报53个时次，预报时长为10天，在 0-72hr内每3hr一个预报结果，72hr-240hr每6hr一个预报结果。
[0202]
将本方法模拟结果与风洞试验结果进行比较，这里用的风洞数据来自于cedval数
据集。通过风洞观测和本方法模拟结构比较，如图9所示，本方法较准确地模拟出街区峡谷风场的关键区域。选取五个不同垂直剖面(x/h＝-2.58，x/h＝-0.5，x/h＝0，x/h＝0.5和 x/h＝2.58)的风洞观测和方法模拟结果，比较全风速u u分量、w分量的垂直变化廓线(图10和图11)。从图中看出，廓线街渠中心区域(x/h＝0.5)变化幅度较大，建筑物高度以下的模式模拟结果存在一定高估，u和w分量的模拟误差较大。分析原因可能是由于街渠中心区域受到气流汇流作用，风速迅速增大，加之墙体对风速的阻挡和拖曳使得此区域风速梯度很大，增加了模拟难度，从而导致模拟风速偏大较多。x/h＝-2.58，x/h＝-0.5， x/h＝0.5和x/h＝2.58四个垂直平面u的模拟效果较好。
[0203]
综上所述，本发明通过利用简单的气象输入以及基本的建筑物信息能够快速模拟出超大城市中超高层建筑物对周围风场环境的影响，计算量小，同时加入坡度因素，模拟结果更加准确。
[0204]
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。