基于N-S方程的计算流体力学算法分析

2022-04-29张荣芳

吕梁学院学报 2022年2期

张荣芳

(山西职业技术学院基础教学部，山西太原 030006)

1 计算流体力学数值模拟软件调研

Fluidyn-PANACHE主要用于对模拟空气运动的Navier-Stokes 方程的解进行求取，其求解的过程主要应用 3D 有限体积法(FVM)[1].Fluidyn-PANACHE还属于计算机模拟模型的一种，主要应用于对和大气过程存在关联的危害以及污染进行诊断[2，3].它使用 3D 计算流体动力学(CFD) 方法来计算复杂的风场、污染物扩散以及火灾和爆炸[4].当与 MM5 等其他诊断风模型结合使用时，可用于流行率和风险的诊断分析.

与高斯模型相比，PANACHE 主要应用确定性的解决方案技术，同时结合物理模型的应用，其中包括较少的经验和对场地、释放场景和污染特征的敏感性较低.PANACHE 内部存在一个3D 自动晶格生成器，它可以创建一个算术大小的有限尺寸网络，适用于障碍物和起伏地形周围的物体.网络不但可以是非结构化的三角形，同时也可以为结构化的矩形.为了达到点聚合的需要一部分区域能够应用嵌套场网络完成.软件中加入了特殊大气边界层(PBL)模型，考虑了地形、植被冠层、城市冠层等所有地形的影响.

IconCFD是基于OpenFOAM开发的CFD程序，可以自动生成收敛效果好的计算模型，修复OpenFOAM的一些已知bug，改进了新的概念模型.另外，在OpemFOAM的基础上，大大提升了页面的友好度和易操作性[5-6].

OpenFOAM是一个完全用C语言编写的面向对象编程的CFD库.它的求解器使用相对有限的体积法，适用于多面体网格，因此可以求解复杂的几何设计外观.得到的问题范围很广，不仅可以得到化学变化、渗流、热对流等复杂的流动性，还可以得到固态动力学模型和电磁场理论等，适用于大规模并行处理变换成网格图和计算.

2 计算流体力学算法分析

本文以高大工程建筑附近的复杂地貌和住宅小区为研究对象.采用非结构化网格图和有限体积法求解N-S方程，模拟势流.模拟湍流流动使用k-ε模型.

2.1 控制方程

CFD以Navier-Stokes方程为基本方程，而且还能够获得外来物种的能量守恒定律方程、浓度值以及质量守恒定律方程[7].鉴于颗粒物和大气气溶胶的传播，还需要用到拉格朗日描述的粒子适应度轨迹方程.对于渗流法，求上式的Leylo均值法.瑞利内应力实体模型主要应用线性涡粘性实体模型方程[1]：

(1)物质守恒浓度

(2)连续性方程

上述公式对全部物质的和进行求解，能够获得连续性方程(质量守恒)

(3)Navier-Stokes方程(动量守恒)

(4)能量守恒

2.2 大气边界层(ABL)方程

空气边界层(ABL)是靠近地面的空气区域.路面基于渗流交换弹丸运动、发热和水分.描述空气最大范围的方程没有考虑空气和地球表面之间的相互作用.这种相互影响体现在渗流健身运动上，其优点是极限小、分格图少，因此必须进行参数化设计(模型).考虑所有地质结构性危害，例如地貌波动、植物群落的冠层和大城市的冠层.

ABL实体模型用于在气象监测和 CFD 推导器所需的初始条件之间创建的接口.其构成主要包括以下两部分[8]：

(1)小气候实体模型，以观测数据信息作为参考，能够完成PBL的主要物理参数的计算；

(2)附着面层的实体模型需要风、温度和渗流的垂直截面.

2.3 湍流模型

本课题选择应用的模型为两方程线性涡黏性k-ε模型.在CFD模拟中主要应用可压缩性的标准高雷诺数格式以及修正浮力[2-3]，如下列方程所示，为通过该模型对湍流动能k的输运方程及其耗散率ε方程进行求解：

k-ε实体模型通过局部渗流特性测量涡流渗流脉冲的长度尺度，因此适用于机械设备剪切(障碍物、地形起伏、树干)和水浮力(可靠性和水浮力/重气).k-ε实体模型是各向异性渗流实体模型.因此，通过该模型测量某一位置的结果，不管是垂直还是水平角度的渗流扩展都是一致的[6].

图1 结构化和非结构化网格

2.4 边界条件

对方程式的解进行求取时，对于主要范围内的地面、周围和障碍物都要求进行边界条件的给定.其中，边界条件的类别，包括出口、入口以及壁面三种，分别为入口、出口和壁面，如表1所示.

表1 边界条件

一般情况下都将顶部边界作为出流边界，而以风向为参考都会将区域的侧边界作为出口边界和入口边界，如图2所示.

图2 流动边界条件

