一种结合实测的风致雪漂移的CFD数值模拟方法

一种结合实测的风致雪漂移的cfd数值模拟方法
技术领域
1.本发明属于土木工程领域，涉及风致雪漂移的cfd数值模拟方法。


背景技术：

2.近些年来，在寒冷地区暴雪给人们的出行，通讯及生产生活带来了极大的不便，同时也极大地威胁了房屋以及公共建筑的安全，其中风致雪漂移容易导致建筑局部位置的荷载过大，而这种建筑上的不均匀雪荷载往往会对结构产生额外的作用，使得结构更容易失稳。因而对风雪的研究以及预防变得越来越重要。
3.在风雪灾害的早期研究中，囿于计算机的条件，国内外的学者们以实测以及实验研究为主，随着社会的不断进步和科学技术的迅猛发展，以计算流体力学为基础的数值模拟研究逐渐被大众所采纳。目前的风致雪漂移数值模拟研究采用单向耦合的假设，参数设置采用经验值，单向耦合即为假设不考虑雪相对空气相的影响，将雪相视为被动运输的固体介质，而雪相的沉积与侵蚀在实际中会与对积雪边界条件产生影响，从而进一步影响流场分布，而该方法本质上为稳定的风场下对雪相的扩散进行计算，这与实际情况大相径庭，降低了模拟的准确性。


技术实现要素：

