CFD稳态数值模拟的建议

近年来，建筑计算风工程的研究和应用得到了很大的进步，但其数值计算的精度非常重要。数值模拟是一种近似解，误差的大小决定了求解的精度，误差主要产生于三个方面：模型误差、离散误差和迭代误差。下面主要根据数值模拟方面的实践和体会，并参考一些资料，从数值模拟计算域的尺寸、计算网格、对流项插值阶数、湍流模型、数值模拟结果的判断等几个方面提出一些建议。同时还选取一个立方体实测模型的结果进行比较验证。

(a) 计算域设置

在对建筑物表面风压进行数值模拟时，是将一个无限大的空间用一个有限的计算域来代替，即在距离建筑较远的地方人为设置几个避免，使求解于封闭，并保证这些壁面设置不会对建筑表面风压数值计算结果产生影响，即求解域的大小不宜太小，但也不宜太大以免增加计算量。从影响建筑物壁面风压考虑，对低矮建筑物（包括大跨度建筑），设h为建筑物的高度，建议入口距建筑物迎风面保证4h~5h的距离，建筑物侧面和顶面距各自流域边界的距离应大于4h。此时最大阻塞率小于3%。但是，高层建筑与低矮建筑物有所区别，因为低矮建筑以顶面绕流为主，而高层建筑则以侧面绕流为主。高层建筑计算域的高度H可小于3h，而计算域的水平宽度B应大于8倍建筑物宽度，此时阻塞率小于5%。背风壁面距出口的距离应使湍流充分发展，所以出口应距建筑物远一些，一般要求9h~10h。若距离太小，出口处有回流，则计算会出现发散。在大尺度建筑物平均风压模拟时，有时也可适当减少背风壁面距离，因为一般远场的网格较粗，湍流耗散较快，并且输运方程中都以对流项为主，较远下游的流动对上游影响较小，所以大多取7h~8h就可基本消除人为设置出口边界的影响。

(b) 计算域网格设置

进行CFD数值模拟计算时，首先要将计算区域离散化，即网格划分，数值计算是在离散网格点上满足流体动力学基本方程，因此网格划分将对数值模拟结果有直接影响。网格划分对计算精度的影响包括网格类型和网格尺度两方面。网格类型分为结构化网格和非结构化网格，一般建筑物体型较为复杂，生成结构化网格比较困难，在实际应用时往往采用非结构化网格，但对于简单形体，应该采用结构化网格，因为当流体流动和网格成一条直线时数值耗散最小。网格尺度的影响主要体现在近壁面网格的密度上，近壁面网格距壁面的距离可以用无量纲距离表征。

$y^+=\frac{\rho u_\tau y}{\mu}$
式中uz为摩擦速度，y为第一排网格节点距壁面的法向距离，ρ和μ分别为流体的密度和动力粘性系数。

由于非平衡壁面函数基于Launder和Spalding提出的平均速度对数率强烈的依赖于压力梯度理论，而壁面速度的对数率仅在y+>30~60范围的边界层中适用，因此，在运用壁面函数来模拟近壁面流动时，网格布置的关键是确定第一排网格的位置，一般原则是第一排网格应布置在y+=30左右而不宜太密。数值风洞的模拟结果很大程度上取决于计算域网格划分的优劣。因此网格应当精细到足够捕捉到剪切层、旋涡等物理现象的特征变化，网格的质量也必须足够的高。在流动变化梯度大的区域，网格的拉伸或压缩率应当小，使得截面误差小，两个相邻单元的扩展系数应当小于1.3。在大型建筑中，流动参数沿壁面法线方向变化剧烈，因此在壁面法向，网格需要加密，以精确求解壁面边界层内的压力。就计算网格形状而言，六面体优于四面体，前者可以引入较小的截断误差，并显出良好的收敛性。然而，在风工程中几何体总是非常复杂的，通常使用四面体网格。但在壁面处网格线必须与壁面正交，因此棱柱体网格和四面体网格应一起使用。如果想预测壁面的压力，至少应布置5层棱柱体网格在壁面处。

用有限体积发进行数值计算时，对流项的插值格式非常重要。一阶迎风格式通常包含较大的数值耗散。目前推荐采用二阶迎风插值格式，其数值耗散明显低于一阶迎风格式，具有较高的精度。对于强烈变化的流场，二阶迎风常常会产生数值振荡，此时可采用接近二阶的混合格式。

(c) 湍流物理模型的选取

