室内空气流动数值模拟的误差预处理法

Error-pretreatment method for numerical simulation of indoor air flow

　　如何快速、准确地模拟和预测工程中需要优化或进行比较的大量工况，一直是CFD（computational fluid dynamics）技术应用于空调通风房间内空气流动的数值模拟仿真存在的问题。作为室内空气流动数值模拟主要组成部分的代数方程求解算法，对计算速度有着很大影响。多重网格法对加快控制室内空气流动的非线性N-S方程的迭代计算收敛速度比较有效^[1]，但是对于室内空气流动数值模拟常用的SIMPLE （semi-implement method of pressure linked equation）算法而言，多重网格法的收敛加速效果并不显著^[1]。基于多重网格法的思想，提出了适于SIMPLE算法的误差预处理方法，以加快迭代计算的收敛速度，更好地适应工程需要。
　　
1 高频和低频误差理论
　　首先简要介绍多重网格法关于高频和低频误差的理论。因为对任意形状计算域中的任意形状的网格，总可以使用保角变换把它们变到矩形计算域上的矩形网格。所以以下讨论均基于矩形计算域和矩形网格。
　　任给定一个划分好网格的计算域，总可以把它划分为A×B×C个区域，使得每个区域上的网格分别在3个方向上是均匀的。设f（x， y， z）为差分后的初始误差，第（u， v， ω）区域的尺寸为2a_u × 2b_v × 2c_ω，网格数为M _u × N _v× K_ω。在每个区域内，3个方向上分别只有周期等于或低于2M _u ，2N _v，2K_ω的分量能够被表示出来。因此：
　　
　　　[A_{i，j，k，u，，v，w} sin +B _{I，j，k，v，w} cos ]·[C_{i，j，k，u，，v，w} sin +D_{i j，k，v，w} cos ]·[E_{i，j，k，u，，v，w} sin +F _{i，j，k，v，w} cos ]， 　　
　　 　A_{i，j，k，u，，v，w} ，B_{i，j，k，u，，v，w} ，C_{i，j，k，u，，v，w} ，D_{i，j，k，u，，v，w} ，E_{i，j，k，u，，v，w} ，F_{i，j，k，u，，v，w} ∈(u ， v ， w)
　　　A_{i，j，k，u，，v，w} ，B_{i，j，k，u，，v，w} ，C_{i，j，k，u，，v，w} ，D_{i，j，k，u，，v，w} ，E_{i，j，k，u，，v，w} ，F_{i，j，k，u，，v，w} = 常数
　　 A_{i，j，k，u，，v，w} ，B_{i，j，k，u，，v，w} ，C_{i，j，k，u，，v，w} ，D_{i，j，k，u，，v，w} ，E_{i，j，k，u，，v，w} ，F_{i，j，k，u，，v，w} (u ， v ， w)
　　 A_{i，j，k，u，，v，w} ，B_{i，j，k，u，，v，w} ，C_{i，j，k，u，，v，w} ，D_{i，j，k，u，，v，w} ，E_{i，j，k，u，，v，w} ，F_{i，j，k，u，，v，w} = 0 　（1）

　　其中：t _{x，u，v，w} =0，1，…，M _u ；t _{y，u，v，w} =0，1，…，N v ；t _{z，u，v，w} =0，1，…，K_w ； 　　
　　根据多重网格理论，I / M _u ，j / N _v，k/K _ω 小的项（称为"低频误差"），收敛速度很慢，是妨碍迭代法快速收敛的主要因素。而如果减小M _u ， N _v， K _ω，即把网格粗化，就可以使误码差频率变高，并使一部分细网格上的低频误差转达化为粗网格上的高频误差。如果交替运用粗细网格进行迭代，利用粗网格快速消除低频误差，在细网格上提高精度，就可以加快收敛速度。
　　
2 误差预处理的基本思想和算法
　　2．1 误差预处理法
　　传统迭代法是在整个计算区域的所有网格上同时进行迭代计算至收敛。但如果能先求出一些点的值，然后再根据它们来预测所有点的值，只要预测方法得当，从误差削减的角度而言，就可以获得比迭代法高得多的效率。可先在粗网格上迭代计算至问题收敛，然后借助该预测方法把粗网格上的值映射到细网格上作为迭代初值，再进行迭代计算直至收敛。该方法的基本思想在于通过在粗网格上用比较少的时间（相对于在细网格上迭代）来获得问题的大致描述，以便通过线性插值大量削减误差，称其为"误差预处理法"。假设已知如图1计算域中点A，B，C，D，E，F，G，H的真实值。各点坐标分别为：
　　A（x - d a ， y + d y ， z + d z）
　　B（x + d x ， y + d y ， z + d z）
　　C（x - d a ， y + d y ， z - d c）
　　D（x + d x ， y + d y ， z - d c）
　　E（x - d a ， y - d b ， z + d z）
　　F（x + d x ， y -d b ， z + d z）
　　G（x - d a ， y - d b ， z - d c）
　　H（x + d x ， y - d b ， z - d c）
　　　　　　　　　
　　　　　　　　　　　　　　　图1 计算示意图

　　由流体流动的输运控制方程可知，点I上的值只受其相邻点上的值的影响，即只受A，B，C，D，E，F，G，H点的影响。因此，可以用这些点上的值来预测一个误差比较小的I的估计值。
　　很自然，用这8个点的加权平均值作为I的值：