4.本发明是要解决现有的数值模拟方法不考虑风雪之间的相对运动、准确性差的技术问题，而提供一种结合实测的风致雪漂移的cfd数值模拟方法。该方法考虑风雪相对运动，从风雪两相流的建模，到参数设置，模拟计算等一系列过程进行了优化，得到风致雪漂移的cfd模拟数值。
5.本发明的结合实测的风致雪漂移的cfd数值模拟方法，按以下步骤进行：
6.步骤一：收集数值模拟需要的基本信息：收集具体研究的建筑物的几何特征、该建筑所处的地理位置、地貌参数、参考高度h0处平均风速v0、风剖面规范、湍动能规范、湍动强度规范、积分湍流尺度规范、湍能耗散率规范；同时测得建筑物所在地积雪的雪颗粒密度、雪颗粒直径ds、休止角、风速、风向及建筑物的初始积雪深度，计算阈值摩擦速度、跃移雪粒速度、地面粗糙度、建筑物墙面粗糙度、雪面粗糙度和雪粒沉降速度wf；
7.步骤二：建立数值模拟需要的网格模型：利用步骤一所收集的建筑物几何信息和初始积雪分布形态进行模型的建立以及网格的划分，得到网格文件；
8.步骤三：编写数值模拟用的空气项边界条件：利用步骤一所收集的建筑物所在地的地理位置信息、参考高度平均风速、风剖面、湍动能、湍动强度、积分湍流尺度、湍流耗散率规范信息编写空气相的边界条件，再在visual studio中以c语言、fluent二次开发接口为基础，编写指数型风剖面、湍动能、湍流耗散率程序，保存c语言文件，并将其与步骤二网格模型文件放在同一文件夹中；
9.步骤四：编写数值模拟用的雪相边界条件：利用步骤一所收集的雪颗粒密度、建筑物所在地的重力加速度、阈值摩擦速度，采用雪浓度经验公式及其公式中所采用的雪相运
动速度、雪颗粒摩擦速度、跃移和悬移层临界高度信息，在visual studio中以c语言、fluent二次开发接口为基础编写雪相体积分数公式，作为雪项的入口边界条件，保存c语言文件，并将其与步骤二模型网格文件、步骤三的c语言文件放在同一文件夹中；
10.步骤五：数值模拟需要的程序接入，具体操作如下：打开ansys-fluent界面选择2d/3d、双精度、并行，在general options-working directory下选择包含步骤二的网格文件、步骤三的c语言文件、步骤四的c语言文件的文件夹目录，选择环境下set up compilation environment for udf接口打开软件；选择ansys-fluent主菜单下的user defined
–
functions
–
compiled，选择步骤三、四的c语言文件进行加载编译，得到空气相及雪相的入口速度及体积分数；
11.步骤六：采用ansys-fluent软件，先设置模型、材料、求解器参数，再导入步骤二所得的网格模型、步骤五所得的空气相及雪相的入口速度及体积分数，进行风雪两相流模拟，得到壁面剪切力、空气相及雪相的速度及体积分数的收敛文件；
12.步骤七：利用步骤六计算出的收敛文件，进行壁面剪切力和速度的后处理：将ansys-fluent计算出的收敛文件以tecplot支持格式保存，用tecplot进行后处理计算，通过提取计算的计算域三维空间中x、y、z方向的速度、压力、壁面剪切力的分量，用tecplot自定义函数specify equations，做出计算域范围内的压力、速度、壁面剪切力的云图、等值线图；
13.步骤八：数值模拟雪沉积侵蚀的后处理计算：基于步骤七计算获得的壁面剪切力，通过tecplot提取模型底面或底面横、纵方向中线上的壁面剪切力，通过公式把壁面剪切力转换为壁面的摩擦速度，其中τ为壁面剪切力，ρ为空气密度；对壁面摩擦速度与阈值磨擦速度进行比较，当壁面摩擦速度大于阈值磨擦速度时，积雪发生侵蚀；当壁面摩擦速度小于阈值磨擦速度时，积雪发生沉积；分别计算侵蚀量、沉积量，并且计算积雪改变厚度，最后将积雪改变厚度与建筑物的初始积雪深度进行叠加，得到积雪的最终分布形态，完成结合实测的风致雪漂移的cfd数值模拟。
14.更进一步地，步骤二中所述的网格文件的求解过程是：在ansys-icem中按照建筑物几何特征确定模型计算域且计算域设置满足阻塞率小于3％，接着在icem中依据模型建立网格模型，为保证计算顺利收敛及模拟的准确性，依据模型大小设置合理的网格尺寸，调控网格偏斜率(skewness)三角形与四面体不大于0.95，调控网格纵横比(aspect ratio)，调控网格压扁程度(squish)在0-1之间，最后以ansys-fluent格式输出求解的网格文件。
15.更进一步地，其中合理的网格尺寸是指：网格模型控制结构化网格质量在0.6～1，非结构化网格质量在0.3～1；网格纵横比为：流动核心区内网格纵横比小于5:1，边界层网格纵横比小于10:1。
16.更进一步地，步骤四中所述的雪浓度经验公式为1990版的pomeroy&gray、1992版的pomeroy&male所提出的雪浓度经验公式。
17.更进一步地，步骤六中所述的风雪两相流模拟，具体的过程如下：
18.1)采用ansys-fluent软件，在general下time一栏选择瞬态transient；
19.2)在模型models下多项流multiphase一栏选择mixture模型，在模型models下viscous一栏选择k-epsilon两方程模型下的realizable k-epsilon模型，在模型models下
多项流multiphase一栏选择phase interaction-slip，设置滑动速度代数方程；
20.3)在材料materials下选择流体fluid，添加空气相，设置viscosity粘度、density密度并且新建材料雪相，设置雪相viscosity粘度与空气相一致并依据步骤一设置density密度、粒径；
21.4)在体网格cell zone conditions下，把网格计算域内设置为所需要的流体fluid、固体solid，在boundary conditions设置边界条件：入口边界为速度入口、计算域两侧及顶部为对称边界、地面及模型表面为无滑移壁面、出口边界为自由出口；
22.5)在boundary conditions-inlet(入口)接入步骤三的空气相边界条件：设置风剖面、湍动能及耗散率，在boundary conditions-inlet(入口)-snow雪相-multiphase下volume fraction接入步骤四的雪相体积分数；fluent-mixture模型采用有限体积法计算，由于风雪流只有两相，且两相的积分分数和为1，因此只需给出雪相的体积分数，空气相由软件计算无需设置；
23.6)在求解器设置solution下的选择方法methods中，选用simple算法；
24.7)在求解器设置solution下的选择monitors-residual中，设置连续性方程、动量方程以及k-epsilon的收敛精度，其中k-epsilon的收敛精度设置为1e-6
精度；
25.8)在求解器设置solution下的初始化initialization中，选用标准初始化standard initialization，在入口初始化；
26.9)最后在run calculation设置计算的时间步长、计算步数，开始计算，得到了壁面剪切力、空气相及雪相的速度及体积分数的收敛文件。
27.更进一步地，步骤八中，侵蚀量的计算公式为：其中a
ero
为常系数，u
*
为入口摩擦速度，u
*t
为阈值速度。
28.更进一步地，步骤八中，沉积量的计算公式为：其中wf中为雪的沉降速度，φs为雪的质量浓度，φs＝ρsf，f为雪的体积分数，ρs为雪颗粒密度。
29.更进一步地，步骤八中，积雪改变厚度的计算公式为：其中qs为q
er
o或q
dep
，γ为积雪的最大体积分数，ρs为雪颗粒密度。
30.本发明的方法运用icem cfd软件进行模型的建立，将实测所测得的立方体及周边积雪分布运用于模型的建立中，考虑风致雪漂移的初始积雪分布，并采用ansys fluent进行数值模拟分析，本方法基于euler-euler模型的两相流理论，结合运用visual studio改进的湍流模型的数值计算开发分析程序，将实测所测得的风速、风向、雪颗粒的物理特性(形状、密度、尺寸)、风雪漂移的阈值速度、雪颗粒的降落速度、风雪漂移的运动模式、休止角等数据运用于模拟中，对建筑的积雪分布进行模拟，完成结合实测的风致雪漂移的cfd数值模拟。将模拟结果与风致雪漂移经典实验进行对照，并同未与实测相结合的方法进行对照，验证了该模拟方法的准确性。本发明考虑了风雪之间的相对运动并与实测进行的结合，提高了风致雪漂移数值模拟的准确性。
31.本发明的结合实测的风致雪漂移的cfd数值模拟方法的整体流程示意如图1所示，本发明的结合实测的风致雪漂移的数值模拟方法具有下列优点：
32.1、考虑了风相与雪相之间的相对滑移速度，通过求解两相之间的相对速度来考虑
两相之间的作用。
33.2、将实测获得的风吹雪的积雪分布形态运用在模型的建立中，考虑了雪相的沉积与侵蚀对计算域内流场分布的影响。
34.3、运用实验工具测得了各项基于当地条件的数据并将其运用在数值模拟中，使得模拟结果与实际更为贴合，更为准确。
35.本发明的方法有助于认识风雪运动的规律，有助于预防建筑倒塌、减少人员伤亡及经济损失。该方法可用于寒冷地区的建筑结构设计领域。
附图说明
36.图1是本发明的结合实测的风致雪漂移的cfd数值模拟方法的整体流程示意图；
37.图2是实施例1的网格划分以及计算域的图；
38.图3是实施例1的经步骤七得到的积雪分布形式图；
39.图4是实施例1得到的单位时间内的积雪通量、经典实验得到的单位时间内的积雪通量及实施例2即未与实测相结合的方法得到的单位时间内的积雪通量的对照图。
具体实施方式
40.用下面的实施例验证本发明的有益效果：
41.实施例1：本实施例为结合实测的风致雪漂移的数值模拟方法，按以下步骤进行：
42.步骤一：收集数值模拟需要的基本信息：收集具体研究的建筑物的几何特征：边长为1m的立方体模型，该建筑所处的地理位置地貌参数为0.16、参考高度10m处平均风速5m/s，风剖面采用公式z0与v0分别表示参考点的高度(1m)及对应的参考平均风速(5m/s)，z表示任一高度；α表示地貌参数，对于b类地貌，α取值为0.16；湍动能采用公式湍能耗散率采用公式其中，c
μ
为经验常数，取0.09；i和l分别为湍流强度和湍流长度尺度，其中：
43.湍动能强度采用公式积分湍流尺度采用公式其中，zg、zb、α'、均为地貌参数，在日本aij规范中，b类地貌对应ⅱ类平地，以上三个地貌参数取值分别为5m、350m和0.15。
44.通过采用分样筛、锥形漏斗、snow fork以及hobo小型气象站仪器测得当地的雪颗粒直径0.00015m、休止角45
°
、雪密度145kg/m3、计算阈值摩擦速度为0.2m/s、跃移雪粒速度0.2m/s、地面粗糙度系数为0.25、雪粒沉降速度0.2m/s。
45.步骤二：建立数值模拟需要的网格模型：在ansys-icem cfd中按照建筑物在步骤一所收集的建筑物几何信息和初始积雪分布形态进行模型的建立以及网格的划分，确定模型计算域，长、宽、高分别为20m、10m、8m，且计算域设置阻塞率为1.25％，接着在icem cfd中依据模型建立网格模型，依据模型大小网格数为360万。网格模型控制结构化网格质量在0.8，非结构化网格质量在0.8，最后以ansys-fluent格式输出求解的网格文件。网格划分以及计算域的图如图2所示；图2中1为对称边界，2为自由发展边界，3为速度入口，4为无滑移壁面。
46.步骤三：编写数值模拟用的空气项边界条件：利用步骤一所收集的建筑物所在地理位置信息、参考高度平均风速、风剖面、湍动能、湍动强度、积分湍流尺度、湍流耗散率规范信息编写空气相的边界条件，在visual studio中以c语言、fluent二次开发接口为基础，编写指数型风剖面、湍动能、湍流耗散率，保存c语言文件，并将其与步骤二网格模型文件放在同一文件夹中；
47.步骤四：编写数值模拟用的雪相边界条件：利用步骤一所收集的雪颗粒密度145kg/m3、当地重力加速度9.81m/s2、阈值磨擦速度0.2m/s，采用pomeroy&gray(1990)、pomeroy&male(1992)所提出的雪浓度经验公式，即公式其中ρs表示雪颗粒密度，g为重力加速度，取9.81m/s2，u
p
为跃移层雪粒平均速度，u
p
＝2.8u
*t
，u
*
为入口摩擦速度；根据beyers等实测拟合得到的摩擦速度与10m高度处风速u
10
的关系确定，雪相运动速度取5m/s、雪颗粒摩擦速度取0.15m/s、u
*t
为阈值速度；取0.15m/s；hs为跃移层和悬移层的临界高度，雪相运动速度5m/s、雪颗粒摩擦速度0.15m/s、在visual studio中以c语言、fluent二次开发接口为基础编写雪相体积分数公式作为雪项的入口边界条件，保存c语言文件，并将其与步骤二模型网格文件、步骤三的c语言文件放在同一文件夹中；
48.步骤五：数值模拟需要的程序接入，具体操作如下：打开ansys-fluent界面选择3d、双精度、并行。在general options-working directory下选择包含步骤二的网格文件、步骤三、四的自定义程序的文件夹目录，选择环境下set up compilation environment for udf接口打开软件。选择ansys-fluent主菜单下user-defined-functions-compiled下选择步骤三、四的c语言文件进行加载编译，得到空气相及雪相的入口速度及体积分数；
49.步骤六：采用ansys-fluent软件，先设置模型、材料、求解器参数，再导入步骤二所得的网格模型、步骤五所得的空气相及雪相的入口速度及体积分数，进行风雪两相流模拟，得到壁面剪切力、空气相及雪相的速度及体积分数的收敛文件，具体如下：
50.1)在general下time一栏选择瞬态transient；
51.2)在模型models下多项流multiphase一栏选择mixture模型、在模型models下viscous一栏选择k-epsilon两方程模型下的可实现realizable k-epsilon模型，在模型models下多项流multiphase一栏选择phase interaction-slip，设置相对滑动速度0.25m/s；
52.3)在材料materials下选择流体fluid，添加空气相，设置viscosity粘度、density密度，并且新建材料雪相，设置雪相viscosity粘度与空气一致，依据步骤一设置density密度145kg/m3(由snow fork雪特性分析仪测得)、粒径0.00015m(由分样筛测得)；
53.4)在体网格cell zone conditions下的网格计算域内设置需要的流体fluid、固体solid，在boundary conditions设置边界条件：入口边界为速度入口、计算域两侧及顶部为对称边界、地面及模型表面为无滑移壁面、出口边界为自由出口；
54.5)在boundary conditions-inlet(入口)接入步骤三的空气相边界条件：设置风剖面、湍动能及耗散率，在boundary conditions-inlet(入口)-snow雪相-multiphase下volume fraction接入步骤四的雪相体积分数；fluent-mixture模型采用有限体积法计算，由于风雪流只有两相，且两相的积分分数和为1，因此只需给出雪相的体积分数，空气相由软件计算无需设置。
55.6)在求解器设置solution下选择方法methods中选用simple算法。
56.7)在求解器设置solution下选择monitors-residual中设置连续性方程、动量方程以及k-epsilon的收敛精度，本实施例设置收敛精度为1e-6
精度。
57.8)在求解器设置solution下初始化initialization中选用标准初始化standard initialization在入口初始化。
58.9)最后计算在run calculation设置计算的步数为1000步，精度为10-6
，开始计算，得到壁面剪切力、空气相及雪相的速度及体积分数的收敛文件。
59.步骤七：利用步骤六计算出的收敛文件，进行壁面剪切力和速度的后处理：将ansys-fluent计算出的收敛文件以tecplot支持格式保存，用tecplot进行后处理计算，通过提取计算的计算域三维空间中x、y、z方向的速度、压力、壁面剪切力的分量，用tecplot自定义函数specify equations，做出计算域范围内的压力、速度、壁面剪切力的云图、等值线图；其中积雪分布形式图如图3所示。
60.步骤八：数值模拟雪沉积侵蚀的后处理计算：基于步骤七计算获得的壁面剪切力：基于步骤七计算获得的壁面剪切力，通过tecplot提取模型底面的壁面剪切力，通过自定义函数把壁面剪切应力转换为壁面的摩擦速度，将壁面摩擦速度与阈值磨擦速度进行比较，当壁面摩擦速度大于阈值磨擦速度时为侵蚀，采用公式计算侵蚀量，其中a
ero
为常系数，取7.0
×
10-4
；当壁面摩擦速度小于阈值磨擦速度时沉积，采用公式计算沉积量，其中wf为雪的沉降速度，φs为雪的质量浓度，φs＝ρsf，f为雪的体积分数，ρs为雪颗粒密度；再采用公式计算积雪改变厚度，其中qs为q
ero
或q
dep
，γ为积雪的最大体积分数，取0.62；最后将积雪改变厚度与步骤一的建筑物的初始积雪深度进行叠加，得到积雪的最终分布形态，完成结合实测的风致雪漂移的cfd数值模拟。
61.实施例2：本实施例在未结合实测、不考虑风相和雪相的相对滑移速度条件下对风致雪漂移的数值进行模拟，具体的方法按以下步骤进行：
62.步骤一：收集数值模拟需要的基本信息：具体研究的建筑物的几何特征：边长为1m
的立方体模型，该建筑所处的地理位置地貌参数为0.16、参考高度10m处平均风速5m/s，风剖面采用公式z与v0分别表示参考点的高度(1m)及对应的参考平均风速(5m/s)；α表示地貌参数，对于b类地貌，α取值为0.16；湍动能采用公式湍能耗散率采用公式其中，c
μ
为经验常数，取0.09；i和l分别为湍流强度和湍流长度尺度，湍动能强度采用公式积分湍流尺度采用公式其中，zg、zb、α'、均为地貌参数，在日本aij规范中，b类地貌对应ⅱ类平地，以上三个地貌参数取值分别为5m、350m和0.15。
63.通过经验采用相关数值：雪颗粒直径0.00015m、不考虑休止角、雪密度150kg/m3、计算阈值摩擦速度为0.2m/s、跃移雪粒速度0.2m/s、地面粗糙度系数为0.25、雪粒沉降速度0.2m/s。
64.步骤二：建立数值模拟需要的网格模型：在ansys-icem cfd中按照建筑物在步骤一所收集的建筑物几何信息和初始积雪分布形态进行模型的建立以及网格的划分，确定模型计算域，长、宽、高分别为20m、10m、8m，且计算域设置阻塞率为1.25％，接着在icem cfd中依据模型建立网格模型，依据模型大小网格数为360万。网格模型控制结构化网格质量在0.8，非结构化网格质量在0.8，最后以ansys-fluent格式输出求解的网格文件。
65.步骤三：编写数值模拟用的空气项边界条件：利用步骤一所收集的建筑物所在地理位置信息、参考高度平均风速、风剖面、湍动能、湍动强度、积分湍流尺度、湍流耗散率规范信息编写空气相的边界条件，以现有技术在visual studio中以c语言、fluent二次开发接口为基础，编写指数型风剖面、湍动能、湍流耗散率，保存c语言文件，并将其与步骤二网格模型文件放在同一文件夹中；
66.步骤四：编写数值模拟用的雪相边界条件：利用步骤一所收集的雪颗粒密度150kg/m3、当地重力加速度9.81m/s2、阈值磨擦速度0.2m/s，采用pomeroy&gray(1990)、pomeroy&male(1992)所提出的雪浓度经验公式，即公式其中ρs表示雪颗粒密度，g为重力加速度，取9.81m/s2，u
p
为跃移层雪粒平均速度，u
p
＝2.8u
*t
，u
*
为入口摩擦速度，根据beyers等实测拟合得到的摩擦速度与10m高度处风速u
10
的关系确定，雪相运动速度取
5m/s、雪颗粒摩擦速度取0.15m/s、u
*t
为阈值速度；取0.15m/s；hs为跃移层和悬移层的临界高度，雪相运动速度5m/s、雪颗粒摩擦速度0.15m/s、在visual studio中以c语言、fluent二次开发接口为基础编写雪相体积分数公式作为雪项的入口边界条件，保存c语言文件，并将其与步骤二模型网格文件、步骤三的c语言文件放在同一文件夹中；
67.步骤五：数值模拟需要的程序接入，具体操作如下：将步骤二的网格文件、步骤三、四的自定义程序存放相同文件夹，打开ansys-fluent界面选择3d、双精度、并行。在general options-working directory下选择上述文件夹目录，选择环境下set up compilation environment for udf接口打开软件。选择ansys-fluent主菜单下user-defined-functions
[0068]-compiled下选择步骤三、四的c语言文件进行加载编译，得到空气相及雪相的入口速度及体积分数；
[0069]
步骤六：将数值模拟软件的模型、材料、求解器参数进行具体设置，采用ansys-fluent导入网格模型、基于步骤五导入的二次开发编程进行风雪两相流模拟，具体如下：
[0070]
1)在general下time一栏选择瞬态transient；
[0071]
2)在模型models下多项流multiphase一栏选择mixture模型、在模型models下viscous一栏选择k-epsilon两方程模型下的可实现realizable k-epsilon模型，在模型models下多项流multiphase一栏选择phase interaction-slip，此时不考虑风相与雪相的相对滑移速度；
[0072]
3)在材料materials下选择流体fluid，添加空气相，设置viscosity粘度、density密度，并且新建材料雪相，设置雪相viscosity粘度与空气一致，依据步骤一设置density密度150kg/m3、粒径0.00015m；
[0073]
4)在体网格cell zone conditions下的网格计算域内设置需要的流体fluid、固体solid，在boundary conditions设置边界条件：入口边界为速度入口、计算域两侧及顶部为对称边界、地面及模型表面为无滑移壁面、出口边界为自由出口；
[0074]
5)在boundary conditions-inlet(入口)接入步骤三的空气相边界条件：设置风剖面、湍动能及耗散率，在boundary conditions-inlet(入口)-snow雪相-multiphase下volume fraction接入步骤四的雪相体积分数；fluent-mixture模型采用有限体积法计算，由于风雪流只有两相，且两相的积分分数和为1，因此只需给出雪相的体积分数，空气相由软件计算无需设置；
[0075]
6)在求解器设置solution下选择方法methods中选用simple算法。
[0076]
7)在求解器设置solution下选择monitors-residual中设置连续性方程、动量方程以及k-epsilon的收敛精度，本实施例设置收敛精度为1e-6
精度。
[0077]
8)在求解器设置solution下初始化initialization中选用标准初始化standard initialization在入口初始化。
[0078]
9)最后计算在run calculation设置计算的计算步数步数为1000步，精度为10-6
，开始计算，得到壁面剪切力、空气相及雪相的速度及体积分数的收敛文件。
[0079]
步骤七：利用步骤六计算出的收敛文件，进行壁面剪切力和速度的后处理：在ansys-fluent计算收敛后文件以tecplot支持格式保存。用tecplot后处理计算，通过提取
计算的x、y、z的速度、压力、壁面剪切力的分量，用tecplot自定义函数specify equations，做出计算域范围为内的压力、速度、壁面剪切力的云图、等值线图。
[0080]
步骤八：数值模拟雪沉积侵蚀的后处理计算：基于步骤七计算获得的壁面剪切力，通过tecplot提取模型底面上壁面剪切力，通过自定义函数把壁面剪切应力转换为壁面的摩擦速度，对壁面摩擦速度与阈值磨擦速度进行比较，当壁面摩擦速度大于阈值磨擦速度时侵蚀，采用公式其中a
ero
为常系数，取7.0
×
10-4
；当壁面摩擦速度小于阈值磨擦速度时沉积，采用公式其中wf中为雪的沉降速度，φs为雪的质量浓度(φs＝ρsf，f为雪的体积分数，ρs为雪颗粒密度)计算侵蚀量、沉积量；基于步骤一的初始积雪深度，采用公式计算积雪改变厚度、叠加积雪改变量(其中，统一用qs代替q
ero
和q
dep
，γ为积雪的最大体积分数，取0.62)，得到积雪的最终分布形态，完成风致雪漂移的数值模拟。
[0081]
将实施例1的模拟结果显示的单位时间内的积雪通量、实施例2在未与实测相结合、不考虑风相和雪相的相对滑移速度的条件下得到的单位时间内的积雪通量及经典实验的单位时间内的积雪通量曲线图进行比较，如图4所示，从图4可以看出，实施例2在未与实测相结合、不考虑风相和雪相的相对滑移速度的条件下得到的单位时间内的积雪通量与经典实验方法得到的曲线相差较大，而实施例1的模拟结果与经典试验的结果能很好地吻合，说明实施例1的模拟方法的准确性高。

