1 引言

尽管目前CFD 技术已经越来越多地应用于指导空调通风房间的气流组织设计及其评价分析,但对于实际工程而言,存在风口入流边界条件描述复杂、难以选择合适的湍流模型、迭代计算耗时等问题,影响了CFD 在实际工程中的应用。对于工程设计而言,人们总希望能在规划设计阶段用尽可能短的时间完成对室内空气流动的数值模拟,以及时指导和优化设计;而且,绝大多数暖通空调工程师并不具备速度很快的超级或大型计算机,他们希望能在普通的微机上较快地完成模拟。为此,本文提出一个简化的风口入流边界条件描述方法———N 点风口动量模型以及加速迭代计算收敛的新数值算法———误差预处理法,并结合一个新的MIT 零方程湍流模型,在保证工程精度要求的前提下,快捷地数值模拟空调通风室内空气流动, 以期推进. CFD在实际工程中的应用^[1]。

2 N点风口动量模型

N 点风口动量模型的基本思想是利用N 个简单开口替代形状复杂的、不同出流方向的送风口,每个开口采用动量模型的原理,既保持所描述风口的外形尺寸和位置不变,同时又保证正确模拟入流质量流量、动量流量以及浮力通量。空调风口射流的浮力通量为①，式中为重力加速度（m/s²）；为入流风量（m³/s），为送风温度T_in和室内平均温度之差（K），为室内平均温度（K）。实际空调通风房间的送风温度基本均匀，视为结合Boussinesq关于密度（温度）变化不大时其变化只作用于浮升力项的假设，可以不变，因此由式①知，只需保证入流质量（体积）流量和实际情况相等即可准确模拟入流浮力通量。难点在于如何准确模拟动量流量。多数实际的风口外形面积和有效通过面积并不一致，仅用简单开口代替实际风口，难以同时准确模拟质量流量和动量流量。N点风口模型采用如下描述方法：不同的出流方向代表不同的入流动量（动量为矢量），送风口出流不同方向中的每一个方向可称为一“束”，每一束为由送风口形成的一股或数股小射流的合成，出流不同方向束数取为N。即N=送风口出流方向束数②；对于一简单开口按如下方式描述入流动量③

式中为入流动量流量（kgm/m²）；为入流质量流量（kg/s）；为风口有效面积（m²）；为风口外型总面积（m²）；为有效面积系数，风口有效面积和外形总面积之比，1对于按式②划分的N个简单开口，只需对每个开口定义入流速度保证质量流量和实际情况相等，然后在计算动量源项时按式③附加送风口的有效面积系数，即实现了对入流动量流量的准确模拟^[4][2]。

3 MIT 零方程湍流模型

这是在室内空气自然对流和混合对流的直接数值模拟DNS结果的基础上提出的湍流模型, 该模型针对房间内非等温流动的Rayleigh 数范围(2.6～3.0×10¹⁰)为涡粘系数正比于流体密度、当地速度和距壁面最近之距离,比例系数由直接数值模拟的结果拟合而得=0.03874Qvρl, (2)其中: v 为当地时均速度, l 为当地距壁面最近的距离。该模型少求解2个微分方程,而仅求解关于质量、动量和能量守恒的5个微分方程, 故计算最省时间。其方程如下^[5][3]：

连续方程

动量方程

能量方程

