基于方法B的有效离散格式

摘要 　本文在方法B网格生成的基础之上，提出了比一般格式更加有效、适应性强的离散格式。作者利用此格式分别对通风房间、带移动顶盖的方型空隙、外掠后台阶的流动进行了数值计算，与 参考文献 提供的值比较，结果基本合理。 　　 关键词 　方法B 离散格式 数值模拟 有效

1 引言 　　 　　在对流体流动、传热及传质问题进行数值模拟时，首先要进行区域离散化。视节点在子区域中位置的不同，可把区域离散化方法分为两类：方法A（外节点法）及方法B（内节点法）。方法A是指先确定节点的坐标再计算相应的界面，使节点位于子区域角顶的方法；方法B是指规定界面位置而后确定节点，使节点位于子区域中心的方法。 　　当网格不均匀时，方法B中节点永远位于控制容积的中心，这样方法B中的节点更能代表控制容积，界面导数的计算精确度也要高。所求解的区域中物性发生阶跃性的变化时，采用方法B时可以较方便地反发生阶跃变化的面作为界面，从而可以避免在同一控制容积内物性发生突变的情形，若采用方法A时实现这一点要困难的多。综述之，一般通风空调工程问题中采用方法B进行区域离散化。本文的出发点也基于此，希望通过常用的方法B离散网格建立有效的、强健的差分格式。 　　 2 对流-扩散方程式的离散及基于方法B的有效差分格式 　　 　　对流换热问题的数学描述包括质量、动量及能量守恒三种方程式，可用统一张量形式可表示为： 　　　　　　　　　　（1） 　　上述方程的四项分别是不稳态项，对流项，扩散项以及源项。据[1]，无论守恒方程式的何种离散或差分格式，必须满足四项基本法则，作者提出的守恒型离散格式亦必须满足这四项原则。现给出一个通用的离散方程式： 　　　　　　　　　　　　　　　　（2） 　　用数值方法求解上述方程（尤其是满足动量守恒的N-S方程）的主要困难来自非线性的一阶导数项及压力梯度项的不妥离散方式。到目前，为解决一些相关问题已提出了许多的差分离散方案：指数、混合、幂函数（乘方格式）、QUICK格式、SHYBRID二阶混合等方案。这些方案，在具体应用中，各有优劣。作者在暖通空调工程数值实验中，基于B方法生成的交错网格推导出一种有效差分格式，特阐述如下： 　　控制方程组（1）中包含有两类变量：矢量与标量。{U，V，W}为矢量，扩散系数、温度、压力、紊流脉动能k、紊流耗散率ε等为标量。这样交错网格的数据存储可表述为：标量数据皆存储在节点位置{P、E、W、N、S等}；而矢量数据皆存储在截面位置{e、w、n、s等}；图1中阴影区为标量控制容积区，标量分布均匀，由节点P等代表。矢量亦均匀分布于各截面，用e等代表。下面皆以二维稳态情形为例进行分析。 　　 　　2．1 标量场 　　首先分析标量场，以φ代表通用标量，由（2）其离散方程的一般形式可写为： 　　　　　　　（3） 　　据幂函数方案AI（I代表e，w，n，s等界面）的取法，上式各系数可表述为： 　　　　　　　　　　　　（4） 　　（4）式中，Peclet数PI=FI/DI，F代表单位截面上对流流量，而D代表单元内物理量的扩散。因为物性参数（标量）亦储存在各个节点上，据前述的四项基本法则，须满足界面上连续性的要求，故可用调和平均来计算界面参数。以为例可由下式求得： 　　　　　　 　　　（5） 　　其他的系数同理可求。 　　 　　2．2 矢量场 　　二维矢量场（即U、V两个速度矢量的分布），而速度的控制容积与标量的不同，其示意图如图2所示。比较图1与图2可以看出：矢量控制容积只是在标量控制容积的基础上沿特定方向发生错移，这也是交错网格的特点之一。以二维稳态离散U方程为例，由N-S方程对上图的控制容积积分可以导出[3]： 　　　　　　　　 　　　　　　　　　　　图1 标量控制容积与离散格式图案 　　　　　　　　 　　　　　　　　　　　　图2 矢量控制容积与离散格式图案 　　 　　[注1]：图2中加"（）"符号的节点一致； 　　[注1]：图1及图2中的 分别等于将要用到的 分别等于将要用到的 。 　　　　 　　同(3)所述，亦采用幂函数格式，上式中的系数可以按（4）表达之，仅需调整I（代表E、W、S、N等节点位置）。从图2可见，Peclet数P=F/D的求解难在界面交点N、S等。界面效点N处流速VN须应用加权平均计算，这样，界面流量FN可表述为： 　　　 　　　　　　　　　　　（6） 　　扩散系数 　　　　　　　　　　　（7） 　　式（7）中系数是标量，而处为界面位置，不存储标量，则可通过周围四个节点上值的面积加权求解（若三维情形，应为体积加权），这样较容易满足上述的四项基本法则。所以： 　　　　 　　　　（8） 　　联立求解，则DN的计算可进一步简化。同理，可以导出DS的表达式： 　　　　 　　　　　　　　　　（9） 　　　 　　　　　　（10） 　　至此，已得到整个U方程系数的求解方法，对于其他矢量方程的系数，同理可求。 　　 　　2．3 新格式与各种差分格式的特性比较及分析 　　对于一维稳态无源项的对流扩散问题，其守恒方程为： 　　　　 　　　　　　　　　　 （11） 　　假定边界条件为： 　　通过计算中间点x=0.5处Φ值，比较了理论解、一阶上风格式、中心差分格式及本文新的格式：中心差分格式随着Pe数的增加而产生发散现象；一阶上风格式和本文提出的差分格式则与理论解有较为接近的值。特别当Pe数增大时，本文的差分格式的更接近理论解。