　　P （I）= M_A f (A) + M _B f (B) + M_C f (C) + M_D f (D) + M_E f (E) + M_F f (F) + M_G f (G) + M_H f (H)，　　　 (2)

　　其中，P（I）表示预测值，f( * ) 表示相应点的真实值，M x表示加权系数：
　　由守恒律，有：

　　M_A + M _B + M_C + M_D + M_E + M_F + M_G + M_H = 1 　　　　　　　　　　(3)

　　离I比较近的点对I的影响应该比较大，即较近的点加权系数应较大。
　　结合式（3），可以根据面积律导出加权系数^[1]：

　　M_A =(d x × d b × d c )/[(d x + d a )(d y + d b )(d z + d c )]
　　M_B =(d a × d b × d c )/[(d x + d a )(d y + d b )(d z + d c )]
　　M_C =(d x × d b × d z )/[(d x + d a )(d y + d b )(d z + d c )]
　　M_D =(d a × d b × d z )/[(d x + d a )(d y + d b )(d z + d c )]
　　M_E =(d x × d y × d c )/[(d x + d a )(d y + d b )(d z + d c )]　　　　(4)
　　M_F =(d a × d y × d c )/[(d x + d a )(d y + d b )(d z + d c )]
　　M_G =(d x × d y × d z )/[(d x + d a )(d y + d b )(d z + d c )]
　　M_H =(d a × d y × d z )/[(d x + d a )(d y + d b )(d z + d c )]

　　将式（4）代入式（2），得

　　P（ x ， y ， z ）= [(d a × d b × d c )f (x + d x，y +d y，z + d z)+(d a × d b × d z )f (x + d x，y +d y，z - d c)+(d a × d y × d c )f (x + d x，y -d b，z + d z)+(d a × d y × d z )f (x + d x，y -d b，z - d c)+(d x × d b × d c )f (x - d a，y +d y，z + d z)+(d x × d b × d z )f (x - d a，y +d y，z - d c)+(d x × d y × d c )f (x - d a，y -d b，z + d z)+(d x × d y × d z )f (x - d a，y -d b， z - d c) ]/[(d x + d a )(d y + d b )(d z + d c )] 　　　　　　　　　　　　　　　　　　　　　　(5)

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

　　类似写出其它项，把它们代入式（5），并考虑到

　　O (d x² ) = O (d a² )，O (d y² ) = O (d b² )，O (d z² ) = O (d c² )　　　　　(7)

　　可以得出： P（x， y， z）-f（x， y， z）= O (d x² + d y² + d z²) 　　　　　　(8)

　　上式表明上述插值方法具有二阶精度，而且它是显式方法，效率远远高于迭代方法。综上所述，误差预处理的算法如图2所示。
　　　　　　　　　　
　　　　　　　　　　　　　图2 误差预处理法的计算流程

　　 图2中当运算进行到限制和插值时即采用上述的预测方法。由此可见，误差预处理法本身不仅是数学方法，而且也蕴涵了一定的物理意义。 　　
　　2．2 误差预处理法计算速度分析
　　如前所述，对一个数值模拟过程，先把细网格合并成粗网格进行计算，得到问题的一个粗略描述，通过它们预测出细网格上的估计值，以其为初值再进行迭代。由于在粗网格上的每一步计算量大大低于在细网格上的计算量，而且预测过程的误差消除效率高于迭代法，故只要将由粗网格迭代和网格值映射所带来的迭代步数的增加控制在一定范围以内，就可以加快迭代收敛速度。
　　对一维问题，由细网格生成的粗网格的数目大约是细网格数止的1/2，对于三维问题，大约为1/8。SIMPLE算法每次迭代的计算量大约正比于网格数，故对实际三维工程问题，细网格上每项次迭代的计算量是粗网格上的8倍。设普通迭代法需要M步收敛，耗时T（M）；误差预处理法在粗网格上迭代Na步，在细网格上迭代N步后收敛，各耗时Ta (Na)，T（N）。

　　如果误差预处理法有效，即　T（M）&> Ta (Na)+T（N）， 　　　　　　　　　(9)

　　而　Ta (Na)+T（N）≈T (Na/8) +T（N）= T (Na/8+N） 　　　　　　　　　(10)

　　应满足　Na/8+N&< M 　　　　　　　　　　　　　　　　　（11）

　　在实践中我们发现，一般的工程问题都是符合这一条件的。另外，由于每次粗网格上的迭代计算量仅约为细网格上的1/8，因此粗网格上的迭代时间约为细网格迭代时间的1/8。所以，采用三层或三层以上的网格所能节省的时间是非常有限的，最多不超过由细网格直接生成的第一层粗网格上的迭代时间，即细网格迭代时间的1/8。故我们认为，在误差预处理法中，采用粗细两重网格已经足够了。
　　误差预处理法比多重网格法更简单，且易实现，它适用于SIMPLE算法，从而避免了多重网格法对基于SIMPLE算法的迭代计算收敛不快的缺点。实际上，任何基于差分议程格式的迭代法都可以运用误差预处理法。根据误差预处理法的算法过程，只要原方法有效准确，误差预处理法就有效准确。
　　

计算方法

网格划分

计算至收敛时间/min

例1

例2

例1

例2

误差预处理法传统迭代法

32×32×32

32×32×32

180

314

36×36×32

36×36×32

312

501