上述各式中, 为 方向的时均速度; 为对于i= 1, 2, 3, x_i 代表三个垂直坐标轴坐标; 为空气密度; 为方向的时均速度;为空气压力; 为等效粘度; 为空气热膨胀系数; 为参考温度; 为空气温度; 为i 方向之重力加速度; 为空气焓值,是空气定压比热 与温度之乘积; 为等效热扩散系数; 为热源;为空气层流动力粘度, 是空气物性参数之一;为空气湍流动力粘度, 与湍流情况有关;为空气导热系数; 为空气定压比热; 为湍流普朗特数.涡粘系数模型的核心就是求解上述的湍流动力粘度(或称湍流粘性系数) , 根据附加求解的微分方程个数可将涡粘系数模型分为零方程模型、一方程模型、二方程模型等.本文选择MIT建筑技术系开发的新零方程模型(下简称M IT 零方程模型). 由于零方程模型耗时少、计算和收敛速度快, 因此很适合一般工程师和设计人员所具有的微型机使用.为了编制通用的计算程序,可以将室内空气流动用通用微分方程表示如下；其中代表通用变量、、、等，、、、分别表示密度、速度矢量、扩散通量和源项。扩散通量由下式确定：；其中表示通用变量的有效交换系数。在直角坐标系下，、、的值如下表所示。其中、、分别表示、、三个方向的速度，和分别表示空气的焓和压力，、、分别表示、、三个方向的重力加速度，为空气的参考密度，、、分别表示层流粘性系数，湍流粘性系数和有效粘性系数，为的当量普朗特数，表示单位体积的发热量。

4 误差预处理法

传统迭代法是在整个计算区域的所有网格上同时进行迭代计算至收敛。但如果能先求出一些点的值, 然后再根据它们来预测所有点的值, 只要预测方法得当, 从误差削减的角度而言, 就可以获得比迭代法高得多的效率。可先在粗网格上迭代计算至问题收敛, 然后借助该预测方法把粗网格上的值映射到细网格上作为迭代初值,再进行迭代计算直至收敛。该方法的基本思想在于通过在粗网格上用比较少的时间(相对于在细网格上迭代) 来获得问题的大致描述,以便通过线性插值大量削减误差,称其为“误差”。假设已知如图1计算域中点A ,B ,C, D,E,F,G,H 的真实值。各点坐标分别为: A (x - da, y + dy , z + dz ) , B (x + dx , y + dy , z + dz ) , C (x - da, y + dy , z - dc) , D (x + dx , y + dy , z - dc) , E (x - da, y - db, z + dz ) , F (x + dx , y - db, z + dz ) , G (x - da, y - db, z - dc) , H (x + dx , y - db, z - dc).

图1　计算域示意图

由流体流动的输运控制方程可知, 点I 上的值只受其相邻点上的值的影响, 即只受A , B , C, D , E , F , G, H 点值的影响。因此,可以用这些点上的值来预测一个误差比较小的I 的估计值。很自然用这8 个点的加权平均值作为I 的值: P(I) =M_Af(A )+M_Bf(B )+M_Cf(C)+M_Df(D )+M_Ef(E )+M_Ff(F)+M_Gf(G)+M_Hf(H ) ④ ；其中: P (I) 表示预测值, f(*)表示相应点的真实值, M_X表示加权系数。由守恒律有M_A+M_B+M_C+M_D+M_E+M_F+M_G +MH = 1⑤ 。离I 比较近的点对I 的影响应该比较大, 即较近的点加权系数应较大。结合式⑤可以根据面积律导出加权系数:

M_A = (dx×db×dc)/ [dx +da] (dy+db) (dz+dc)],

M_B = (da×db×dc)/ [（dx+da）(dy+db) (dz+dc)],

M_C = (dx×db×dz )/ [(dx+da)(dy+db) (dz+dc)],

M_D = (da×db×dz )/ [(dx+da)(dy+db) (dz+dc)],

M_E = (dx×dy×dc)/ [(dx+da)(dy+db) (dz+dc)],

M_F = (da×dy×dc)/ [(dx+da)(dy+db) (dz+dc)],

M_G = (dx×dy×dz )/ [(dx+da)(dy+db) (dz+dc)],

M_H = (da×dy×dz )/ [(dx+da)(dy+db) (dz+dc)]. ⑥

将式⑥代入式④, 得