3 算例应用、比较及分析
　　
　　3．1 单出入口二维通风房间
　　　　　　　　　
　　　　　　　　　　　 　　图3 具有单个出入口的通风房间
　　　　　　　　
　　　　　　　　　　 　　　图4 外掠后台阶的平面图
　　
　　室内流体流态为稳态紊流，采用LRN（Low-Reynolds-Number）模型计算。压力、速度通过SIMPLEC算法耦合求解。入口边界参数给定，出口边界变量的法向梯度为零。以入口速度U为基准速度，以入口尺寸L（1m）为基准尺寸。基于B方法的网格生成为X×Y方向网格数为120×100。之后，作者重新布置网格（80×60），取某垂直截面段（X/L=15）并比较了两者的水平速度分布。如图5可见，两者几乎无差别，这可以说明本文格式应用于工程实际时，对非均匀网格具有较好的适应性。
　　　　　　　　
　　　　　　　　　　　 　图5 两套网格系统截面水平速度的比较
　　
　　对两套网格系统计算过程的收敛史与差分格式未改进的情形进行比较（图6），发现计算精度有较明显的提高，可能的发散趋势被抑制。
　　　　　　
　　　　　　　　图6 迭代过程收敛史的比较（纵坐标即压力方程的余量）
　　　　　　[注]1-未改进的格式；2-格式改进后的第一套网格；3-格式改进后的第二套网格
　　
　　3．2 带移动顶盖的方形空腔流
　　流体流态为稳态层流。以顶盖移动速度U为基准速度，以空腔进深D为基准尺寸。基于B方法的风格生成为X×Y方向网格数为20×20。固体壁面为无滑移条件，因为是层流运动，故动量方程及连续性方程中的流速边界都是如此。
　　　　　　　　
　　　　　　　　　　 图7 截面水平速度的比较（Y/D=0.5）
　　　　　　　　　
　　　　　　　　　　 图8 截面垂直速度的比较（Y/D=0.5）
　　
　　取截面（Y/D=0.5）上的速度矢量与文献[9]比较，由图7及图8可发现两者相差甚微，再次证明本文改进格式的有效性。笔者在计算过程中还比较了一般程度消耗CPU时间一本文所花的时间相差不大，可能是因为物理模型规模较小所致。但是对于规模较大、风格划分不均匀的情形，本文的离散方程具有不可替代的优势。
　　
　　3．3 外掠后台阶的流动
　　如图4所示，H为后台阶高度，并作为基准尺寸：入口及出口尺寸如图4所示；未流速度U取为1m/s，并作为基准速度。来流温度T分布均匀，作为基准温度；壁面为均匀热流条件。基于B方法的网格生成为X×Y方向网格数为60×50。
　　流体流态为稳态二维紊流，服从Bousinesq假设。对此开口流场，出口处采用自由流速边界条件，压力修正方程界为Dirichlet条件。壁面附近采用壁面函数法加以修正。中心轴线上各个变量采用第二类边界条件。
　　图9是截面（X/H=10）上的流速分布，与文献[4]结果相吻合。另外，作者还与文献[5]比较了不同截面上N_u/N_uH数，如图10所示，两者较接近。
　　　　　　　　　
　　　　　　　　图9 截面水平流速分布的比较（X/H=15 ， Re=1.0×10³）
　　　　　　　　
　　　　　　　　　　图10 外掠后台阶局部Nu数分布的比较（Re=1.7×10³）
　　
4 结论及展望
　　
　　本文提出的差分格式从本质上将讲是对幂函数格式的修正及改进，其中应用了中心差分、调和平均等思想，从而达到提高计算精度、减少发散的目的。但从一些运算结果来看，格式并非十分理想，譬如对向后台阶流动（文中未列出），涉及对流换热、壁面边界等问题，结果与文献[2]提供的数据有一定差距；尽管上述运算较原来的算法迭代次数有所减少，但计算时间并未缩短，并且计算机内存容量要求变有提高，这些问题有待今后改进。
　　