在建筑计算风工程中，采用比较多的湍流物理模型有RNG κ-ε、SST κ-ω、SSG-RSM、BSL-RSM四种，它们分别是标准κ-ε、κ-ω模型、雷诺应力模型的改进模型，在分离流、旋涡等复杂流场模拟中有较强功能。通过对多个大跨度屋盖的数值模拟与风洞试验对比，建议在模拟建筑物表面平均风压时，一般采用SST κ-ω模型或BSL-RSM雷诺应力模型。当建筑物体型比较复杂或计算机配置较低时，优先采用SST κ-ω模型，当建筑物体型比较规则或计算机配置较高时，可采用BSL-RSM雷诺应力模型。

(d) 稳态数值模拟的边界条件

1）入口

在数值风洞入口处，通常指定等效的大气层边界。平均风速通过对数率或指数率的剖面得到，对来流风方向地面的粗糙长度z0，指数率的平均风剖面为

$\frac{\bar{u}(z)}{\bar{u_b}}={(\frac{z}{z_b})}^\alpha$
其中，ub是z=zb处的参考速度，α是对应于不同地面粗糙度类别的指数。

入口处湍流强度边界条件可以通过湍流黏性或湍动能与湍流耗散率来描述。对两方程湍流物理模型来说，湍动能和湍流耗散率可以通过等效边界层假设给出：

$k\left(z\right)=1.5[I(z)u(z)]^2$
$\epsilon\left(z\right)=\frac{C_\mu^{0.75}{k(z)}^{1.5}}{l}$
式中，Cu=0.09，I(z)是入口处的湍流强度剖面，可以参考日本或澳大利亚规范等国规范选取，uz是平均流速，l是入口处的湍流特征尺度，可以选择为建筑物的高度。

2）计算域顶壁和侧壁面的边界

这些边界位置设置为自由滑移壁面。

3）出口

由于出流接近完全发展，采用完全发展出流边界条件，即流场任意物理量沿出口法向的梯度为零。

4）建筑物壁面和地面

采用无滑移壁面。为保证壁面附近边界层的对数运算不出错，壁面附近的网格离开壁面最近的计算节点应当置于离壁面距离大于壁面粗糙高度。

(e) 收敛性及结果初步判定

迭代计算以相邻两步之间的残差为收敛标准。残差的量级一般在第一次迭代后即可得到，最终的残差应当趋于零。一般情况下，残差下降4个量级即可认为结果收敛。同时可以监测感兴趣的位置，记录当地的流动变量。如果这些变量趋于常数，这个区域的解可以认为已经收敛。在数值风洞稳态模拟中，计算结果一般要满足以下几个要求：

（1）来流远场风剖面的保持性：平均风速剖面需保持不变，湍流强度也达到近似保持的要求。

（2）最后结果必须采用高精度的离散格式计算。

（3）无滑移壁面上得到的风速必须为零或接近零，使物理边界条件应用准确。

（4）最终解是网格无关解，可通过网格加密或网格优化得到网格无关解。

(f) CFD模拟的不确定性因素

数值模拟的结果不仅与湍流模型紧密相关，还很大程度上依赖于网格设计和数值格式。而后者在使用相同的程序的情况下，不可避免掺杂了人为因素。Crown等. (1997)总结了由欧盟多家研究机构参与的有关建筑绕流和扩散合作项目的部分结果，并通过“EUP项目”，分析比较了不同机构的网格划分、差分格式及湍流模型引起的差别，得出结论认为，计算结果很显著的依赖于网格设计参数、空间离散格式以及湍流模型。也就说明了，在实践中，数值模拟所面临的不确定因素，在某些情况下，都不会比由于“用户因素”所造成的差异更大。计算风工程中什么样的误差是可以接受的？这个问题没有一个明确的标准。理想状况是获得和场地实测数据完全一致的结果，但目前还难做到。数值模拟希望尽可能减小误差，使得CFD计算结果与场地实测数据尽量接近。但不同的研究者给出的可接受的误差范围也不同，并且可接受的误差范围还和应用的对象有关。目前，在计算流体动力学所依赖的理论本身和数值计算技术没有重大突破之前，期望CFD的结果与实测结果吻合非常好，既不现实也容易引起误导。处于大气边界层的钝体建筑结构，周围的流场高度湍流化，在时间和空间上呈现随机脉动的特征，因此，一般使用统计意义上的量来描述流动特征。对于这样复杂的现象进行数值模拟，不可避免的存在误差。