CFD数值模拟的系统误差反馈及其实现

赵福云 汤广发 刘娣 张泠

简介： 全文回顾了近年来CFD领域数值计算误差及不确定度的研究进展，剖析了系统误差反馈的基本程序及CFD在通风空调领域的应用过程中所遇到的特殊困难，并提出了CFD数值模拟的系统误差反馈思路。

关键字：计算流体动力学 系统误差 通风空调 　在实验流体力学及传热学方面，实际例子很多，比较常见且广为引用的有如下几种[25][26][27][28][29]：

1 长方形截面90º弯头中压力系数的分布

2 圆管突然收缩区的轴向流速分布

3 Benard流动

4 竖直狭缝中二维自然对流

5 方形空腔内弧在平板的自然对流

6 三维顶盖驱动分层流

显然这类实验也是广泛针对工业管道流动及换热的特点进行设计或布置的，具有明显的局限性，即不能直接为建筑环境领域的CFD研究者所采用。

近些年来，在建筑环境领域也悄然兴起了室内气流组织实验、自然通风建筑原型比例实验、通风测试、室外环境测试等浪潮，美国、日本等国家的建筑环境的CFD专家更是乐衷于此[3][30][31]。但是这些测试或实验普遍不属规范性实验，并且实验数据难以共享，更为严峻的是，对流体流动及传热计算所需的边界条件缺乏详细说明或规定。总之，建筑环境领域可靠的通风气流组织标准化测试实验还有很长一段路要走。

3 通风空调领域CFD系统误差反馈的研究

文献[3]已提及在建筑环境领域，美国学者Q.Chen等已经开始了建筑环境领域的CFD的验证与考核研究，其基本思路局限于二维或三维的自然对流、强迫对流或混合对流等具体实例的计算比较，实验手段及程序缺乏规范或程序，还没有上升到理性研究高度。另外，前文业已陆续讲到CFD在通风空调领域应用中的系统误差研究所存在的不足，现从系统误差反馈的全局出发，继续深入探讨该问题。

3.1 流体流动及换热的新特点及困难

在通风空调领域，我们普遍遇到的是建筑室内外空气的流动状况及热分布规律问题，其中主要的难点在于：

1 高大空间(热边界、速度矢量边界等条件难以确定)，多风口多障碍物(包括出入风口参数的确定)，几何外形复杂(如圆形屋顶)；

2 自然对流与强迫对流并存(如置换通风)，温差及浓差驱动流并存(双扩散效应，如烟气流动、壁面有机物材料扩散(VOC)等)；

3 热辐射(室内辐射采暖、灯光照明，室外太阳辐射等)及墙体耦合换热；

4 室外及室内气流的相互干扰，活动人员及出入建筑空调区域的设备等对气流组织的干扰。

1933年，英国人Thom首次用数值方法求解了二维粘性流体偏微分方程，CFD由此诞生，但是针对流体流动及换热模型的讨论就从来没有休止过。同样地，对于建筑环境流体应用研究为主的领域，模型问题更是如此，并且在上述新的特点下又有了新的困难。综合目前的相关文献，有如下内容：

a 湍流模型不完备。近年来，美国学者Q.Chen等提出了新的零方程模型，主要用于计算室内气流计算，其可信度及适用范围值得考核；日本学者Murakami[31]等考虑了室外气流模拟时标准k-ε模型的不足及其改进方案，但局限性依然存在；当室内气流受温差影响较大时，必须对湍流模型进行修正，Markatos等人早在20世纪八十年代初就已开始研究，但结果及适用范围仍不尽人意。

b 壁面边界的处理没有规范。采用湍流模型的最大问题之一是壁面附近层流粘性支层及过渡区的处理，目前通行的做法是在壁面附近添加大量的网格，采用低雷数湍流模型，使模型由固壁表面一直用到湍流核心，但耗费的计算时间是显而易见的；若采用不断尝试各种壁面网格间距的办法，便使得计算过程及实践变得复杂，且引入了人的主观因素；另外，过去应用较多的壁面函数法，大量的经验常数已经使人们怀疑它的生命力。

