一种基于CFD与GIS的超大城市风环境模拟方法

一种基于cfd与gis的超大城市风环境模拟方法
技术领域
1.本发明属于城市风环境模拟领域，具体涉及一种基于cfd与gis的超大城市风环境模拟方法。


背景技术：

2.城市风环境是近地面自然风在城市下垫面影响下形成的风速与风向的分布情况，城市风环境影响人们的体感舒适度、城市空气质量等。随着城市规模的扩大，尤其对于城区常住人口在1000万以上的超大城市，密度与高度较大的建筑群对周围的风速产生了显著的削弱作用，不仅容易造成严重的城市风环境问题，还可能降低城市内居民的热舒适度，并可能使大气污染物在城区聚集，从而降低空气质量，因此，需要对超大城市内的风速与风向进行高空间分辨率数值模拟从而为改善城市风环境提供科学依据。
3.现有技术中，通常采用实地观测与数值模拟相结合的方法来获取城市风场。实地观测方法由于城市中的观测站点分布稀疏，且需要耗费大量人力和物力进行风速、风向测量，虽然能获取时间连续的离散观测值，但是很难获得区域或者空间连续的风场。中尺度模式可以模拟城市风场，但也受限于模式机理无法获取百米尺度以下的城市风场，无法满足获取几十米空间分辨率的风场来探究建筑三维空间形态对城市风环境影响的需求。此外，城市中大量三维空间形态复杂的建筑对风环境模拟带来较大的计算成本。


技术实现要素：

