房间内强自然对流的两方程模型模拟研究

高乃平，董 昆，朱 彤，向 阳

（1.同济大学 机械工程学院，上海201804；2.北京西门子西伯乐斯电子有限公司，北京100094）

房间内的空气流动是含有射流、绕流、汇流等多种流态的复杂、不稳定流动.通风房间内的气流场数值模拟自1976年丹麦Aalborg大学Nielsen教授的博士论文开始，已经开展了近35年的研究.大量的模拟结果表明：工程上常用的Standard（标准）k－ε模型和RNG（重整化）k－ε模型基本可以准确预测室内的速度场、温度场和污染物浓度场［1］.

然而室内发生火灾时，火源处的强热量释放速率导致室内湍流流动的一个显著特点是温度变化大，密度变化相应也很大.而且，火源驱动的强自然对流从本质上来说是一种非定常流动［2］，但在进行建筑物内火灾烟气迁移规律的数值模拟分析时，为了节省计算量和计算时间，常采用稳态模拟的方法.稳态模拟能否真实反映强自然对流的时间平均值需要深入研究.

在现有的计算能力下，强热源自然对流流场模拟常用的湍流模型是k－ε模型和大涡模拟［3］.两方程k－ε模型是通过求解附加的湍流脉动动能（k）方程和其耗散率（ε）方程，计算湍流黏性μt进而使得系综平均的Navier－Stokes控制方程封闭；大涡模拟则是利用滤波函数，把湍流运动分解成所谓的大涡和小涡，大涡直接求解而小涡通过亚格子模型进行模化求解.大涡模拟在网格划分时需要直接求解约80%的涡旋能量，这就限制了其网格尺度不能太大.同时，时间步长也应远远小于大涡的特征时间.因此，虽然大涡模拟能够提供诸如湍流脉动速度等瞬时信息，但通常计算同一个工况时其计算量约比k－ε模型多一个数量级或者更多［4］，这在一定程度上限制了其在火灾工程领域的广泛应用.

采用两方程k－ε模型研究强热源上方的热羽流动特性在文献中已有不少报道.文献［5］采用Standardk－ε模型模拟410～2 240kW的池火焰热羽流，并与经验公式计算所得热羽流体积流量进行对比，发现不同的经验公式结果相差较大，模拟值与部分经验公式计算结果吻合度较好.文献［6］认为k－ε会过高估计热羽流中心轴线的速度和温度，最终导致模拟得出的热羽形状比实际显得更为粗大.Kato等［7］发现模拟得到的浮升力大于实验测量值.然而，k－ε模型对整个房间内的流场预测效果的分析还比较少.

本文采用3种两方程k－ε模型，对房间内强自然对流流场进行稳态数值模拟，对比模拟结果与激光多普勒风速仪（LDV）和12μm的热电偶所测得的速度和温度分布，旨在分析3种常用的k－ε模型能否准确模拟室内大温度差和密度差的强自然对流特征，为其在火灾工程领域的应用提供基础.

1 模型小室实验

1.1 小室几何模型

实验的房间模型如图1所示［8］，与国际标准化组织（ISO：Room fire test）所用的火灾房间十分相似，但尺寸缩小为长1.8m×宽1.2m×高1.2m.

图1 火灾实验房间模型 （单位：mm）Fig.1 Model of fire test room （unit：mm）

四周墙壁采取绝热措施，通过地面、侧墙和屋顶散失的热量小于总发热量的2%.考虑到真实火源放热及火焰形态不稳定，并且会产生很多化学物质（燃烧产物），因此实验中采用电加热板来代替火源，以得到精确、均匀且稳定的热量输出.电加热板的表面辐射发射率为0.27.

1.2 实验工况

加热板输入功率为12kW·m－2，分别测量了3种工况：工况一，热源位于房间地面正中，即图1中的面①，总热量输入为1 080W；工况二，热源位于门对面的地板和墙壁正中，即图1中面②和③，总热量输入为5 400W；工况3，热源位于房间角落处，即图1中面④，⑤和⑥，总热量输入为9 720W.