P (x ,y ,z ) = [(da×db×dc)f(x+dx ,y+dy ,z+dz)+(da×db×dz )f(x+dx ,y+dy , z - dc)+(da×dy×dc)f(x+dx,y-db, z+dz)+(da×dy×dz)f(x+dx,y -db,z-dc)+(dx×db×dc)f(x-da,y+dy,z+dz)+(dx×db×dz)f(x-da,y+dy, z-dc)+(dx×dy×dc)f(x -da,y-db, z + dz )+(dx×dy×dz )f(x-da,y-db,z-dc)]/[(dx+da)(dy+db)(dz+dc)]. ⑦

显然, 实际物理问题的解都是简单函数。故由Taylor 公式, 有f (x +dx ,y +dy ,z + dz )=f (x ,y ,z ) + f_x(x ,y ,z )dx + f_y(x ,y ,z )dy +f_z(x ,y ,z )dz + O(dx²+dy ²+dz ²). ⑧

类似写出其它项, 把它们代入式⑦,并考虑到O (dx ²)=O(da²) ,O (dy ²)=O(db²) ,

O(dz²)=O(dc²),可以得出P (x,y,z )-f (x,y,z )=O (dx ²+dy ² +dz²).此式表明上述插值方法具有二阶精度, 而且它是显式方法, 效率远远高于迭代方法。综上所述,误差预处理法的算法如图2所示。

图2 中当运算进行到限制和插值时即采用上述的预测方法。由此可见, 误差预处理法本身不仅是数学方法, 而且也蕴涵了一定的物理意义^[7]。

5 数值计算方法

上面各微分方程相互耦合,具有很强的非线性特征,一般只能采用数值方法求解.在求解区域生成网格后,利用有限容积法可将个微分方程离散为以下形式:

⑨

式中a为离散方程的系数,为个网格节点的变化值,为离散方程的源项.下标、、、、、和分别表示本网格、东边网格、西边网格、北边网格、南边网格、上面网格和下面网格处的值.所有的动量方程、能量方程均可离散成上述形式。根据质量守恒方程，离散后得到下列形式的压力修正方程：

⑩

式中为压力修正方程的系数，为修正压力。于是用CFD方法模拟通风空调内温度分布的步骤如下：

Ⅰ 给所有变量附初值；

Ⅱ 求解离散的动量方程⑨，得到各坐标轴方向上的速度、和；

Ⅲ 求解压力修正方程⑩，得到修正的压力，用修正压力修正估计压力；

Ⅳ 用修正压力对步骤Ⅱ中得到的速度进行修正；

Ⅴ 求解离散的能量方程⑨得到空气的焓或温度；

Ⅵ 判敛，如收敛则结束，否则转步骤Ⅱ^{[6] [8]}。

6 结束语

根据该算法，清华大学建筑技术科学系自行开发了STACH-3可以很好的模拟出空调房间的速度分布和温度分布。以上模型相对其他模型具有简便易形的特点且可以在较短的时间内模拟室内空气的速度场和温度场，很适合工程师在实践中运用其解决实际问题。

参考文献

1 赵彬,李先庭,彦启森.室内空气流动的简捷数值模拟方法, 技术交流园地,2003,33(3);102-104

2 赵彬,李先庭,彦启森.室内空气流动数值模拟的风口模型综述, 国外技术介绍, 2000,30(5);33-37

3 赵彬,李先庭,彦启森.用零方程湍流模型模拟通风空调室内的空气流动, 清华大学学报(自然科学版),2001,41(10);109-113

4 赵彬,李先庭,彦启森.孔板空调风口送风射流的数值模拟,力学与实践,2002,24(1);18-21

5 赵彬,李先庭,彦启森.用CFD方法改进室内非等温送风气流组织设计, 应用基础与工程科学学报, 2000,8(4);376-386.

6 赵彬,李先庭,彦启森.用CFD方法指导通风空调设计,制冷与空调,2001,1(5);11-15

7 赵彬,李先庭,彦启森践,室内空气流动数值模拟的误差预处理, 清华大学学报2001,41(10);114-117.

8 赵彬,李先庭,胡斌,彦启森.计算流体力学在暖通空调工程中的应用,计算机应用与信息技术.

余亮,男,1982.1出生,在读硕士