暖通空调气流组织数值模拟的特殊性
清华大学 赵 彬☆ 李先庭 彦启森
摘要 从计算流体力学CFD 应用于暖通空调工程的角度, 阐述了通风空调房间内气流组织数值计算的特殊性和难点, 对风口模型、湍流模型、辐射模型、热源模型、计算收敛速度以及复杂物理条件的描述等作了分析和介绍, 并提出了相应的对策或解决思路。
关键词 计算流体力学 通风空调房间 气流组织 数值模拟
Pa rti c ul a rit y of n u m e ri c a l si m ul a ti o n f or a ir distri b uti o n
i n v e ntil a t e d a n d a ir c o n diti o n e d r o o m s
By Zhao Bin ★, Li X ianting and Y an Qisen
Abs t r a ct Elaborates t he p articularit y and difficulties of t he numerical simulation f or air dist ribution by means of comp utational fluid dynamics (C FD ) applied in HVA C engineering. A nalyses air op ening model , turbulence model , radiation model and heat source model , and p resents t he calculation of convergence sp eed and t he descrip tion of complicated p hysical conditions. Puts f orward corresp onding counter measures and solutions f or t he difficulties.
Keywor ds comp utational fluid dynamics , ventilated and air conditioned room , air dist ribution , numerical simulation
★Tsinghua University , Beijing , China
0 引言
调领域数值计算的特殊性。为将CFD 技术这一强有力的工具应用于暖通空调气流组织设计, 还需要对这些问题进行研究。下面将结合笔者和其他研究者的研究成果和经验, 逐一介绍之, 以期得到同行的指正。1 风口模型
1. 1 风口模型的重要性和不足
自1974年丹麦学者P. V. Nilsen 首次将计算流体力学CFD 技术应用于暖通空调室内空气流动的数值模拟以来, 众多学者在该领域付出了辛勤的劳动, 而大家只有一个奋斗目标:利用数值模拟技术快速、准确地解决暖通空调实际工程中的流动和传热问题。随着计算机技术、数值计算技术以及湍流理论的发展, 人们在暖通空调工程的CFD 应用研究领域取得了长足的进步。在短短的20多年内, 涌现出了各式各样关于空调通风房间室内空气流动数值模拟的算法、湍流模型以及对一些实际应用中的特殊问题的解决办法, 暖通空调的数值计算正在蓬勃发展当中。但是,CFD 技术本身还不完善, 比如人们对湍流现象的认识尚未完全清楚; 而且, 暖通空调实际工程问题比较复杂, 如送风口的入流边界条件描述, 辐射换热和热源的考虑, 实际人员、设备等边界条件的描述等; 另外, 在实际应用中还存在要求计算速度快、能在较短时间内模拟大量工况等问题。以上各种问题共同构成了暖通空
风口入流边界条件的描述是当前通风空调房间气流组织数值模拟的难题之一。从数学角度而言, 对于CFD 所求解的描述室内空气湍流输运的非线性偏微分方程组, 边界条件对计算结果的准确性有着极大的影响; 从物理角度而言, 空调通风房间的空气流动就是由风口射流引起的, 射流动量流量和质量流量对室内空气分布情况起关键作用。
☆赵彬, 男,1974年12月生, 博士研究生, 工学博士, 讲师
100084北京清华大学建筑学院建筑技术科学系
(010) 62779995
E 2mail :[email protected]. cn 收稿日期:20010103修回日期:2004
09
因此建立风口模型来合理描述房间入流边界条件
是准确、快速地模拟通风空调房间气流组织的关键。不同于CFD 技术在其他领域的应用, 暖通空调工程领域中室内空气流动数值模拟的风口入流边界条件是很复杂的。为了正确预测室内空气分布情况, 我们需知道送风参数的详细情况以正确描述数值模拟时的入流边界条件。而实际的送风口几何形状却很复杂, 如散流器、百叶风口和孔板风口等, 且送风射流雷诺数通常大于30, 为湍动射流[1]。这样的入流条件极其复杂, 难以简单描述。为了了解这些风口的送风参数, 在用CFD 方法进行数值模拟时计算区域内的网格必须划分至mm 甚至更小的量级; 而为了模拟室内空气流动情况, 反映室内湍流大涡旋对流动的影响, 计算区域必须达到m 的量级。这样, 对于实际应用中的三维问题而言, 必将导致计算区域内的网格节点数目巨大, 超过目前一般计算机的容量, 从而使得计算量无法为一般的工程应用所接受。因此, 在对通风空调房间进行室内空气流动的数值模拟时, 必须采用较少的计算网格来描述复杂的风口入流边界条件以适应工程应用中快速计算的要求, 同时又不失一定的准确度, 于是, 对风口模型的研究就成为CFD 技术应用于暖通空调领域内的独特之处, 同时也是暖通空调气流组织数值模拟需要研究的一个非常重要的问题。
自20世纪70年代起, 就有学者提出用风口模型来解决室内空气流动数值模拟的送风口边界条件描述问题[2]。经过20多年的发展, 逐渐形成了两类风口模型:间接描述类风口模型和直接描述类风口模型。间接描述类风口模型包括盒子模型[3]、指定速度模型[3～4]和主流区模型[5]等, 这类模型把送风口入流边界条件描述转化到包围送风口的一个盒子的边界上, 借助实验值或者射流公式描述边界条件, 避免了直接描述送风口出口处复杂的流动情况, 但这种方法需要实验测量, 非常麻烦, 且无法用于设计阶段的预测; 另外, 还受射流公式不可靠、受射流是否受限或等温的影响。直接描述类风口模型包括基本模型[6]、动量模型[7], 这类模型直接在送风口处描述入流边界条件, 具有不受送风射流等温与否、受限与否的限制等优点, 但是基本模型对于开口有效面积系数(开口有效面积与总面积之比) 较小的送风口误差较大; 动量模型使用方便, 但对于如何简化合
并不同出流方向的射流并不明确。1. 2 对策或解决思路
针对现有风口模型的不足, 笔者提出了N 点风口模型以适应工程中快速简捷模拟风口入流边界条件的需要。N 点风口模型在保证能准确模拟入流质量流量、动量流量以及浮力通量的前提下, 用N 个简单开口替代不同出流方向的送风口以实现其入流边界条件的简化描述, 属于直接描述的风口模型, N 的取值根据风口射流动量数进行简化合并, 有关详细情况可参见文献[8]。该风口模型对百叶风口常见不同出流条件、孔板这类有效面积小的风口都能获得满意的模拟结果, 下面举例说明之。
图1为某百叶送风口出流方向张角为60°
的
图1 百叶送风口60°张角出流的N 点模型(N =3)
N 点风口模型, 这是一种百叶风口比较常见的出
流条件。图2a
为此时该百叶风口在出流流速为
图2 百叶送风口60°张角出流下射流主体段
特性模拟值和实验值对比
2, 5, 8, 10m/s 时的主体段轴心速度衰减和实验
值的对比, 可见N 点风口模型模拟结果和实验值
吻合得很好; 图2b 为四种流速下的N 点风口模型模拟所得主体段断面量纲一的速度分布与实验结果的对比, 可见模拟所得结果符合送风口湍流射流主体段断面流速自相似分布的特点。图中v m 为射流轴心速度, v o 为射流出口速度, r 为断面上某考察点距轴心的径向距离, r 0. 5为断面上流速为轴心速度0. 5倍之点距轴心的径向距离。尽管N 点风口模型能解决常见的百叶送风口和多孔板送风口的出流边界条件描述问题, 比较适用于工程应用, 但是对于更为复杂的送风口边界条件尚存在较大误差, 需要寻找更通用的风口模型, 或者新的室内空气流动的模拟思路, 如分段模拟室内空气流动的思路[9]。2 湍流模型
2. 1 气流组织数值计算湍流模型的特点及不足
由于通风空调房间内空气流动多为湍流流动, 而人们对湍流机理的认识至今还不完全, 对于工程应用, 目前多为半经验的、唯象的模拟, 故此对湍流的模拟始终成为室内空气流动数值模拟的一个主要研究方向。限于目前的计算机能力, 工程中最常采用的是涡黏系数模型EVM (eddy viscosity
model ) 中的K ε2
两方程模型或其变形。而K 2ε模型对于等温流动情形能模拟得很好, 但是对于空调通风房间内的非等温、混合对流的情形却有较大误差[10～11]。Nielsen 于1998年指出, 对于热羽流、贴壁射流、温度分层流动等需要不同的湍流模型进行模拟方能取得满意的精度[12], 而这些流动形式往
往并存于实际的空调房间, 采用K ε2
模型势必导致很大的误差。而一些高级的湍流模型, 如低雷诺
数K ε2模型、代数应力模型、雷诺应力模型等, 可以获得较好的精度, 但是其计算量已相当于复杂二阶封闭模型或大涡模拟的计算量, 无法为工程应用所接受[11]。
实际应用中往往需要对各种不同空调通风气流组织工况进行对比以优化设计, 因此希望在规划设计阶段就能对这些工况进行分析评价以作出选择。这就要求CFD 程序能快速模拟室内空气分布情况, 而采用复杂模型目前还难以满足这一要求。于是, 在现阶段, 为快速对通风空调房间气流组织进行数值模拟, 必须采用较简单的湍流模型, 但又
应该能保证工程上的准确度。
2. 2 对策或解决思路
近年来, 随着湍流实验和模拟技术的发展, 一些学者根据实验数据或直接数值模拟DNS (directly numerical simulation ) 的计算结果, 提出满足工程应用要求的快速、简单的零方程模型, 针对性地对所关心的问题进行模拟。尽管零方程模型比复杂模型简单, 但是在专门的领域内却能获得比复杂模型更与实验数据吻合的结果。对于暖通空调领域的室内空气流动模拟而言,Xu 等人在M IT 对室内空气自然对流和混合对流流动研究基础上提出一方程模型和零方程模型, 取得了令人满意的准确度, 且计算速度很快[11]。作为解决以上问题的一个对策, 这里简要介绍笔者对M IT 零方程湍流模型的应用情况。
M IT 零方程湍流模型是在室内空气自然对流
和混合对流的DNS 结果基础上提出的[11], 该模型针对房间内非等温流动的瑞利数Ra 范围((2. 6～3. 0) ×1010) , 认为涡黏系数正比于流体密度、当地速度和距壁面最近之距离, 比例系数由DNS 的结果拟合而得。由于该模型少求解两个微分方程, 而仅求解关于质量、动量和能量守恒的5个微分方程, 故计算最省时间。下面举一较为复杂的非等温算例。
H. B. Awbi 为了验证房间通风情况数值计算的结果, 于1989年作了如下实验[13]。在一间长×宽×高为4. 2m ×4. 2m ×2. 8m 的小屋内非等温送风, 采用顶送风, 送风口为宽0. 024m 的条缝风口, 回风口位于小屋右下角, 高为0. 05m 。斜向下45°送风, 送风量为0. 06m 3/s , 送风温度为15. 5℃, 回风温度为22℃, 送、回风温差为6. 5℃, 室内由电加热器模拟均匀分布热负荷19W/m 2。房间结构如图3所示
。
图3 实验及计算用房间示意图
通常人们关心的只是工作区的温度和速度值,
故实验所测点为房间正中(长度方向一半处, x =2. 1m ) 高度依次为0. 15,0. 6,1. 2,1. 8m 处的温度和速度值, 利用前述模型进行计算, 所得结果与实验值的比较如图4所示。各点速度差值均在
图4 实验与计算值对比
0108m/s 以下, 温度差也很小, 最大差值仅为0. 8
℃, 计算值和测量值吻合较好。在一台P Ⅲ500, 128M 内存的个人计算机上, 对于本次计算采用的
网格数(24×16×14) , 计算约10min 即获收敛。
更多的算例表明, 该零方程湍流模型对于室内空气流动的数值模拟有着速度快、能满足工程上精度要求的特点[14], 值得在工程应用中尝试。而简单的、适用性更广的湍流模型也需要进一步研究, 这可借助更具广泛性的直接数值模拟结果或实验数据来进行。
3 计算收敛速度
3. 1 计算速度对实际工程应用的意义和存在的问
题
暖通空调中的实际流动和传热问题比较复杂, 可能需要较多的网格才能反映真实的物理特点, 从而带来计算耗时太长的问题。而计算速度快对于工程应用而言, 具有重大的意义。由前所述, 实际应用中往往需要比较分析多个工况以得出最优方案, 这就要求数值计算能较快获得模拟结果, 以保证空调气流组织规划设计阶段的周期满足工程要求。前面的讲述表明, 这可从风口模型、湍流模型等方面作出改进。另外, 作为通风空调房间气流组织数值模拟主要组成部分的代数方程求解算法, 对计算速度有着很大的影响。近年来, 也有研究者对此作了研究, 如文献[15]中提出利用多重网格法进
行室内空气流动的数值模拟。该文所举算例表明,
采用多重网格法可以将迭代计算收敛时间缩短1/3左右。但是该文献只对一个比较简单的等温流动算例进行了验证, 对于空调工程中常见的非等温流动以及内部带有家具和人员等复杂流动的算例没有说明。最主要的是, 由于目前暖通空调领域内的数值模拟通常采用有限容积法离散, 并用SIMPL E 算法求解, 具有明显的物理意义, 但不存在多重网格法所要求的数学上的迭代残差, 故多重网格法并不适用。3. 2 对策
鉴于上述问题, 笔者提出了比多重网格法更简单易行、更适合SIMPL E 算法的代数方程求解方法———误差预处理法, 以获得加快迭代计算收敛速度的效果。误差预处理法以多重网格法的高频和
低频误差理论为基础, 在粗网格上预处理低频误
差, 再在细网格上进一步消除高频误差, 由于较难消除的低频误差在粗网格上只需很短的时间就可计算收敛而消除, 故节约了计算时间。该方法从控制体节点物理量与周围物理量之间的关系出发, 用面积率表示不同粗细网格之间节点值的关系, 具有明显的物理意义, 且不会改变原来方程的守恒特性[16]。仍以2. 2节中图3所示算例进行比较。误差预处理法只是针对代数方程求解算法作出修改, 因此求解结果变化不大, 与图4所示类似, 和实验值吻合得很好, 这里不再列举其结果。而计算时间和传统迭代方法的时间比较见表1(所用计算机为PIII500,128M 内存的个人计算机) 。
表1 误差预处理法和传统迭代法计算时间比较
网格划分
计算至收敛时间/min
误差预处理法粗:18×18×16细:36×36×32314传统迭代法
36×36×32
501
文献[16]中对室内等温和非等温流动的模拟
算例都表明, 误差预处理法结合SIMPL E 算法进行通风空调房间内空气流动数值模拟, 可使迭代计算的收敛速度加快1/3～1/2左右。因此, 误差预处理法可满足实际工程快速计算的要求。4 辐射和内热源模型4. 1 辐射换热和内热源对室内空气分布的重要性和数值模拟存在的问题
CFD 求解的是控制流体流动和传热的N 2S 方
程, 解决典型的对流扩散问题, 但并没有考虑壁面
之间的相互辐射。通常空调房间内各壁面的温度相
差不大, 通过壁面之间的相互辐射, 改变壁面温度, 壁面最终再通过对流换热的形式对室内空气分布的影响是很小的, 所以在实际工程中这种考虑方法可行。但是对于采用辐射吊顶、地板供暖等空调方式以及散热器供暖的空间, 不同壁面之间温度相差较大, 不考虑壁面之间辐射换热将无法准确给出影响室内空气分布的边界条件。目前比较常见的是实测壁面温度的值, 利用辐射网络法计算各个壁面之间的辐射换热量, 提供CFD 计算的边界条件。但是对于需要模拟预测的工况, 往往其暖通空调系统还没有建设完成, 无法事先进行测量。
另外, 数值模拟气流组织时, 一般认为内部热源的所有发热量以对流换热的形式散失到室内, 这主要是基于空调负荷计算的思路考虑的。但是, 与空调负荷计算的集总参数考虑方法不同, 利用CFD 对室内空气场分布进行数值模拟预测的目标是得到气流组织的场分布特性, 此时热源的对流散热和辐射散热对室内空气流动的影响是不同的:对流散热是通过微元体之间的温差, 形成空气微团之间的热交换和质交换, 从而带动室内的空气流动; 而辐射散热是热源通过辐射传热影响壁面温度, 再以壁面与室内空气对流换热的形式影响室内的空气场的流动。进一步而言, 热源辐射散热量的方向性对室内空气分布也有影响, 不同热源在不同方向的辐射散热量也是有所区别的。
综上可见, 辐射换热和内热源传热对室内空气流动情况都有着重要的影响, 尤其对于辐射换热比例大的用散热器供暖、以辐射吊顶或地板供暖空调的房间, 以及内热源比例很大的空间, 如大礼堂、影剧院、工厂车间等更是如此。而目前对通风空调房间气流组织的数值计算还没有或很少考虑这些问题。4. 2 解决思路
对于辐射换热问题, 可以通过对壁面的各项传热量, 如与空气的对流换热量、墙壁之间的辐射换热量、外墙传热量等, 依据能量守恒建立辐射计算模型, 据此计算墙壁温度或热流, 作为CFD 数值计算的边界条件, 可根据温度的不均匀性将墙壁划分为数个网格进行数值求解。而CFD 数值计算的结果对墙壁换热量又有影响, 辐射换热计算和CFD 计算是互相耦合的, 因此可迭代计算, 从而避免需
要实验测定边壁温度作为边界条件; 对于热源换热问题, 收集、调研各种热源的对流、辐射发热量比例, 研究热源传热机理, 建立常见热源的热源模型, 如对流辐射散热比例、辐射散热的方向性等, 在模拟计算中再分别处理两种热源形式。比如, 笔者采用这种方式对某灯光、人员密集的大会堂的气流组织进行模拟以进行改造, 取得了较为满意的结果[17]。由于内热源传热含有辐射换热, 因此在辐射模型中也应该考虑之, 再与CFD 方法耦合计算。图5
所示的算法可将辐射换热和内热源传热的影
图5 辐射换热、热源传热与CFD 的联合算法
响全面地考虑到气流组织的数值计算中, 可作为解
决室内空气流动数值模拟的辐射换热和内热源传热影响问题的主要思路。
但是, 上述方法还存在很多难点, 比如辐射模型中墙壁对流换热量的表面传热系数的确定、耦合计算带来的计算时间过长等问题, 还需要进一步的研究。5 复杂几何物理条件的描述5. 1 实际工程问题的复杂性
实际工程中的气流组织数值计算除了存在以上特殊性或难点之外, 往往还会遇到计算区域或物体几何形状不规则的情况, 这种复杂的几何条件很难直接用通常采用的直角坐标或柱坐标系进行描述, 目前的做法多是在上述坐标系下用阶梯网格逼近复杂的几何形状, 但带来计算量大、计算结果并不令人满意的问题; 另外, 对于体育场馆、大会堂等大空间, 人员、灯光以及座椅设备过多, 难以一一模拟, 而它们作为室内的热源或阻挡物, 对气流组织有着重要的影响。目前采取的对策多将其作为整体处理, 即按局部的集总参数法考虑, 以保证室内工作区定性的空气分布特性, 但这样难以反映它们的传热传质特性, 无法正确模拟各个人员、设备之间的微气候。
5. 2 解决思路或对策
对于几何形状复杂的区域或物体, 目前可用贴
体坐标、非结构化网格等方法进行离散和描述, 但是随之而来的网格生成、计算时间、收敛问题等仍是影响实际工程应用的难点。近来, 有人提出用有限元方法对复杂形状的房间及其内部的物体进行模拟[18], 并且与面向对象的方法结合起来, 取得了一定的成果, 值得进一步推广使用; 而对于具有众多人员、灯光、设备的空间的物理条件描述问题, 仍没有有效的解决办法。目前有一种类比多孔介质研究方法的思路值得借鉴[19], 但是有关的传热、传质特性尚待进一步研究, 离实际的工程应用还有一定的距离。6 小结
风口模型、湍流模型、计算收敛速度、热源和辐射模型以及复杂几何物理条件的描述都是暖通空调气流组织数值模拟的特殊问题。针对以上问题, 笔者分别提出了N 点风口模型, 采用M IT 零方程湍流模型、误差预处理法等方法, 通过有关算例表明这些方法都是颇有成效的; 也提出了辐射换热、热源传热和CFD 耦合计算的思路来考虑辐射和热源对气流组织的影响, 但是诸如墙壁表面传热系数的取值等问题还需要研究; 同时, 也指出了面向对象有限元方法、多重介质类比方法等可能有效解决实际问题复杂几何物理条件的途径, 作为进一步研究的思路。
总之, 尽管在暖通空调气流组织数值计算中还存在诸多困难, 但人们正一步步向利用CFD 技术解决复杂工程问题这个目标迈进, 随着众多学者的努力, 我们相信, 对暖通空调工程复杂的气流组织进行快速、准确的数值计算, 不会是高不可攀的梦想, 这个堡垒必将被征服, 让我们和所有同仁一道, 继续不懈地努力。参考文献
1 余常昭. 紊动射流. 北京:高等教育出版社,19932 Nielsen P V , Restivo A , Whitelaw J H. The velocity
characteristics of ventilated rooms.
Trans ASME J of Fluids Engineering , 1978,100(9) :291
298
3 Nielsen P V. Description of supply openings in numerical
models for room air distribution. In :ASHRAE Trans. 1992,98. 963
971
4 G osman A D , Nielsen P V , Restiv o A , et al. The flow
properties in rooms with small ventilation openings. Trans ASME J of Fluids Engineering , 1980,102(3) :321
323
5 Huo Y , Haghighat F , Zhang J S , et al. A systematic
approach to describe the air terminal device in CFD simulation for room air distribution analysis. Building and Environment , 2000, 35(6) :563
576
6 Heikkinen J. Modeling of a supply air terminal for room
air flow simulation. In :Proceedings of the 12th AIVC Conference on Air Movement and Ventilation Control within Buildings , Ottawa , Canada , 1991. 213230
7 Chen Qingyan , Moser Alfred. Simulation of a multi ple 2
nozzle diffuser.
In :Proceedings of the 12th AIVC
Conference on Air Movement and Ventilation Control within Buildings , Ottawa , Canada , 1991. 113
8 赵彬. 室内空气流动数值模拟的风口模型研究及应用.
[博士学位论文].北京:清华大学,2001
9 赵彬, 李先庭, 彦启森. 室内空气流动数值模拟的风口模型综述. 暖通空调, 2000,30(5) :33
37
10Chen Qingyan , Moser A , Huber A.
Prediction of buoyant , turbulent flow by a low 2Renolds 2number k ε2
model. In :ASHRAE Trans. 1990,96. 564
573
11Xu Weiran. New turbulence models for indoor air flow
simulation :[Ph D thesis].Massachusetts :MassachusettsInstitute of Technology , 199812Nielsen P V.
The selection of turbulence models for
prediction of room airflow , In :ASHRAE Trans. 1998,
104. 1119
1126
13Awbi H B. Application of computational fluid dynamics in
room ventilation. Building and Environment , 1989, 24(1) :73
84
14赵彬, 李先庭, 彦启森. 用零方程湍流模型模拟通风空
调室内空气流动. 清华大学学报(自然科学版) ,2001,
41(10) :109
113
15邓启红, 汤广发. 多重网格法在通风空调气流数值模拟
中的应用. 应用基础与工程科学学报,1998,6(1) :48
54
16赵彬, 李先庭, 彦启森. 室内空气流动数值模拟的误差
预处理法. 清华大学学报(自然科学版) ,2001,41(10) :
114
117
17赵彬. 某国际会展中心大会堂气流组织计算分析. 见:
清华同方人环工程公司报告. 1999
18邓启红, 汤广发. 面向对象有限元分析———求解暖通领
域复杂传热与流动问题的有效途径. 见:全国暖通空调制冷2000年学术文集. 248
254
19Li Yuguo. CFD in building service engineering. 见:清华
大学建筑学院建筑技术科学系交流报告. 1999
暖通空调气流组织数值模拟的特殊性
清华大学 赵 彬☆ 李先庭 彦启森
摘要 从计算流体力学CFD 应用于暖通空调工程的角度, 阐述了通风空调房间内气流组织数值计算的特殊性和难点, 对风口模型、湍流模型、辐射模型、热源模型、计算收敛速度以及复杂物理条件的描述等作了分析和介绍, 并提出了相应的对策或解决思路。
关键词 计算流体力学 通风空调房间 气流组织 数值模拟
Pa rti c ul a rit y of n u m e ri c a l si m ul a ti o n f or a ir distri b uti o n
i n v e ntil a t e d a n d a ir c o n diti o n e d r o o m s
By Zhao Bin ★, Li X ianting and Y an Qisen
Abs t r a ct Elaborates t he p articularit y and difficulties of t he numerical simulation f or air dist ribution by means of comp utational fluid dynamics (C FD ) applied in HVA C engineering. A nalyses air op ening model , turbulence model , radiation model and heat source model , and p resents t he calculation of convergence sp eed and t he descrip tion of complicated p hysical conditions. Puts f orward corresp onding counter measures and solutions f or t he difficulties.
Keywor ds comp utational fluid dynamics , ventilated and air conditioned room , air dist ribution , numerical simulation
★Tsinghua University , Beijing , China
0 引言
调领域数值计算的特殊性。为将CFD 技术这一强有力的工具应用于暖通空调气流组织设计, 还需要对这些问题进行研究。下面将结合笔者和其他研究者的研究成果和经验, 逐一介绍之, 以期得到同行的指正。1 风口模型
1. 1 风口模型的重要性和不足
自1974年丹麦学者P. V. Nilsen 首次将计算流体力学CFD 技术应用于暖通空调室内空气流动的数值模拟以来, 众多学者在该领域付出了辛勤的劳动, 而大家只有一个奋斗目标:利用数值模拟技术快速、准确地解决暖通空调实际工程中的流动和传热问题。随着计算机技术、数值计算技术以及湍流理论的发展, 人们在暖通空调工程的CFD 应用研究领域取得了长足的进步。在短短的20多年内, 涌现出了各式各样关于空调通风房间室内空气流动数值模拟的算法、湍流模型以及对一些实际应用中的特殊问题的解决办法, 暖通空调的数值计算正在蓬勃发展当中。但是,CFD 技术本身还不完善, 比如人们对湍流现象的认识尚未完全清楚; 而且, 暖通空调实际工程问题比较复杂, 如送风口的入流边界条件描述, 辐射换热和热源的考虑, 实际人员、设备等边界条件的描述等; 另外, 在实际应用中还存在要求计算速度快、能在较短时间内模拟大量工况等问题。以上各种问题共同构成了暖通空
风口入流边界条件的描述是当前通风空调房间气流组织数值模拟的难题之一。从数学角度而言, 对于CFD 所求解的描述室内空气湍流输运的非线性偏微分方程组, 边界条件对计算结果的准确性有着极大的影响; 从物理角度而言, 空调通风房间的空气流动就是由风口射流引起的, 射流动量流量和质量流量对室内空气分布情况起关键作用。
☆赵彬, 男,1974年12月生, 博士研究生, 工学博士, 讲师
100084北京清华大学建筑学院建筑技术科学系
(010) 62779995
E 2mail :[email protected]. cn 收稿日期:20010103修回日期:2004
09
因此建立风口模型来合理描述房间入流边界条件
是准确、快速地模拟通风空调房间气流组织的关键。不同于CFD 技术在其他领域的应用, 暖通空调工程领域中室内空气流动数值模拟的风口入流边界条件是很复杂的。为了正确预测室内空气分布情况, 我们需知道送风参数的详细情况以正确描述数值模拟时的入流边界条件。而实际的送风口几何形状却很复杂, 如散流器、百叶风口和孔板风口等, 且送风射流雷诺数通常大于30, 为湍动射流[1]。这样的入流条件极其复杂, 难以简单描述。为了了解这些风口的送风参数, 在用CFD 方法进行数值模拟时计算区域内的网格必须划分至mm 甚至更小的量级; 而为了模拟室内空气流动情况, 反映室内湍流大涡旋对流动的影响, 计算区域必须达到m 的量级。这样, 对于实际应用中的三维问题而言, 必将导致计算区域内的网格节点数目巨大, 超过目前一般计算机的容量, 从而使得计算量无法为一般的工程应用所接受。因此, 在对通风空调房间进行室内空气流动的数值模拟时, 必须采用较少的计算网格来描述复杂的风口入流边界条件以适应工程应用中快速计算的要求, 同时又不失一定的准确度, 于是, 对风口模型的研究就成为CFD 技术应用于暖通空调领域内的独特之处, 同时也是暖通空调气流组织数值模拟需要研究的一个非常重要的问题。
自20世纪70年代起, 就有学者提出用风口模型来解决室内空气流动数值模拟的送风口边界条件描述问题[2]。经过20多年的发展, 逐渐形成了两类风口模型:间接描述类风口模型和直接描述类风口模型。间接描述类风口模型包括盒子模型[3]、指定速度模型[3～4]和主流区模型[5]等, 这类模型把送风口入流边界条件描述转化到包围送风口的一个盒子的边界上, 借助实验值或者射流公式描述边界条件, 避免了直接描述送风口出口处复杂的流动情况, 但这种方法需要实验测量, 非常麻烦, 且无法用于设计阶段的预测; 另外, 还受射流公式不可靠、受射流是否受限或等温的影响。直接描述类风口模型包括基本模型[6]、动量模型[7], 这类模型直接在送风口处描述入流边界条件, 具有不受送风射流等温与否、受限与否的限制等优点, 但是基本模型对于开口有效面积系数(开口有效面积与总面积之比) 较小的送风口误差较大; 动量模型使用方便, 但对于如何简化合
并不同出流方向的射流并不明确。1. 2 对策或解决思路
针对现有风口模型的不足, 笔者提出了N 点风口模型以适应工程中快速简捷模拟风口入流边界条件的需要。N 点风口模型在保证能准确模拟入流质量流量、动量流量以及浮力通量的前提下, 用N 个简单开口替代不同出流方向的送风口以实现其入流边界条件的简化描述, 属于直接描述的风口模型, N 的取值根据风口射流动量数进行简化合并, 有关详细情况可参见文献[8]。该风口模型对百叶风口常见不同出流条件、孔板这类有效面积小的风口都能获得满意的模拟结果, 下面举例说明之。
图1为某百叶送风口出流方向张角为60°
的
图1 百叶送风口60°张角出流的N 点模型(N =3)
N 点风口模型, 这是一种百叶风口比较常见的出
流条件。图2a
为此时该百叶风口在出流流速为
图2 百叶送风口60°张角出流下射流主体段
特性模拟值和实验值对比
2, 5, 8, 10m/s 时的主体段轴心速度衰减和实验
值的对比, 可见N 点风口模型模拟结果和实验值
吻合得很好; 图2b 为四种流速下的N 点风口模型模拟所得主体段断面量纲一的速度分布与实验结果的对比, 可见模拟所得结果符合送风口湍流射流主体段断面流速自相似分布的特点。图中v m 为射流轴心速度, v o 为射流出口速度, r 为断面上某考察点距轴心的径向距离, r 0. 5为断面上流速为轴心速度0. 5倍之点距轴心的径向距离。尽管N 点风口模型能解决常见的百叶送风口和多孔板送风口的出流边界条件描述问题, 比较适用于工程应用, 但是对于更为复杂的送风口边界条件尚存在较大误差, 需要寻找更通用的风口模型, 或者新的室内空气流动的模拟思路, 如分段模拟室内空气流动的思路[9]。2 湍流模型
2. 1 气流组织数值计算湍流模型的特点及不足
由于通风空调房间内空气流动多为湍流流动, 而人们对湍流机理的认识至今还不完全, 对于工程应用, 目前多为半经验的、唯象的模拟, 故此对湍流的模拟始终成为室内空气流动数值模拟的一个主要研究方向。限于目前的计算机能力, 工程中最常采用的是涡黏系数模型EVM (eddy viscosity
model ) 中的K ε2
两方程模型或其变形。而K 2ε模型对于等温流动情形能模拟得很好, 但是对于空调通风房间内的非等温、混合对流的情形却有较大误差[10～11]。Nielsen 于1998年指出, 对于热羽流、贴壁射流、温度分层流动等需要不同的湍流模型进行模拟方能取得满意的精度[12], 而这些流动形式往
往并存于实际的空调房间, 采用K ε2
模型势必导致很大的误差。而一些高级的湍流模型, 如低雷诺
数K ε2模型、代数应力模型、雷诺应力模型等, 可以获得较好的精度, 但是其计算量已相当于复杂二阶封闭模型或大涡模拟的计算量, 无法为工程应用所接受[11]。
实际应用中往往需要对各种不同空调通风气流组织工况进行对比以优化设计, 因此希望在规划设计阶段就能对这些工况进行分析评价以作出选择。这就要求CFD 程序能快速模拟室内空气分布情况, 而采用复杂模型目前还难以满足这一要求。于是, 在现阶段, 为快速对通风空调房间气流组织进行数值模拟, 必须采用较简单的湍流模型, 但又
应该能保证工程上的准确度。
2. 2 对策或解决思路
近年来, 随着湍流实验和模拟技术的发展, 一些学者根据实验数据或直接数值模拟DNS (directly numerical simulation ) 的计算结果, 提出满足工程应用要求的快速、简单的零方程模型, 针对性地对所关心的问题进行模拟。尽管零方程模型比复杂模型简单, 但是在专门的领域内却能获得比复杂模型更与实验数据吻合的结果。对于暖通空调领域的室内空气流动模拟而言,Xu 等人在M IT 对室内空气自然对流和混合对流流动研究基础上提出一方程模型和零方程模型, 取得了令人满意的准确度, 且计算速度很快[11]。作为解决以上问题的一个对策, 这里简要介绍笔者对M IT 零方程湍流模型的应用情况。
M IT 零方程湍流模型是在室内空气自然对流
和混合对流的DNS 结果基础上提出的[11], 该模型针对房间内非等温流动的瑞利数Ra 范围((2. 6～3. 0) ×1010) , 认为涡黏系数正比于流体密度、当地速度和距壁面最近之距离, 比例系数由DNS 的结果拟合而得。由于该模型少求解两个微分方程, 而仅求解关于质量、动量和能量守恒的5个微分方程, 故计算最省时间。下面举一较为复杂的非等温算例。
H. B. Awbi 为了验证房间通风情况数值计算的结果, 于1989年作了如下实验[13]。在一间长×宽×高为4. 2m ×4. 2m ×2. 8m 的小屋内非等温送风, 采用顶送风, 送风口为宽0. 024m 的条缝风口, 回风口位于小屋右下角, 高为0. 05m 。斜向下45°送风, 送风量为0. 06m 3/s , 送风温度为15. 5℃, 回风温度为22℃, 送、回风温差为6. 5℃, 室内由电加热器模拟均匀分布热负荷19W/m 2。房间结构如图3所示
。
图3 实验及计算用房间示意图
通常人们关心的只是工作区的温度和速度值,
故实验所测点为房间正中(长度方向一半处, x =2. 1m ) 高度依次为0. 15,0. 6,1. 2,1. 8m 处的温度和速度值, 利用前述模型进行计算, 所得结果与实验值的比较如图4所示。各点速度差值均在
图4 实验与计算值对比
0108m/s 以下, 温度差也很小, 最大差值仅为0. 8
℃, 计算值和测量值吻合较好。在一台P Ⅲ500, 128M 内存的个人计算机上, 对于本次计算采用的
网格数(24×16×14) , 计算约10min 即获收敛。
更多的算例表明, 该零方程湍流模型对于室内空气流动的数值模拟有着速度快、能满足工程上精度要求的特点[14], 值得在工程应用中尝试。而简单的、适用性更广的湍流模型也需要进一步研究, 这可借助更具广泛性的直接数值模拟结果或实验数据来进行。
3 计算收敛速度
3. 1 计算速度对实际工程应用的意义和存在的问
题
暖通空调中的实际流动和传热问题比较复杂, 可能需要较多的网格才能反映真实的物理特点, 从而带来计算耗时太长的问题。而计算速度快对于工程应用而言, 具有重大的意义。由前所述, 实际应用中往往需要比较分析多个工况以得出最优方案, 这就要求数值计算能较快获得模拟结果, 以保证空调气流组织规划设计阶段的周期满足工程要求。前面的讲述表明, 这可从风口模型、湍流模型等方面作出改进。另外, 作为通风空调房间气流组织数值模拟主要组成部分的代数方程求解算法, 对计算速度有着很大的影响。近年来, 也有研究者对此作了研究, 如文献[15]中提出利用多重网格法进
行室内空气流动的数值模拟。该文所举算例表明,
采用多重网格法可以将迭代计算收敛时间缩短1/3左右。但是该文献只对一个比较简单的等温流动算例进行了验证, 对于空调工程中常见的非等温流动以及内部带有家具和人员等复杂流动的算例没有说明。最主要的是, 由于目前暖通空调领域内的数值模拟通常采用有限容积法离散, 并用SIMPL E 算法求解, 具有明显的物理意义, 但不存在多重网格法所要求的数学上的迭代残差, 故多重网格法并不适用。3. 2 对策
鉴于上述问题, 笔者提出了比多重网格法更简单易行、更适合SIMPL E 算法的代数方程求解方法———误差预处理法, 以获得加快迭代计算收敛速度的效果。误差预处理法以多重网格法的高频和
低频误差理论为基础, 在粗网格上预处理低频误
差, 再在细网格上进一步消除高频误差, 由于较难消除的低频误差在粗网格上只需很短的时间就可计算收敛而消除, 故节约了计算时间。该方法从控制体节点物理量与周围物理量之间的关系出发, 用面积率表示不同粗细网格之间节点值的关系, 具有明显的物理意义, 且不会改变原来方程的守恒特性[16]。仍以2. 2节中图3所示算例进行比较。误差预处理法只是针对代数方程求解算法作出修改, 因此求解结果变化不大, 与图4所示类似, 和实验值吻合得很好, 这里不再列举其结果。而计算时间和传统迭代方法的时间比较见表1(所用计算机为PIII500,128M 内存的个人计算机) 。
表1 误差预处理法和传统迭代法计算时间比较
网格划分
计算至收敛时间/min
误差预处理法粗:18×18×16细:36×36×32314传统迭代法
36×36×32
501
文献[16]中对室内等温和非等温流动的模拟
算例都表明, 误差预处理法结合SIMPL E 算法进行通风空调房间内空气流动数值模拟, 可使迭代计算的收敛速度加快1/3～1/2左右。因此, 误差预处理法可满足实际工程快速计算的要求。4 辐射和内热源模型4. 1 辐射换热和内热源对室内空气分布的重要性和数值模拟存在的问题
CFD 求解的是控制流体流动和传热的N 2S 方
程, 解决典型的对流扩散问题, 但并没有考虑壁面
之间的相互辐射。通常空调房间内各壁面的温度相
差不大, 通过壁面之间的相互辐射, 改变壁面温度, 壁面最终再通过对流换热的形式对室内空气分布的影响是很小的, 所以在实际工程中这种考虑方法可行。但是对于采用辐射吊顶、地板供暖等空调方式以及散热器供暖的空间, 不同壁面之间温度相差较大, 不考虑壁面之间辐射换热将无法准确给出影响室内空气分布的边界条件。目前比较常见的是实测壁面温度的值, 利用辐射网络法计算各个壁面之间的辐射换热量, 提供CFD 计算的边界条件。但是对于需要模拟预测的工况, 往往其暖通空调系统还没有建设完成, 无法事先进行测量。
另外, 数值模拟气流组织时, 一般认为内部热源的所有发热量以对流换热的形式散失到室内, 这主要是基于空调负荷计算的思路考虑的。但是, 与空调负荷计算的集总参数考虑方法不同, 利用CFD 对室内空气场分布进行数值模拟预测的目标是得到气流组织的场分布特性, 此时热源的对流散热和辐射散热对室内空气流动的影响是不同的:对流散热是通过微元体之间的温差, 形成空气微团之间的热交换和质交换, 从而带动室内的空气流动; 而辐射散热是热源通过辐射传热影响壁面温度, 再以壁面与室内空气对流换热的形式影响室内的空气场的流动。进一步而言, 热源辐射散热量的方向性对室内空气分布也有影响, 不同热源在不同方向的辐射散热量也是有所区别的。
综上可见, 辐射换热和内热源传热对室内空气流动情况都有着重要的影响, 尤其对于辐射换热比例大的用散热器供暖、以辐射吊顶或地板供暖空调的房间, 以及内热源比例很大的空间, 如大礼堂、影剧院、工厂车间等更是如此。而目前对通风空调房间气流组织的数值计算还没有或很少考虑这些问题。4. 2 解决思路
对于辐射换热问题, 可以通过对壁面的各项传热量, 如与空气的对流换热量、墙壁之间的辐射换热量、外墙传热量等, 依据能量守恒建立辐射计算模型, 据此计算墙壁温度或热流, 作为CFD 数值计算的边界条件, 可根据温度的不均匀性将墙壁划分为数个网格进行数值求解。而CFD 数值计算的结果对墙壁换热量又有影响, 辐射换热计算和CFD 计算是互相耦合的, 因此可迭代计算, 从而避免需
要实验测定边壁温度作为边界条件; 对于热源换热问题, 收集、调研各种热源的对流、辐射发热量比例, 研究热源传热机理, 建立常见热源的热源模型, 如对流辐射散热比例、辐射散热的方向性等, 在模拟计算中再分别处理两种热源形式。比如, 笔者采用这种方式对某灯光、人员密集的大会堂的气流组织进行模拟以进行改造, 取得了较为满意的结果[17]。由于内热源传热含有辐射换热, 因此在辐射模型中也应该考虑之, 再与CFD 方法耦合计算。图5
所示的算法可将辐射换热和内热源传热的影
图5 辐射换热、热源传热与CFD 的联合算法
响全面地考虑到气流组织的数值计算中, 可作为解
决室内空气流动数值模拟的辐射换热和内热源传热影响问题的主要思路。
但是, 上述方法还存在很多难点, 比如辐射模型中墙壁对流换热量的表面传热系数的确定、耦合计算带来的计算时间过长等问题, 还需要进一步的研究。5 复杂几何物理条件的描述5. 1 实际工程问题的复杂性
实际工程中的气流组织数值计算除了存在以上特殊性或难点之外, 往往还会遇到计算区域或物体几何形状不规则的情况, 这种复杂的几何条件很难直接用通常采用的直角坐标或柱坐标系进行描述, 目前的做法多是在上述坐标系下用阶梯网格逼近复杂的几何形状, 但带来计算量大、计算结果并不令人满意的问题; 另外, 对于体育场馆、大会堂等大空间, 人员、灯光以及座椅设备过多, 难以一一模拟, 而它们作为室内的热源或阻挡物, 对气流组织有着重要的影响。目前采取的对策多将其作为整体处理, 即按局部的集总参数法考虑, 以保证室内工作区定性的空气分布特性, 但这样难以反映它们的传热传质特性, 无法正确模拟各个人员、设备之间的微气候。
5. 2 解决思路或对策
对于几何形状复杂的区域或物体, 目前可用贴
体坐标、非结构化网格等方法进行离散和描述, 但是随之而来的网格生成、计算时间、收敛问题等仍是影响实际工程应用的难点。近来, 有人提出用有限元方法对复杂形状的房间及其内部的物体进行模拟[18], 并且与面向对象的方法结合起来, 取得了一定的成果, 值得进一步推广使用; 而对于具有众多人员、灯光、设备的空间的物理条件描述问题, 仍没有有效的解决办法。目前有一种类比多孔介质研究方法的思路值得借鉴[19], 但是有关的传热、传质特性尚待进一步研究, 离实际的工程应用还有一定的距离。6 小结
风口模型、湍流模型、计算收敛速度、热源和辐射模型以及复杂几何物理条件的描述都是暖通空调气流组织数值模拟的特殊问题。针对以上问题, 笔者分别提出了N 点风口模型, 采用M IT 零方程湍流模型、误差预处理法等方法, 通过有关算例表明这些方法都是颇有成效的; 也提出了辐射换热、热源传热和CFD 耦合计算的思路来考虑辐射和热源对气流组织的影响, 但是诸如墙壁表面传热系数的取值等问题还需要研究; 同时, 也指出了面向对象有限元方法、多重介质类比方法等可能有效解决实际问题复杂几何物理条件的途径, 作为进一步研究的思路。
总之, 尽管在暖通空调气流组织数值计算中还存在诸多困难, 但人们正一步步向利用CFD 技术解决复杂工程问题这个目标迈进, 随着众多学者的努力, 我们相信, 对暖通空调工程复杂的气流组织进行快速、准确的数值计算, 不会是高不可攀的梦想, 这个堡垒必将被征服, 让我们和所有同仁一道, 继续不懈地努力。参考文献
1 余常昭. 紊动射流. 北京:高等教育出版社,19932 Nielsen P V , Restivo A , Whitelaw J H. The velocity
characteristics of ventilated rooms.
Trans ASME J of Fluids Engineering , 1978,100(9) :291
298
3 Nielsen P V. Description of supply openings in numerical
models for room air distribution. In :ASHRAE Trans. 1992,98. 963
971
4 G osman A D , Nielsen P V , Restiv o A , et al. The flow
properties in rooms with small ventilation openings. Trans ASME J of Fluids Engineering , 1980,102(3) :321
323
5 Huo Y , Haghighat F , Zhang J S , et al. A systematic
approach to describe the air terminal device in CFD simulation for room air distribution analysis. Building and Environment , 2000, 35(6) :563
576
6 Heikkinen J. Modeling of a supply air terminal for room
air flow simulation. In :Proceedings of the 12th AIVC Conference on Air Movement and Ventilation Control within Buildings , Ottawa , Canada , 1991. 213230
7 Chen Qingyan , Moser Alfred. Simulation of a multi ple 2
nozzle diffuser.
In :Proceedings of the 12th AIVC
Conference on Air Movement and Ventilation Control within Buildings , Ottawa , Canada , 1991. 113
8 赵彬. 室内空气流动数值模拟的风口模型研究及应用.
[博士学位论文].北京:清华大学,2001
9 赵彬, 李先庭, 彦启森. 室内空气流动数值模拟的风口模型综述. 暖通空调, 2000,30(5) :33
37
10Chen Qingyan , Moser A , Huber A.
Prediction of buoyant , turbulent flow by a low 2Renolds 2number k ε2
model. In :ASHRAE Trans. 1990,96. 564
573
11Xu Weiran. New turbulence models for indoor air flow
simulation :[Ph D thesis].Massachusetts :MassachusettsInstitute of Technology , 199812Nielsen P V.
The selection of turbulence models for
prediction of room airflow , In :ASHRAE Trans. 1998,
104. 1119
1126
13Awbi H B. Application of computational fluid dynamics in
room ventilation. Building and Environment , 1989, 24(1) :73
84
14赵彬, 李先庭, 彦启森. 用零方程湍流模型模拟通风空
调室内空气流动. 清华大学学报(自然科学版) ,2001,
41(10) :109
113
15邓启红, 汤广发. 多重网格法在通风空调气流数值模拟
中的应用. 应用基础与工程科学学报,1998,6(1) :48
54
16赵彬, 李先庭, 彦启森. 室内空气流动数值模拟的误差
预处理法. 清华大学学报(自然科学版) ,2001,41(10) :
114
117
17赵彬. 某国际会展中心大会堂气流组织计算分析. 见:
清华同方人环工程公司报告. 1999
18邓启红, 汤广发. 面向对象有限元分析———求解暖通领
域复杂传热与流动问题的有效途径. 见:全国暖通空调制冷2000年学术文集. 248
254
19Li Yuguo. CFD in building service engineering. 见:清华
大学建筑学院建筑技术科学系交流报告. 1999