1.3 测量方法和测点布置

速度测量采用多普勒激光测速仪，测量时间长度为1～3min，采样数据为1 000～6 000个.温度测量采用探头直径为12μm的镍铝热电偶，采样频率为200Hz.每个工况重复测量3次取平均值.

测点的布置如图2所示，沿X方向布置7个点（AA′线），沿Z方向布置5个点（BB′线）.其中AA′线上的第4点和BB′线上的第2点重合.每个点代表一组（x，z）坐标，变化高度（即y值），测得该点不同高度上的速度值和温度值.

图2 测量点布置 （单位：mm）Fig.2 Measuring point arrangement（unit：mm）

2 数值模拟方法

2.1 湍流模型

本文考察的两方程模型是Standardk－ε模型，RNGk－ε模型和Realizable（可实现）k－ε模型.Standardk－ε模型的控制方程为

式中：ρ为密度；t为时间；u为速度矢量；Γφ为变量φ的扩散系数；Sφ为变量φ的源项.表1给出φ取不同值时，公式（1）分别表达质量守恒、动量守恒、能量守恒方程和k与ε的传输方程.

表1 Standard k－ε模型的控制方程［9］Tab.1 Governing equations of standard k－εmodel

表1中，P表示压力；u，v和w分别为x，y和z方向的速度；μ，μt和μeff分别为动力黏性系数、湍流黏性系数和有效黏性系数；Pr为普朗特数；σk和σε分别为k和ε的湍流普朗特数；σT为湍流普朗特数；为湍流应力或速度梯度引起的k产生项，表示由浮升力引起的k产生项.

湍流黏性系数μt计算公式如下：

RNGk－ε模型中ε方程的源项多了一项Rε，其源项为项可更好地模拟强旋流动，计算方法为

式中：η＝Sk／ε；η0＝4.38；β＝0.012.其他各常数的取值为：Cμ＝0.084 5，C1ε＝1.42，C2ε＝1.68.

Realizablek－ε模型中ε方程的源项为ρC1Sε－其他各常数的取值为：C1ε＝1.44，C2＝1.9，σk＝1.0，σε＝1.2.

湍流黏性系数仍由公式（2）确定，但是Cμ不再是常数，而是k，ε、平均切应力等变量的函数.

表1中空气密度的变化采用下式计算：

式中：Pop为房间开口处的压力；R为气体常数.

2.2 几何条件与数值方法

由于所模拟分析的房间开口处的流动是双向的，有外界冷气流从开口下部的流入与室内热气流从开口上部的流出，因此门处边界条件较难确定.模拟时定义一个更大的空间作为室外环境（见图3），避免了直接给定门处的边界条件.通过网格独立性分析，总网格数量为108万.

图3 数值模拟几何模型Fig.3 Geometric model for numerical simulation

电加热板表面设为定热流边界条件，墙壁设为绝热边界条件.计算中，采用单精度的分离隐性算法器，压力与速度的耦合运用SIMPLEC算法，压力离散格式为体积力加权算法，3个坐标方向的速度方程和k，ε方程的对流项离散采用2阶上风差分，扩散项的离散采用2阶中心差分.辐射模型采用离散坐标模型.

3 结果分析

分析房间内的温度分布时，取相对温度，即把房间内各点的绝对温度减去环境温度（环境温度为25℃）.

3.1 工况一结果分析

工况一中热源位于房间地面正中（图1中面①，热量输入为1 080W），模拟得到垂直于z轴的房间中心截面的速度、温度和湍流动能的分布分别如图4，5和6所示.