4.鉴于现有技术中的上述缺陷或不足，本发明旨在提供一种基于cfd与gis的超大城市风环境模拟方法，用以实现超大城市几十米空间分辨率风环境的数值模拟。为了实现上述目的，本发明实施例采用如下技术方案：
5.一种基于cfd与gis的超大城市风环境模拟方法，包括如下步骤：
6.步骤s1，获取待研究城市的建筑物数据，并采用地理信息系统gis软件中的矢栅转换工具将建筑物数据中的建筑矢量多边形转换为建筑立面栅格与建筑高度栅格；
7.步骤s2，通过gis软件中的栅格运算工具与重分类工具将建筑立面栅格与建筑高度栅格转换为简化建筑密度栅格与简化建筑容积高度栅格，并进一步将简化建筑密度栅格与简化建筑容积高度栅格转换为建筑块栅格；
8.步骤s3，在cad软件中将建筑块栅格转化为三维建筑块实体，以刻画建筑物的三维空间形态；
9.步骤s4，在计算流体力学cfd软件中建立cfd数值模拟的计算域，在所述计算域中使用所述三维建筑块实体刻画建筑物的三维空间形态，并使用非结构化网格对计算域进行离散化；
10.步骤s5，在cfd软件中设置湍流模型参数、边界条件与求解器，对计算域进行初始化并完成数值求解，模拟当前城市风场。
11.上述方案中，步骤s1的建筑立面栅格中，使用二元值a与b表示栅格对应的像元是
否被建筑覆盖；建筑高度栅格中，使用建筑高度属性表示建筑的高度，所述建筑高度属性使用大于等于0的实数表示。
12.上述方案中，步骤s2包括：
13.步骤s21，建立覆盖建筑立面栅格或建筑高度栅格的连续格网，该格网由形状大小一致的单元构成，单元的边长r由建筑密度与城市道路宽度决定；采用栅格运算工具对单元内的建筑密度与建筑容积高度进行分区统计，生成建筑密度栅格与建筑容积高度栅格；
14.步骤s22，使用gis栅格运算工具对建筑密度栅格进行重分类，使用二元值c与d分别表示不含建筑覆盖的区域与含建筑物覆盖的区域，得到简化建筑密度栅格；采用相等间隔分类方法对建筑容积高度栅格进行重分类得到简化建筑容积高度栅格；使用栅格计算工具将简化建筑密度栅格和简化建筑容积高度栅格转换为建筑块栅格。
15.上述方案中，所述采用栅格运算工具对单元内的建筑密度与建筑容积高度进行分区统计，采用建筑覆盖率c与建筑容积高度hvol来表示单元内的建筑密度与建筑容积高度，计算方法如公式(1)和(2)所示：
[0016][0017][0018]
式(1)和(2)中，s
c，i
表示第i个边长为r的单元的面积，sj表示s
c，i
范围内第j个建筑要素的立面面积，hj表示sj的高度属性。
[0019]
上述方案中，步骤s3包括：
[0020]
步骤s31，使用gis软件中的栅矢转换工具将建筑块栅格转换为建筑块矢量，栅格属性值被转换为建筑块的高度属性；所述的建筑块栅格使用简化后的建筑容积高度值表示建筑块的高度，使用像元值0表示除建筑块以外的其他区域；
[0021]
步骤s32，使用cad软件将建筑块矢量转换为三维建筑块实体；所述转换方法为：使用gis软件将建筑块矢量导出为建筑块三维多段线，并使用cad软件中的拉伸工具将所述建筑块三维多段线转换为建筑块三维实体，拉伸高度为建筑块矢量的高度属性值，从而将建筑块矢量转换为三维建筑块实体。
[0022]
上述方案中，步骤s4包括：
[0023]
步骤s41，依据待研究城市与三维建筑块实体建立cfd数值模拟的计算域，建立计算域的方法为：在cad软件中建立计算域三维实体；所述计算域的边界包括底边界、侧边界、顶边界与建筑表面；使用布尔运算工具将计算域三维实体与三维建筑块实体求减运算；
[0024]
步骤s42，对计算域进行非结构化网格划分，包括面网格尺寸与体网格尺寸的设置；将建筑表面网格的最大尺寸设置为单元边长r，将最小面网格的尺寸限制在2/3r；底边界表面网格的最大尺寸设置为100m以下，底边界、侧边界、顶边界与建筑表面的面网格增长率为1.05以下，以保证核心计算域至侧边界网格尺寸的均匀过度。
[0025]
上述方案中，所述计算域的范围包括水平范围和垂直范围，水平范围h
hor
与垂直范围h
ver
分别由待研究城市的水平范围h
city
与三维建筑块实体高度的最大值h
max
决定，且h
hor
应为h
city
的五倍以上，h
ver
至少达到待研究城市的边界层高度；所述核心计算域为水平范围为h
city
的cfd数值模拟核心区域。
[0026]
上述方案中，步骤s5包括：
[0027]
步骤s51，在cfd软件中设置湍流模型及相关模拟参数；
[0028]
步骤s52，设置cfd数值模拟的边界条件；
[0029]
步骤s53，设置cfd数值模拟的求解器与cfd数值模拟的初始化方法，并对计算域进行初始化与迭代求解。
[0030]
上述方案中，所述湍流模型为基于雷诺平均方法的湍流模型。
[0031]
上述方案中，所述边界条件包括：顶边界为对称边界或无滑移壁面，底边界与建筑表面为无滑移壁面，出口边界为自由流，侧边界为速度入口，速度入口采用指数风剖面定义风速，计算方法如公式(3)所示：
[0032][0033]
式(3)中，z0为参考风速v0的高度，α为与地貌类别与海拔高度相关的地面风压高度变化系数；风速入口还使用x，y，z三个方向上的分量u，v，w定义入口风向。
[0034]
本发明实施例所提供的技术方案具有如下有益效果：
[0035]
本发明实施例所提供的基于cfd与gis的超大城市风环境模拟方法，实现了几十公里城市范围内百米以下尺度的城市风场模拟，能够以几十米的空间分辨率刻画三维建筑物与周围风速的相互关系，为研究城市建筑三维空间形态对高密度城市风环境的影响、规划高密度城市通风廊道进而改善城市通风状况、提高行人热舒适度、城市污染源扩散等提供技术支撑，从而缓解城市热岛效应。
[0036]
当然，实施本发明的任一产品或方法并不一定需要同时达到以上所述的所有优点。
附图说明
[0037]
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来说，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。
[0038]
图1是本发明实施例所述基于cfd与gis的超大城市风环境模拟方法流程图；
[0039]
图2是本发明一个具体实施例中上海市中心城区的建筑局部三维场景示意图；
[0040]
图3是本发明一个具体实施例中上海市中心城区的简化后建筑块局部三维场景示意图；
[0041]
图4是本发明一个具体实施例中上海市中心城区的计算域示意图；
[0042]
图5是本发明一个具体实施例中上海市中心城区非结构化网格划分局部示意图；
[0043]
图6是本发明一个具体实施例中上海市中心城区的近地面10m剖面风速云图。
具体实施方式
[0044]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。
需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征也可以相互组合。
[0045]
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。在本发明的描述中，术语“第一”、“第二”、“第三”、“第四”等仅用于区分描述，而不能理解为只是或暗示相对重要性。
[0046]
为解决风环境评估中计算流体力学(cfd)数值模拟范围难以覆盖整个城市、对建筑影响下的风环境进行数值模拟时建筑物数量较大的问题，本发明提出了一种基于cfd与gis的超大城市风环境模拟方法，通过gis手段有效简化建筑要素，实现对超大城市几十公里范围、几十米空间分辨率风环境的数值模拟，模拟过程覆盖整个城市范围，提供的基于cfd与gis的超大城市风环境模拟方法，实现了几十公里城市范围内几十米空间分辨率的城市风场模拟，为研究城市建筑三维空间形态对高密度城市风环境的影响提供技术支撑。
[0047]
参见图1，本发明实施例所提供的基于cfd与gis的超大城市风环境模拟方法，包括如下步骤：
[0048]
步骤s1，获取待研究城市建筑物数据，并采用gis软件中的矢栅转换工具将建筑物数据中的建筑矢量多边形转换为建筑立面栅格与建筑高度栅格。
[0049]
本步骤中，所述失栅转换的目的是：栅格数据具有高效的空间叠置运算能力，使用栅格数据能够提高建筑块的生成效率，使大量建筑物高效简化为三维实体。所述建筑数据中包括具有建筑高度属性的建筑矢量多边形数据。所述建筑要素，包括地理坐标、形态与高度。
[0050]
本步骤中，所述待研究城市为常住人口数量大于1000万的超大城市。
[0051]
优选地，所述gis软件为商用gis软件arcgis 10.2系列中的arcmap，所述的矢栅转换工具为arcmap中的“要素转栅格”。
[0052]
本步骤中，所述建筑立面栅格中，使用二元值a与b表示栅格对应的像元是否被建筑覆盖；所述建筑高度栅格中，使用建筑高度属性表示建筑的高度，所述建筑高度属性使用大于等于0的实数表示。所述二元值a与b可以选择0和1或其他相区别的两个不同实数。优选地，建筑立面栅格使用像元值0与1表示建筑物立面的水平位置，像元值1表示该像元上被建筑覆盖，像元值0则表示该像元上没有建筑覆盖。
[0053]
步骤s2，通过gis软件中的栅格运算工具与重分类工具将建筑立面栅格与建筑高度栅格转换为简化建筑密度栅格与简化建筑容积高度栅格，并进一步将简化建筑密度栅格与简化建筑容积高度栅格转换为建筑块栅格。
[0054]
本步骤具体包括：
[0055]
步骤s21，建立覆盖建筑立面栅格或建筑高度栅格的连续格网，该格网由形状大小一致的单元构成，单元的边长r由建筑密度与城市道路宽度决定；采用栅格运算工具对单元内的建筑密度与建筑容积高度进行分区统计，生成建筑密度栅格与建筑容积高度栅格。本步骤中，所述单元为形状规则的多边形，其边长r能够反映建筑的空间分辨率。
[0056]
优选地，所述gis栅格运算工具为arcmap中的“分区统计”；
[0057]
优选地，采用建筑覆盖率c与建筑容积高度h
vol
来刻画单元内的建筑密度与建筑容积高度，其计算方法如公式(1)和(2)所示：
[0058][0059][0060]
式(1)和(2)中，s
c，i
表示第i个边长为r的单元的面积，sj表示s
c，i
范围内第j个建筑要素的立面面积，hj表示sj的高度属性。
[0061]
步骤s22，使用gis栅格运算工具对建筑密度栅格进行重分类，使用二元值c与d分别表示不含建筑覆盖的区域与含建筑物覆盖的区域，得到简化后的建筑密度栅格；采用相等间隔分类方法对建筑容积高度栅格进行重分类得到简化建筑容积高度栅格；使用栅格计算工具将简化后的建筑密度栅格和简化建筑容积高度栅格转换为建筑块栅格。所述二元值c与d可以选择0和1或其他两个相区别的不同实数。
[0062]
优选地，所述gis栅格运算工具为arcmap中的“重分类”；
[0063]
优选地，将建筑密度栅格重分类为值为1的建筑块所覆盖区域与值为0的其他区域，分类结果为简化的建筑密度栅格；
[0064]
优选地，对建筑容积高度栅格采用对建筑高度属性的相等间隔分类方法将建筑容积高度分为多级，生成简化的建筑容积高度栅格。重分类后的值为每一个分级中的中位数。所述相等间隔分类方法能够保证建筑块在垂直方向上的空间分辨率相对均匀。所述分类级数与建筑高度的最大、最小值与待研究城市建筑块在水平方向上的空间分辨率相关。
[0065]
步骤s3，在cad软件中将建筑块栅格转化为三维建筑块实体，以刻画建筑的三维空间形态。
[0066]
本步骤具体包括：
[0067]
优选地，所述cad软件为商用软件公司autodesk出品的autocad 2018；
[0068]
步骤s31，使用gis软件中的栅矢转换工具将建筑块栅格转换为建筑块矢量，栅格属性值被转换为建筑块的高度属性；所述的建筑块栅格使用简化后的建筑容积高度值表示建筑块的高度，使用像元值0表示除建筑块以外的其他区域。
[0069]
优选地，所述栅矢转换工具为arcmap中的“栅格转面”；
[0070]
所述建筑块栅格使用简化后的建筑容积高度值表示建筑块的高度，使用像元值0表示除建筑块以外的其他区域。
[0071]
步骤s32，使用cad软件将建筑块矢量转换为三维建筑块实体，所述转换方法为：使用gis软件将建筑块矢量导出为建筑块三维多段线，并使用cad软件拉伸工具将之拉伸为三维实体，拉伸高度为建筑块矢量的高度属性值。
[0072]
步骤s4，在计算流体力学cfd软件中建立cfd数值模拟的计算域，在所述计算域中使用所述三维建筑块实体刻画建筑物的三维空间形态，并使用非结构化网格对计算域进行离散化。
[0073]
本步骤具体包括：
[0074]
步骤s41，依据待研究城市与三维建筑块实体建立cfd数值模拟的计算域，建立计算域的方法为：在cad软件中建立计算域三维实体；使用布尔运算工具将计算域三维实体与三维建筑块实体求减运算。
[0075]
所述计算域的范围包括水平范围和垂直范围，水平范围h
hor
与垂直范围h
ver
分别由
待研究城市的水平范围h
city
与建筑块高度属性的最大值h
max
决定。h
hor
应为h
city
的五倍以上，hver应至少达到待研究城市的边界层高度。计算域的边界分为底边界、侧边界、顶边界与建筑表面；所述核心计算域为水平范围为h
city
的cfd数值模拟核心区域。
[0076]
步骤s42，对计算域进行非结构化网格划分。本步骤中，网格类型为非结构化网格，能够通过设置网格参数的手段有效捕捉城市内气流的复杂流动，与使用结构化网格相比能有效减少网格数量以降低cfd数值模拟的运算成本。
[0077]
优选地，所述非结构化网格的参数包括面网格尺寸与体网格尺寸的设置。建筑表面网格的最大尺寸设置为单元边长r，最小面网格的尺寸限制在2/3r；底边界表面网格的最大尺寸设置为100m以下，建筑表面、底边界、侧边界、顶边界的面网格增长率为1.05以下，以保证核心计算域至侧边界网格尺寸的均匀过度。
[0078]
步骤s5，在cfd软件中设置湍流模型参数、边界条件与求解器，对计算域进行初始化并完成数值求解，模拟当前城市风场。
[0079]
本步骤具体包括：
[0080]
步骤s51，在cfd软件中设置湍流模型及相关模拟参数。
[0081]
优选地，本步骤中湍流模型为基于雷诺平均方法的rng k-epsilon湍流模型。
[0082]
步骤s52，设置cfd数值模拟的边界条件。
[0083]
优选地，本步骤中边界条件包括：顶边界为对称边界或无滑移壁面；底边界与建筑表面为无滑移壁面；侧边界为速度入口，速度入口采用指数风剖面来定义风速，计算方法如公式(3)所示：
[0084][0085]
式(3)中，z0为参考风速v0的高度，α为与地貌类别与海拔高度相关的地面风压高度变化系数；侧边界即风速入口还使用x，y，z三个方向上的对应分量u，v，w来定义入口风向。
[0086]
步骤s53，设置cfd数值模拟的求解器与cfd数值模拟的初始化方法，并对计算域进行初始化与迭代求解。
[0087]
优选地，本步骤中所述求解器为压力耦合式求解器，初始化方法为混合初始化与多重网格初始化。
[0088]
下面通过一个具体的实施例，对本发明作进一步详细说明。本实施例以上海市中心城区为例，进行风场模拟。
[0089]
上海市是典型的超大城市。上海市政府将上海市外环线s20以内区域定义为上海市中心城区，七普所调查的各区常住人口中，城镇人口约为2220.94万人，约占总人口数的89.3％；中心城区人口约占总人口数的26.9％，符合超大城市的指标。本实施例以上海市中心城区作为待研究城市。
[0090]
采用本发明的基于cfd与gis的超大城市风环境模拟方法进行风场模拟，包括以下步骤：
[0091]
步骤s1，获取待上海市中心城区的建筑物数据，并在arcmap中使用矢栅转换工具将该数据转化为建筑立面栅格与建筑高度栅格。
[0092]
本步骤通过权衡栅格像元大小与栅格文件大小，将建筑立面栅格与建筑高度栅格的像元大小设置为1m。
[0093]
可选地，当待研究城市建筑矢量的层数已知而层高未知时，建筑数据的高度属性可采用“建筑层高
×
3m”表示，处理后得到的上海市中心区的建筑局部三维场景如图2所示。
[0094]
步骤s2、通过gis软件的栅格运算与重分类工具将建筑立面栅格与建筑高度栅格转换为简化建筑密度栅格与简化建筑容积高度栅格，并将之转换为建筑块栅格。采用栅格运算工具arcmap中的“分区统计”，计算建筑密度与建筑容积高度；采用栅格运算工具arcmap中的“重分类”，生成简化建筑密度栅格与简化建筑容积高度栅格；采用栅格运算工具arcmap中的“栅格计算器”生成建筑块栅格。
[0095]
优选地，采用arcmap中的“泰森多边形”在城市地图中建立连续格网。所述的单元区域为所建立的单位正六边形格网，使用正六边形格网的优越性在于：该格网中每个正六边形要素与相邻要素有且只有一条长为r的邻边，在三维实体中则是一个面，能够避免后续建立计算域时出现布尔减运算错误。
[0096]
优选地，当选取r为30m时，将上海市中心城区平均宽度为30m的一级、二级道路与建筑物分离，也能够将建筑物与河流、大片公共空间分离。
[0097]
优选地，建筑密度栅格重分类阈值为0.4，建筑容积高度栅格被重分类为相等间隔的71级，每级的间隔为5m；得到的建筑块局部三维场景如图3所示。
[0098]
步骤s3，在cad软件中将建筑块转化为三维建筑实体，以几十米的空间分辨率刻画建筑的三维空间形态；建筑块栅格由简化后的建筑密度栅格与简化后的建筑高度栅格相乘得到；将arcmap导出的建筑块三维多段线转换为三维建筑块实体的具体步骤为：
[0099]
步骤s31，使用gis软件将建筑块矢量导出为建筑块三维多段线；
[0100]
步骤s32，使用autocad 2018软件拉伸工具将建筑块三维多段线转换为三维实体，向z轴负方向拉伸高度为建筑块矢量的高度属性值。
[0101]
步骤s4，建立cfd模拟的计算域，如图4所示，并在cfd软件中使用非结构化网格对计算域进行离散化。优选地，所述cfd软件采用西门子公司的cfd软件starccm 2020。
[0102]
优选地，计算域为一个底面直径h
hor
约60km，高度h
ver
为1km的圆柱形实体，便于定义入口风向；该设计不仅权衡了网格数量与计算域垂直高度对流场产生的影响，还在侧边界与核心计算域之间预留了较大的湍流发展距离使流场风速达到稳定，核心计算域的水平范围h
city
覆盖了直径约1000km上海市中心城区，是模拟上海市中心城区风场的核心区域。
[0103]
优选地，所述非结构化网格的参数包含：将建筑表面的表面网格最大尺寸设置为3.1节中的r，即30m，将最小面网格的尺寸限制在20m；底边界表面网格的最大尺寸设置为100m，建筑表面、底边界的表面网格增长率为1.025，以保证核心计算域至侧边界网格尺寸的均匀过度；侧边界、顶边界的面网格增长率为1.05。在此基础上采用四面体体网格工具生成体网格，并对底边界与建筑表面附近的区域划分棱柱层网格以较好捕捉近地面的复杂湍流。四面体体网格的增长率为1.15，在底边界与建筑表面上设置的棱柱层数分别为8层与5层，棱柱层总厚度均为30m，相邻棱柱层网格的增长率均为1.06。计算域非结构化网格划分局部如图5所示：建筑表面被整齐的三角形网格所填充，非建筑区域的面网格较光顺，尺寸具有显著差异，降低了数值模拟的计算成本。
[0104]
步骤s5，在cfd软件中设置湍流模型参数、模拟的边界条件与求解器，对计算域进行初始化并完成数值求解，模拟当前城市风场。
[0105]
本步骤中，cfd软件为商用软件ansys fluent19.0。
[0106]
优选地，所述边界条件包括侧边界、顶边界、底边界与建筑表面。侧边界为速度入口，需要设置入口风速与入口风向。由于使用气象站10m处的风速值计算边界条件，参考高度z0为10m，α为与地貌类别与海拔高度相关的地面风压高度变化系数，α为0.3。
[0107]
优选地，据上海市宝山、浦东、虹口气象观测站2018-2020年的全年风速观测数据，设置v0＝3m/s，风向为东南进行了模拟，在连续运算了6小时、迭代了800步以后残差下降至10-4
，近地面10m处水平风速剖面如图6所示。该城市的中心区位于上海市内环以内，建筑的三维空间形态复杂，对周围的风速具有显著的削弱作用，使中心区的风速显著小于中心区以外的风速。
[0108]
由以上技术方案可以看出，本发明实施例所提供的基于cfd与gis的超大城市风环境模拟方法，实现了几十公里城市范围内复杂建筑影响下的城市风场模拟，接近超大城市风场的真实环境，不仅准确表达了风速与建筑物间的相互作用，同时模拟周期短，计算成本低，有效地节约了模拟资源。
[0109]
以上描述仅为本发明的较佳实施例以及对所运用技术原理的说明，并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的优选实施例。本领域技术人员应当理解，本发明中所涉及的发明范围，并不限于上述技术特征的特定组合而成的技术方案，同时也应涵盖在不脱离所述发明构思的情况下，由上述技术特征或其等同特征进行任意组合而形成的其它技术方案。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。