c 辐射模型及计算方法亟待改进。从散见的国际传热及流体杂志，可以发现近些年来，科研及工程技术人员对辐射现象及其计算愈加重视，但是考虑复杂几何形体的空间辐射计算，却很少见，尤其是考虑室外太阳辐射得热的流体流动及传热计算。国际领域盛行的Montle-Carlo法或DOM法不仅限于解决小型空间的辐射，并且计算的收敛性能也待改进。

3.2 数值模拟的基本思想

鉴于HVAC领域通风气流组织的实际情况，一般采用RANS方法(以标准k-ε双方程模型为代表)建模、FVM(有限容积法)离散(对流项采用QUICK 等格式，扩散项采用具有二阶精度的中心差分格式)、TDMA+SOR法求解离散代数方程、SIMPLE类压力速度预估校正思想耦合等组合过程对所选实际的通风气流组织进行数值计算。故主要的分析和研究对象是：湍流模型误差、对流项离散格式误差(雷诺数偏高时)、代数方程迭代求解误差、SIMPLE类方法预估校正的收敛性、边界条件及初始条件的实施、网格分辨率(网格雷诺数的影响)、时间截差(非稳态问题)等几个方面，除湍流模型误差、初边界条件误差、物性参数误差等属考核阶段外，其他皆归于验证阶段的分析和探讨。

4 CFD数值模拟系统误差反馈思路 4.1 系统误差反馈之一：验证

从普通雷诺时均控制方程组可以发现，湍流流动的动量输运方程与层流流动的动量输运方程除了湍流扩散粘性系数与源项有差别外，其他并无差异。在假设湍流扩散粘性系数已知的情况下，完全可以采用层流的验证方法进行，即可以获得相同的误差动量输运方程及质量方程。即：

ErrorMomentum, Turbulent=ErrorMomentum, Laminar+ErrorSource

ErrorMomentum, Laminar=ErrorConvecion term+ErrorDiffusion term+ErrorGrid+ErrorSIMPLE Iteration

采用前述的单网格后向误差估计方法或更加高效可行的方法，对采用FVM离散了的层流控制方程，得出其误差控制方程组(包括对流项离散格式的截断误差分析)，再根据网格分布、动量边界等已知的信息获得ErrorMomentum, Laminar。显然，收敛的计算结果所获得的湍流扩散粘性系数是没有任何问题；但若是在耦合湍流动能及湍流动能耗散率方程计算过程中，并未获得收敛解，则谈论其反馈的意义也就消失了。这样，问题就简化成了湍流粘性系数对湍流流动的动量输运误差方程的影响了。一种思路：将此误差因素计入下面的考核之中，及湍流模型误差之中。一种思路：在作更多的迭代计算，考虑迭代误差的影响因素，但同时需保证获得收敛解或迭代不完全的解。

4.2 系统误差反馈之二：考核

复杂湍流流动的验证是很困难的，从3.2节可以发现，因为湍流粘性系数的原因，仍需湍流模型进行较核，即湍流模型的考核是必然的；另外，模型本身的不完善，也促使了这项内容的诞生。据流体流动的实际特征，以下三步是比需的考核项目：

Step1 经典数值计算或实验模型：顶盖驱动流、封闭空腔、台阶流等

考核源程序和计算思想的正确性。

Step2 二维或三维混合对流数值实验，ReCase=ReReal，RaCase=RaReal，PrCase=PrReal

考核该数值模型的计算稳定性，以及对于高大空间或复杂流动的适应性(建议采用三套单网格系统或多重网格系统，逐步加密，以考核网格分布特性及固壁边界的流态转化)；保证网格雷诺数的相近，以实现小空间流场计算往大空间流场计算的可持续性。

Step3 标准建筑环境气流室(建议参照文献[24]提出的实验数据记录格式提供)