由图4可以看出，3种湍流模型得出的流场分布定性是一致的.室外空气从开口的下部流入房间，沿着地板流向房间深处.遇到地面热源加热之后，气流没有立即上行，而是由惯性的驱动流到开口对面的侧墙后向天花板垂直上升.上行气流未到达天花板高度即往出口方向运动，最终沿开口上半部流出.在房间的中截面下部形成一个顺时针的大涡.房间中部和顶部的气流速度很小，不超过5cm·s－1.最大气流速度出现在门附近，约为40cm·s－1.开口的中和面高度约是开口高度的1／2.由图5可以看出，房间中水平方向温度近乎均匀，垂直方向则产生下低上高的温度分层，在开口顶部附近的温度梯度较大.房间顶部空气温度水平约比底部高30℃.

对比3种模型，Standardk－ε模型得到的开口冷入流速度衰减最快，RNGk－ε模型最慢，而Realizablek－ε模型居中.RNGk－ε模型40℃的等温线进入房间的深度最深.流场分布上，RNG和Realizablek－ε模型得到的结果类似，但两者与Standardk－ε模型的结果有一定差异.图6显示，加热板是引起湍动的原因，流体经过加热板后湍流脉动加强.Standardk－ε模型预测的k值明显较其他2种模型要大，因此其得到的湍流黏性（μt＝cμk2／ε）相应也大，所以其速度衰减快，冷流体入射深度最浅.

文献［4］在对不同湍流模型（从零方程模型直至大涡模拟）对比的研究中发现：在模拟自然对流时，RNGk－ε模型预测温度场的效果最佳，预测速度场的效果一般，预测湍流脉动特征的效果较差；在模拟具有强浮升力的流动时，RNGk－ε模型预测温度场和速度场的效果均较好，其对温度场的模拟精度甚至超过大涡模拟，但预测湍流脉动特征的效果仍然较差.严格来说，Standardk－ε模型是一种高雷诺数模型，对于室内近层流化的低雷诺数区域并不适用.其在模拟立方体的绕流和射流等流态时，会过高地估计湍流动能，得到较大的湍流黏性，最终会导致立方体绕流的剥离区域小，后方逆流弱，主射流变得粗短等与实验不符的现象［10］.

图7—11对比了测量值与3种模型的模拟结果，图中的方框代表实验测量值.

图7 模拟值与实验值对比（温度场）Fig.7 Comparison of simulation value and experimental value（temperature field）

可以看到，温度分布方面，3种模型的差异很小，并且十分接近实验数据.误差比较大的地方出现在靠近房间顶部的位置，模拟得到的温度均偏高，部分原因可能是由于实验中屋顶并不能做到完全绝热所导致.在速度预测方面，大部分点与实验吻合良好，对水平速度的预测优于对垂直速度的预测.部分位置，如A－A′面的测点7，测点6和测点5，误差较大.考虑到室内流动是3维的，且流态复杂，流速很低（相当一部分区域速度低于10cm·s－1），准确模拟存在难度.总体而言，Realizablek－ε模型的结果稍优于其他2种模型；两方程模型对温度场的模拟效果要比速度场更加准确.

3.2 工况二和工况三结果分析

篇幅所限，图12和13分别显示了工况二和工况三Realizablek－ε模型计算得到的z轴中心截面的速度和温度分布.

与工况一类似，工况二和工况三中心截面处空气从开口下部进入房间流至对面墙壁后上行，在未到达屋顶前便返回至开口流出，这是由于房间上部产生很强的温度分层的阻碍作用导致.截面上最大速度出现在地面附近.温度分布方面，随着输入热量的增加，垂直温度梯度变大，开口顶部处的温度梯度也明显增大.

通过模拟值与实验值的对比可以发现［11］，模拟温度值十分接近实验值，且模拟精度高于工况一.在靠近天花板的上部空间，模拟值与实验值仍然吻合较好，这是因为工况二和工况三的输入热量远远高于工况一，屋顶结构的散热效应所起的误差作用相对较小.速度模拟值偏差比较大的地方出现在房间中部，因为该位置受冷空气的水平方向运动和浮升力驱动的竖直方向运动相互作用，流态比较复杂.

3.3 热源引起的通风量对比

开口处的通风量是指房间内外空气交换的速率，表2为各工况的通风量对比.