(1)入口.对于入口边界位置，要求给定温度变量，湍流变量，速度变量以及物种浓度变量.而根据局部范围内的内部通过计算就能够获取压力.以设定的物种背景浓度作为参考，能够完成物种浓度的设置.如下所示，为设置其他变量的方式.

设置入口边界位置的速度变量方式如下所示：

u=Ucosαv=Usinαw=0

式中，U代表的是高度为z位置的风速；α=-φ-90，该角度为屏幕x轴和风向两者之间的夹角；φ代表的是相对正北方的风速方向.以选择的垂直速度剖面作为参考，能够设置高度z位置的风速.风速设置完毕之后，又可以以垂直速度剖面作为参考，完成入口边界温度的设置.

(2)出口.虽然已经将某个边界位置设定为出口边界，然而因为冠层或者障碍物产生干扰，通过该位置风仍然能够进入该范围之内.以此为基础PANACHE能够对一部分变量进行指定，也能够通过内部推导的方式完成一部分变量的推导.

(3)壁面.PANACHE能够对湍流边界层内部固体壁面的拖曳力进行计算，计算过程中主要借助于壁面函数实现，而该函数的获得主要以平衡条件作为前提条件，通过湍流边界层的Navier-Stokes方程进行获取.

2.5 CFD求解器

获得非线性偏微分，例如质量守恒定律、动能和各种浓度值.在三维时间和空间范围内能够通过 CFD 获取器完成操作方程的获取.而且能够通过热对流扩散方程表示此操作方程，使用以下一般方法：

(I) (II) (III) (IV)

式中，ГØ为Ø的交换系数(黏度，热扩散系数等)；Ø为待求变量(动量、温度、浓度等)；SØ为准的源项(污染物排放等).(I)为时间微分项，(II)为对流项，(III)为耗散项，(IV)为源项.

要求立即获得上述类型的操作方程的结果在于为前三项选择的离散变量文件格式.随着离散变量元素逐渐变少(元素总量逐渐增多)，计算误差会增加.然而随着元素总量的提高，需要的会计成本也会有所提高，其中包括耗费的时间、运行内存以及CPU.所以，就需要在计算效率以及精度两者之间进行权衡.在上面的运算方程中，主要以泰勒级数作为参考，完成微分算子离散变量的执行.

(1)时间差分.通常选择应用近似离散的时间序列进行时间差分，如tn(n=0，1，2，…)，时间间隔Δtn=tn+1-tn为时间步长，在tn时刻Øn代表的是对变量的差分近似，微分项∂Ø/∂t近似用一阶差分Øn+1-Øn/tn表达.

(2)空间差分.通过借助于控制体形成的三维网格从而实现空间离散的.为了保障局部范围内查分变量的守恒大致上不同单元的有限差分都选择应用整体平衡方式.PANACHE求解速度向量笛卡尔坐标系下分量的Navier-Stokes方程.在相同的控制量之内能够求解获取包括温度、压力、湍流等在内的全部变量.

选择应用二阶精度方式实现耗散项的离散；一般情况下，求解的精度取决于对流项的离散，如果选择应用的是有限体积方式离散的情况下，可以通过如下所示的公式表示对流项：

3 案例分析

图1和图2分别给出利用CFD进行建模求解[4-5]的例子，对CFD模拟功能设置边界条件和求解类型进行分析.

求解类型：

(1)稳态模型求解；

(2)可进行等温或非等温的模型求解；

(3)可进行多组分传递模型的求解；

(4)可进行拉格朗日粒子(一种颗粒污染物，无相变)模型的求解；

(5)可考虑浮升力的影响；

(6)可进行多种湍流模型的求解，计算使用k-OmegaSST湍流模型，湍流模型参数采用默认值.

模型可实现的边界设置类型：

1)定值边界类型(如定速度、温度等)；

2)进口边界的廓线分布类型(如进口速度、温度以及湍流强度等)；

3)温度边界分区域设置(如地面按照坐标范围设置不同的温度)；

4)可考虑粗糙度的影响(如地面粗糙度)；

5)可设置粒子进口边界(如质量流量、颗粒尺寸分布等).

图3 计算域内部几何模型

图4 局部地面表面网格

为了保证计算结果的可靠性，需要进行计算收敛性的判断.收敛性的判断可采用以下方式进行.

(1)参考计算结果的残差残差的大小并不能决定计算结果的好坏，它只是迭代计算过程的相对误差，每种软件都有不同的定义方式.因此残差仅仅作为参考；

(2)检查监控点物理量的变化是否达到稳定由于实际的物理过程是非稳定的，在达到收敛后，直接监控得到的物理量也会呈周期性波动.因此，我们将一定迭代步数内的物理量进行平均，对平均后的物理量进行监控.例如，图5给出了某监测点CO的质量分数平均值随迭代步数的变化.考虑到前期的计算并不收敛，我们设置从第800迭代步开始进行平均.可以看出，在大约1500迭代步后CO的质量分数的平均值基本稳定.这可以作为判断收敛的标准之一；