技术特征：
1.一种结合实测的风致雪漂移的cfd数值模拟方法，其特征在于该方法按以下步骤进行：步骤一：收集数值模拟需要的基本信息：收集具体研究的建筑物的几何特征、该建筑所处的地理位置、地貌参数、参考高度h0处平均风速v0、风剖面规范、湍动能规范、湍动强度规范、积分湍流尺度规范、湍能耗散率规范；同时测得建筑物所在地积雪的雪颗粒密度、雪颗粒直径d
s
、休止角、风速、风向及建筑物的初始积雪深度，计算阈值摩擦速度、跃移雪粒速度、地面粗糙度、建筑物墙面粗糙度、雪面粗糙度和雪粒沉降速度w
f
；步骤二：建立数值模拟需要的网格模型：利用步骤一所收集的建筑物几何信息和初始积雪分布形态进行模型的建立以及网格的划分，得到网格文件；步骤三：编写数值模拟用的空气项边界条件：利用步骤一所收集的建筑物所在地的地理位置信息、参考高度平均风速、风剖面、湍动能、湍动强度、积分湍流尺度、湍流耗散率规范信息编写空气相的边界条件，再在visual studio中以c语言、fluent二次开发接口为基础，编写指数型风剖面、湍动能、湍流耗散率程序，保存c语言文件，并将其与步骤二网格模型文件放在同一文件夹中；步骤四：编写数值模拟用的雪相边界条件：利用步骤一所收集的雪颗粒密度、建筑物所在地的重力加速度、阈值摩擦速度，采用雪浓度经验公式及其公式中所采用的雪相运动速度、雪颗粒摩擦速度、跃移和悬移层临界高度信息，在visual studio中以c语言、fluent二次开发接口为基础编写雪相体积分数公式，作为雪项的入口边界条件，保存c语言文件，并将其与步骤二模型网格文件、步骤三的c语言文件放在同一文件夹中；步骤五：数值模拟需要的程序接入，具体操作如下：打开ansys-fluent界面选择2d/3d、双精度、并行，在general options-working directory下选择包含步骤二的网格文件、步骤三的c语言文件、步骤四的c语言文件的文件夹目录，选择环境下set up compilation environment for udf接口打开软件；选择ansys-fluent主菜单下的user defined
–
functions
–
compiled，选择步骤三、四的c语言文件进行加载编译，得到空气相及雪相的入口速度及体积分数；步骤六：采用ansys-fluent软件，先设置模型、材料、求解器参数，再导入步骤二所得的网格模型、步骤五所得的空气相及雪相的入口速度及体积分数，进行风雪两相流模拟，得到壁面剪切力、空气相及雪相的速度及体积分数的收敛文件；步骤七：利用步骤六计算出的收敛文件，进行壁面剪切力和速度的后处理：将ansys-fluent计算出的收敛文件以tecplot支持格式保存，用tecplot进行后处理计算，通过提取计算的计算域三维空间中x、y、z方向的速度、压力、壁面剪切力的分量，用tecplot自定义函数specify equations，做出计算域范围内的压力、速度、壁面剪切力的云图、等值线图；步骤八：数值模拟雪沉积侵蚀的后处理计算：基于步骤七计算获得的壁面剪切力，通过tecplot提取模型底面或底面横、纵方向中线上的壁面剪切力，通过公式把壁面剪切力转换为壁面的摩擦速度，其中τ为壁面剪切力，ρ为空气密度；对壁面摩擦速度与阈值磨擦速度进行比较，当壁面摩擦速度大于阈值磨擦速度时，积雪发生侵蚀；当壁面摩擦速度小于阈值磨擦速度时，积雪发生沉积；分别计算侵蚀量、沉积量，并且计算积雪改变厚度，最后将
积雪改变量与建筑物的初始积雪深度进行叠加，得到积雪的最终分布形态，完成结合实测的风致雪漂移的cfd数值模拟。2.根据权利要求1所述的一种结合实测的风致雪漂移的cfd数值模拟方法，其特征在于步骤二中所述的网格文件的求解步骤是：在ansys-icem中按照建筑物几何特征确定模型计算域且计算域设置满足阻塞率小于3％，接着在icem中依据模型建立网格模型，再依据模型大小设置合理的网格尺寸，调控网格偏斜率为三角形与四面体不大于0.95，调控网格纵横比，调控网格压扁程度在0-1之间，最后以ansys-fluent格式输出求解的网格文件。3.根据权利要求2所述的一种结合实测的风致雪漂移的cfd数值模拟方法，为保证计算顺利收敛及模拟的准确性，其特征在于所述的合理的网格尺寸是指：网格模型控制结构化网格质量在0.6～1，非结构化网格质量在0.3～1。4.根据权利要求2所述的一种结合实测的风致雪漂移的cfd数值模拟方法，其特征在于所述的网格纵横比为：流动核心区内网格纵横比小于5:1，边界层网格纵横比小于10:1。5.根据权利要求1、2、3或4所述的一种结合实测的风致雪漂移的cfd数值模拟方法，其特征在于步骤六中所述的风雪两相流模拟，具体的过程如下：1)采用ansys-fluent软件，在general下time一栏选择瞬态transient；2)在模型models下多项流multiphase一栏选择mixture模型，在模型models下viscous一栏选择k-epsilon两方程模型下的realizable k-epsilon模型，在模型models下多项流multiphase一栏选择phase interaction-slip，设置滑动速度代数方程；3)在材料materials下选择流体fluid，添加空气相，设置viscosity粘度、density密度并且新建材料雪相，设置雪相viscosity粘度与空气相一致并依据步骤一设置density密度、粒径；4)在体网格cell zone conditions下，把网格计算域内设置为所需要的流体fluid、固体solid，在boundary conditions设置边界条件：入口边界为速度入口、计算域两侧及顶部为对称边界、地面及模型表面为无滑移壁面、出口边界为自由出口；5)在boundary conditions-inlet接入步骤三的空气相边界条件：设置风剖面、湍动能及耗散率，在boundary conditions-inlet-snow雪相-multiphase下volume fraction接入步骤四的雪相体积分数；fluent-mixture模型采用有限体积法计算，由于风雪流只有两相，且两相的积分分数和为1，因此只需给出雪相的体积分数，空气相由软件计算无需设置；6)在求解器设置solution下的选择方法methods中，选用simple算法；7)在求解器设置solution下的选择monitors-residual中，设置连续性方程、动量方程以及k-epsilon的收敛精度，其中k-epsilon的收敛精度设置为1e-6
精度；8)在求解器设置solution下的初始化initialization中，选用标准初始化standard initialization，在入口初始化；9)最后在run calculation设置计算的时间步长、计算步数，开始计算，得到了壁面剪切力、空气相及雪相的速度及体积分数的收敛文件。6.根据权利要求1、2、3或4所述的一种结合实测的风致雪漂移的cfd数值模拟方法，其特征在于步骤八中，侵蚀量的计算公式为：其中a
ero
为常系数，u
*
为入口摩擦速度，u
*t
为阈值速度。7.根据权利要求1、2、3或4所述的一种结合实测的风致雪漂移的cfd数值模拟方法，其
特征在于步骤八中，沉积量的计算公式为：其中w
f
中为雪的沉降速度，φ
s
为雪的质量浓度，φ
s
＝ρ
s
f，f为雪的体积分数，ρ
s
为雪颗粒密度。8.根据权利要求1、2、3或4所述的一种结合实测的风致雪漂移的cfd数值模拟方法，其特征在于步骤八中，积雪改变厚度的计算公式为：其中q
s
为q
ero
或q
dep
，γ为积雪的最大体积分数，ρ
s
为雪颗粒密度。

技术总结
一种结合实测的风致雪漂移的CFD数值模拟方法，它涉及风致雪漂移的CFD数值模拟方法。它是要解决现有的数值模拟方法不考虑风雪之间的相对运动、准确性差的技术问题。本方法：基于现场实测测得的初始积雪分布形态建模及实测雪颗粒数据，考虑风相与雪相的相对滑移速度，采用编程软件Visual Studio编写UDF导入有限元数值模拟分析软件Ansys Fluent对结构的积雪不均匀分布进行模拟，运用软件Tecplot对模拟数据进行后处理，建立了一种结合实测的考虑风雪相对运动的风致雪漂移的CFD数值模拟方法。本发明提高了数值模拟方法的准确性，可用于寒冷地区的建筑结构设计领域。于寒冷地区的建筑结构设计领域。于寒冷地区的建筑结构设计领域。


技术研发人员：李方慧 魏治平 向本军 王程 苏政忠
受保护的技术使用者：黑龙江科技大学
技术研发日：2022.02.21
技术公布日：2022/5/25