对工况一而言，Realizablek－ε模型得到的通风量最大，与实验值也最接近.Standardk－ε得到的通风量最小，可能的原因是其为高雷诺数模型，在开口处的低流速 （即低雷诺数）区域，其夸大了湍流黏性，开口处平均流速偏小.

从工况一到工况三，随着加热量的增加，通风量也相应地不断加大，但这一增加不是线性的.工况二的加热量约是工况一的3倍，但通风量是工况一的1.6倍.工况三的加热量约是工况一的8～9倍，但通风量仅是工况一的1.8倍.通风量不仅仅是加热量的函数，也与热源形态和位置密切相关.对比模拟值和实验值，3个工况的模拟值均小于实验值，且加热量越小，相对误差越大.部分原因是实验的测量误差，表3给出了实验时的能量平衡关系，通过门处散出的热量均大于加热量，且工况一的误差最大.实验中门处散热量的确定方法是：在门处布置若干测点，测量各点的速度和温度，再通过面积加权计算得出散热量.理论上而言，需要布置无限多的测量点才能准确得出散热量的值.

表2 模拟与实验通风量对比Tab.2 Ventilation rate comparison of simulation and experiment

表3 实验中3个工况的能量平衡关系Tab.3 Energy balance relations of three operating conditions in the experiment

4 结论

（1）房间内存在垂直方向的温度分层，下部温度低，上部温度高.加热量越大，垂直温度分层的梯度也越大.室外气流从开口下部进入房间，沿地板附近流向对面墙壁，然后上行，大部分气流在到达屋顶之前返回开口流出.

（2）总体而言，Realizablek－ε模型的预测效果稍好于Standardk－ε和RNGk－ε模型.两方程k－ε模型预测温度场的准确度比预测低速自然对流速度场要高.

（3）在房间的中央位置，水平方向的室外流入气流与浮升力气流相互作用导致流态复杂，两方程k－ε模型在模拟这一现象时效果不佳.

［1］CHEN Qingyan.Comparison of differentk－εmodels for indoor airflow computations［J］.Numerical Heat Transfer，Part B：Fundamentals，1995，28：353.

［2］Heskestad G.SFPE handbook of fire protection engineering［M ］.3rd ed.Quincy： National Fire Protection Association，2002.

［3］Novozhilov V.Computational fluid dynamics modeling of compartment fires［J］.Progress in Energy and Combustion Science，2001，27：611.

［4］ZHANG Zhao，ZHANG Wei，ZHAI Zhiqiang，et al.Evaluation of various turbulence models in predicting airflow and turbulence in enclosed environments by CFD：part－ 2：comparison with experimental data from literature［J］.HVAC ＆R Research，2007，13（6）：871.

［5］LI Junmei，LI Yanfeng，CHOW Wanki，et al.Numerical studies on fire－induced thermal plumes［J］.Journal of Thermal Science，2005，14（4）：374.

［6］Nam S，Bill R G.Numerical simulation of thermal plume［J］.Fire Safety Journal，1993，21：231.

［7］Kato S，Murakami S，Ohira N.CFD analysis of thermal plume by means ofk－εmodel［C］／／Summaries of Technical Papers of Annual Meeting.［S.l.］：Architecture Institute of Japan，1996：517－518.

［8］Murakami S，Kato S，Yoshie R.Measurement of turbulence statistics in a model fire room by LDV ［J］.Ashrae Transactions，1995，101：287.

［9］Launder B E，Spalding D B.Lectures in mathematical models of turbulence［M］.London：Academic Press，1972.

［10］Murakami S，Mochida A.3－D numerical simulation of airflow around a cubic model by means of the k －εmodel［J］.Journal of Wind Engineering and Industrial Aerodynamics，1988，31：283.

［11］高乃平，朱彤，董昆，等.火灾基准实验CFD模拟研究报告［R］.北京：西门子研究中心，2010.GAO Naiping，ZHU Tong，DONG Kun，et al.Numerical study on the natural convection in a benchmark fire room［R］.Beijing：Research Report to Siemens Research Center，2010.