考核该数值模型在建筑环境流体计算领域的可行性，由此推广至建筑环境的其他应用。

4.3 CFD计算反馈文件(Documentation)

将模型建立、计算对象、初边界条件、验证过程及结果、考核内容及分析结论整理成文(暂未考虑辐射)：

F1 RANS方法建模，采用的湍流模型(零方程、一方程、两方程等)；

F2 计算对象为高大空间、室内置换通风、室外流场分布或烟气流动等；

F3设置边界条件，实验测得参数应注明仪器的精度，尤其对湍流模型的壁面边界条件、热流或物料流边界需要加以仔细说明；

F4验证过程采用的验证方法、误差方程、湍流粘性系数的计算等需要详尽说明，验证得到的误差及不确定度分析需要详细说明；

F5考核过程所采用的经典数值计算或实验模型及计算比较结果、三套网格系统下二维或三维混合对流数值实验及比较结果、标准建筑环境气流室的模拟及比较结果。

需要说明的是，若F4所得结果不能令人满意(置信度取95％，不确定度为5%)，进行第F5步是毫无意义的，也是徒劳的。由此可见，对离散方程的离散误差、尤其是对流项的截断误差分析显得极为重要，它关系到CFD数值模拟及其反馈计算的全局。

需要补充的是，对于辐射的摄入，可以有两种办法，一是采用独立网格系统对辐射单独计算，然后耦合计算结果，需单独列出辐射能的误差分析方程；二是对同一套网格系统，将辐射能耦合入能量方程，不需单独列出辐射能的误差分析方程，但边界上要注意区分一般热流与净辐射得热。

5 参考文献

[1] 陶文铨. 数值传热学的近代进展. 科学出版社

[2] R.W.Douglass, J.D.Ramshaw. Perspective: Future Research Directions in Computational Fluid Dynamics. Journal of Fluid Engineering, Vol116, P212-215, 1994

[3] Qingyan Chen, Zhiqiang Zhai. How Realistic is CFD As a Tool For Indoor Environment Design and Studies Without Experiment. The 4th International Symposium on HVAC, Beijing, China, October 9-11, 2003, P62-77

[4] Fred Stern, et al. Verification and Validation of CFD Simulations. IIHR Report No.407, Iowa Institute of Hydraulic Research, University of Iowa, Sep.1999.

[5] Guide for the Verification and Validation of Computational Fluid Dynamics Simulations, AIAA-G077-1998.

[6] Best Practice Guidelines, Version 1.0, ERCOFTAC Special Interest Group on “Quality and Trust in Industrial CFD”, January 2000.

[7] ITTC Quality Manual

[8] Roache P.J. Verification and Validation in Computational Science and Engineering. Hermosa Publishers, 1998.

[9] Eca L., Hoekstra M. On the Numerical Verification of ship Stern Flow Calculations. 1st Marnet CFD Workshop, Barcelona 1999.

[10] Eca L., Hoekstra M. On the Numerical Verification Procedures in Computational Fluid Dynamics. 2st Marnet CFD Workshop, Copenhagen 1999.

[11] Jasak H. Error Analysis and Estimation in the Finite Volume Method with Applications to Fluid Flow. PhD Thesis, Imperial College, University of London, 1996.

[12] OdenJ.T. Prudhomme S. New Approaches to Error Estimation and Adaptivity for Navier-Stokes Equation. Tenth International Conference on Finite Elements in Fluids, Tucson, January 1998.

[13] Zurigat Y H, Ghajar A J. Comparative study of weighted upwind and second order upwind difference scheme. Numer Heat Transfer Part B, 1990. 18: 61-80.

[14] Braga W. On the use of some weighted upwind schemes for strongly convection flows. Numer Heat Transfer, Part B, 1990. 18: 43-60.

[15] Ethier C R, Steinman D A. Exactly fully 3D Navier-Stokes solution for benchmarking. Int J Numer Methods Fluids, 1994. 19(5): 369-376.