参考文献

　　
　　1 Patankar. S.V. Numerical heat transfer and fluid flow. Hemisphere， Newyork， 1980.
　　2 Gooray. A. M. Turbulent Heat Transfer Computations for Rearward-facing Steps and Sudden pipe Expansions. [J]. Heat transfer， 1985，1:70～76.
　　3 Larnder. B. E. The numerical computation of turbulent flows. Comp. Meth. Appl. Mech. Engng 3， 269 (1974)
　　4 Gartling. D. K. A test problem for outflow boundary conditions -flow over a backward-facing step. Int Numer Meth Fluids， 1990， 11: 953--967
　　5 Gooray. A. M and watkins. C. B. Turbulent Heat Transfer Computations for Rearward-facing Steps and Sudden pipe Expansions. Heat transfer， 1985， 1:70～76.
　　6 Chiang， T. P. W. h. Sheu， S. K. Wang. Side wall effects on the structure of laminar flow over a plane-symmetric sudden expansion. Computers && Fluids， 2000， 29: 467～492
　　7 T. JONGEN and Y.P.MARX. Design of an unconditionally stable positive scheme for the k-εand two-layer turbulence models. Computers && Fluids，
　　8 G. D. mallinson and G. DE vahl davis. Three-dimensinal natural convection in a box: a numerical study， J. Fluid Mech. 1977， 83: 1～31
　　9 Latimer. B. R， Pollard. A. Comparison of pressure-velocity coupling solution algorithm. Numer Heat Transfer， 1985， 8: 635～652
　　10 Latimer. B. Rodi and G. Scheuerer， Turbulence models for Near-wall and low-Reynolds number flows: A Review. AIAA J .23， 1308～1319 (1985)
　　11 K. ABE and T. KONDOH， A new turbulence model for predicting fluid flow and heat transfer in separating and reattaching flows---2. Thermal field calculations， Int. J. Heat Mass Transfer 38， 1467～1481 (1995)
　　12 T.P. Chinag， Tony W. h .Sheu. Side wall effects on the structure of laminar flow over a plane-symmetric sudden expansion. Computers && Fluids， 2000， 29: 467～492转贴于