[16] De Vahl Davis G, Jones I P. Natural convection in a square cavity: a comparison exercise. In: Lewis R W, Morgan K, Schrefler B A, eds. Numerical methods in thermal problems. Pineridge Press, 1981, 552-572.

[17] Ghia U, et al. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method. J comput Phys, 1982. 48: 387-411.

[18] Barakos G. Misoulis E. Natural convection flow in a square cavity revisited: laminar and turbulent models with wall functions. Int J Numer Methods Fluids, 1994. 18(7): 695-719.

[19] Ooserlee C W, Wesseling P, Segal A, Brakkee E. Benchmark solutions for the incompressible Navier-Stokes equations in general co-ordinates on staggered grids. Int J Numer Methods Fluids, 1993. 17: 301-321.

[20] Demirdzic I, Lilek Z, Peric M. Fluid flow and heat transfer test problems for non-orthogonal grids: benchmark solutions. Int J Numer Methods Fluids, 1992. 15: 339-354.

[21] Gartling D K. A test problem of outflow boundary conditions—flow over a backward facing step. Int J Numer Methods Fluids, 1990. 11: 953-967.

[22] Leone Jr J M. Open boundary condition symposium benchmark solution: stratified flow over a backward-facing step. Int J Numer Methods Fluids, 1990. 11: 969-984.

[23] Tang L Q, et al. Transient solutions for three-dimensional lid-driven cavity flows by a least squares finite element method. Int J Numer Methods Eng, 1995. 21: 413-432.

[24] Faeth G M, Samuelsen G S. Fast reduction nopremixed combustion. Progr Energy Combustion, 1986. 12(4): 305-370.

[25] Kim W J, Patel V C. An experimental study of boundary-layer flow in a curved rectangular duct. In: FED Vol. 146, Data for validation of CFD codes. ASME, 1993. 13-28.

[26] Buckle U, Durst F. Investigation of laminar flow in a pipe with sudden contraction of cross sectional area. In: FED Vol. 146, Data for validation of CFD codes. ASME, 1993. 61-78.

[27] Fahadieh R, Tankin R S. Interferometric study of two-dimensional Benard convection cells. J Fluids Mech. 1974, 66:739-752.

[28] Veit C M, Arpaci V. Stability of natural convection in a vertical slot. J Fluids Mech, 1969. 36: 1-5.

[29] Mohamad A A, Viskanta R. Flow and heat transfer in a lid-driven cavity filled with a stably stratified fluid. Appl Math Modelling, 1995. 19: 465-472.

[30] J.W.Weathers, J.D.Spitler. A comparative study of room airflow: numerical prediction using computational fluid dynamics and full-scale experimental measurements. ASHRAE Transactions, 1993, 99: 144-157.

[31] Murakami Shuzo, et al. Numerical prediction of flow around a building with various turbulence models: comparison of EVM, ASM, DSM, and LES with wind tunnel tests. ASHRAE Transactions, 1996, 10: 741-753.

[32] Demuren A O, Wilson R V. Estimating uncertainty in computations of two-dimensional separated flows. ASME J Heat Transfer, 1994, 116: 216-220.

[33] Duraiswami R, Prosperetti A. Orthogonal mapping in two dimensions. J Comput Phys, 1992. 98: 254-268.

[34] Huang Huaxiong, Prosperetti A. Effect of grid orthorgonality on the solution accuracy of the two-dimensional convection-diffusion equation. Numer Heat Transfer, Part B, 1994. 26: 1-20.

[35] Roache P J. Prospective: A method for uniform reporting of grid refinement studies. ASME J Fluids Eng, 1994. 116: 405-413.

[36] Ferziger J H. Estimation and reduction of numerical error. FED, 1993. 158: 1-7.

[37] Brian Van S. A posteriori error estimation for finite volume simulations of fluid flow. Master Thesis (University of Waterloo), July 8,1996