, 张敬奎21. 东北大学 材料电磁过程研究教育部重点实验室，辽宁 沈阳 110819;
2. 武汉科技大学 省部共建耐火材料与冶金国家重点实验室，湖北 武汉 430081;
3. 大连理工大学 能源与动力工程学院，辽宁 大连 116024
收稿日期：2015-11-16
基金项目：国家自然科学基金资助项目 (11402180)。
作者简介：陈元元 (1980-), 女, 河南周口人, 东北大学博士研究生;
李本文 (1965-), 男, 湖南澧县人, 大连理工大学教授, 博士生导师。
摘要：采用Chebyshev配置点谱方法对局部热平衡状态下多孔介质方腔内的自然流动进行了模拟，使用Chebyshev-Gauss-Lobatto配置点对无量纲化的控制方程进行了空间上的离散，离散方程组采用高效矩阵对角化方法进行了求解.将所得结果与已有文献进行了对比，计算结果吻合良好.为验证该数值方法的精度，构造了一个精确解对该方法的求解误差进行了测试，结果表明，Chebyshev配置点谱方法具有很高的计算精度.最后，在验证程序正确性的基础上，研究了Ra对流场、温度场及努塞尔数的影响.
关键词：Chebyshev配置点谱方法多孔介质自然流动局部热平衡数值模拟
Simulation of Natural Convection in Square Porous Cavity Based on Chebyshev Spectral Method
CHEN Yuan-yuan1,2, LI Ben-wen3
, ZHANG Jing-kui21. Key Laboratory of Electromagnetic Processing of Materials, Ministry of Education, Northeastern University, Shenyang 110819, China;
2. The State Key Laboratory of Refractories and Metallurgy, Wuhan University of Science and Technology, Wuhan 430081, China;
3. School of Energy & Power Engineering, Dalian University of Technology, Dalian 116024, China
Corresponding author: LI Ben-wen, professor, E-mail: heatli@dlut.edu.cn
Abstract: The Chebyshev spectral collocation method (CSCM) was used to simulate natural convection in the square porous cavity under a local thermal equilibrium condition. Using the Chebyshev-Gauss-Lobatto collocation points, the dimensionless govering equations can be spatially discredited. The matrix diagonalization method was used to solve these discredited equations. The simulation results were in well agreement with the existing data. Meanwhile, the exact solutions were further constructed to valid the method and the CSCM effectiveness was proved. Finally, the effects of Rayleigh number on the streamlines, the isotherms and the Nusselt number were studied based on the CSCM simulations.
Key Words: Chebyshev spectral collocation methodporous medianatural convectionlocal thermal equilibriumnumerical simulation
多孔介质内的自然流动现象广泛存在于自然界和工程技术领域，如：建筑物隔热、地下污染物迁移、谷物干燥、太阳能收集、填充床内催化反应放射性废物容器冷却等.二维多孔介质方腔内的自然对流模型是研究多孔介质热对流的基本模型.其特点是方腔左右壁面分别维持恒定的一高一低的温度，上下壁面为绝热壁面.由于此模型边界条件简单，因此常用来验证不同数值解法的优劣性.
近几十年来，不同的学者运用了不同的方法对多孔介质内自然对流问题进行了研究.文献[1-2]使用交替隐式有限差分方法研究了多孔介质倾斜方腔和梯形腔内自然对流，并且应用有限容积法研究了局部非热平衡下多孔介质方腔内的自然流动[3]；文献[4]使用Galerkin有限元法研究了多孔介质波纹曲面腔内的自然流动；文献[5]应用稳定的有限元二阶投影法研究了随机构造的多孔介质内的自然流动；文献[6]开发了一种模拟多孔介质内自然流动的热格子-波尔兹曼模型等.不同的研究者采用不同的数值方法对物理问题进行研究，而据作者了解，目前并没有人采用谱方法对多孔介质内自然对流进行模拟.
谱方法的主要优点是精度高，是与有限差分、有限元并驾齐驱的求解偏微分方程的三大数值方法之一.其对光滑函数具有指数性逼近的谱精度, 可以以较少的网格点得到较高的精度, 具有谱解析性和全域性.因此，谱方法在流体流动[7]、辐射换热[8]、磁流体[9]、量子力学[10]等方面有了广泛的应用.本文主要应用Chebyshev配置点谱方法来模拟局部热平衡状态下多孔介质方腔内的自然流动，验证谱方法在该问题上的适用性及其精度，并在此基础上分析无量纲参数Ra对流动和温度分布的影响.
1 物理问题与数学模型本文所研究的物理模型与坐标系统如图 1所示，该模型中，二维方腔的左侧垂直壁面为维持恒定温度Th的热壁面，右侧垂直壁面为维持恒定温度Tc的冷壁面，且Th>Tc，上下壁面为绝热壁面.同时假设：流动为稳定的、不可压缩的层流流动，用Darcy模型描述多孔介质内动量方程；多孔介质内固体与流体处于局部热平衡状态；多孔介质结构均匀，具有各向同性渗透率；除密度以外，多孔介质各物性参数认为是常数.
图 1(Fig. 1)
 | 图 1 物理模型与坐标系统Fig.1 Physical model and coordinate system |
根据以上假设，描述图 1所示物理问题的无量纲控制方程、初始条件和边界条件如下：
 | (1) |
 | (2) |
 | (3) |
 |
2 Chebyshev配置点谱方法为了使数值求解过程更为稳定，在应用Chebyshev配置点谱方法求解控制方程及边界条件 (1)~(3) 时，在方程 (1) 上增加非稳态项 (?Ψ/?τ)[11].采用Chebyshev-Gauss-Lobatto (CGL) 配置点对无量纲化的控制方程进行空间上的离散.
2.1 数值求解过程对控制方程中非稳态项采用一阶有限差分格式进行离散，时间步长为Δτ.方程中的Laplacian项在时间上采用隐式 (n+1时层)，其他项采用显式 (n时层) 进行处理.在Chebyshev配置点谱方法的理论中，无量纲流函数和温度可用配置点上的值由Lagrange插值多项式进行近似[12]，将插值式代入控制方程，并引入初始条件与边界条件可得
 | (4) |
 | (5) |
 | (6a) |
 | (6b) |
 | (6c) |
 | (6d) |
 | (6e) |
 | (6f) |
由于A, B, P和Q矩阵的特征值和特征向量均为实数，可采用高效的矩阵对角化方法[13]对方程组进行求解.求解过程中，时间步长Δτ=10-5，计算终止判据为前后时层待求参数之差小于10-8.
2.2 网格独立解测试对于二维多孔介质方腔，竖直壁面的平均努塞尔数为
 | (7) |
表 1(Table 1)
 表 1 不同网格节点数下平均努塞尔数Table 1 Average Nusselt number for different grid numbers
| 表 1 不同网格节点数下平均努塞尔数 Table 1 Average Nusselt number for different grid numbers | ||||||||||||||||||||||||||||||
2.3 有效性验证表 2对当前数值算法所求出的竖直壁面努塞尔数与其他研究者所求出的结果进行了比较，从表中看出，当前数值解法与其他解法所得结果误差均小于3%.
表 2(Table 2)
表 2 平均努塞尔数的比较Table 2 Comparison of average Nusselt number
| 表 2 平均努塞尔数的比较 Table 2 Comparison of average Nusselt number | |||||||||||||||||||||||
2.4 误差测试为了测试本数值方法的精确性，构造了一个精确解对其进行了误差测试.
精确解的形式如式 (8) 所示，精确解对应的控制方程和边界条件详见作者前期工作[16].
 | (8) |
 |
图 2和图 3分别为当X等于0.5时，不同的Y对应节点及Y等于0.5，不同的X对应节点上数值计算解与精确解之间的误差.从图中可以看出，其流函数最大误差为1.76×10-12，而温度最大误差为2.79×10-11.由此表明，采用谱方法计算二维多孔介质方腔的自然流动问题具有很高的计算精度.
图 2(Fig. 2)
 | 图 2 无量纲流函数数值解与精确解之间的E2误差Fig.2 E2 errors of stream function between exact solutions and CSCM results |
图 3(Fig. 3)
 | 图 3 无量纲温度数值解与精确解之间的E2误差Fig.3 E2 errors of dimensionless temperature between exact solutions and CSCM results |
3 无量纲参数Ra的影响不同Ra下，流场和温度场的分布如图 4所示，从图中可以看出，随着Ra的增大，流场内中心涡的形状发生变化，且流动逐渐加强；温度场内等温线发生偏移，等温线上部向冷壁面集聚，等温线下部则向热壁面集聚，致使等温线聚集的部位传热加强.
图 4(Fig. 4)
 | 图 4 流线和等温线分布Fig.4 Distribution of stream lines and isotherms (a)—Ra=10，流线；(b)—Ra=10，等温线; (c)—Ra=100，流线；(d)—Ra=100，等温线; (e)—Ra=1 000，流线；(f)—Ra=1 000，等温线. |
不同Ra下，竖直壁面平均努塞尔数的变化趋势如图 5所示.从图中可以看出，随着Ra的增大，平均努塞尔数增大，且二者在测试范围内 (10≤Ra≤1 500) 呈现线性变化的规律.
图 5(Fig. 5)
 | 图 5 Ra对竖直壁面平均努塞尔数的影响Fig.5 Effect of Ra on average Nusselt number at vertical walls |
4 结论1) 本文构造了一个精确解对误差进行了测试，结果表明本方法具有很高的计算精度.
2) 在验证数值方法精度的基础上，分析了无量纲参数Ra对多孔介质内流场、温度场及竖直壁面平均努塞尔数的影响.发现，随着Ra的增加，腔体内流动得到加强；垂直壁面上的平均努塞尔数随着Ra的增大而增大，在10≤Ra≤1 500区间内，二者呈现线性变化的规律.
参考文献
| [1] | Baytas A C, Pop I. Free convection in oblique enclosures filled with a porous medium[J].International Journal of Heat and Mass Transfer, 1999, 42(6): 1047–1057.DOI:10.1016/S0017-9310(98)00208-7 |
| [2] | Baytas A C, Pop I. Natural convection in a trapezoidal enclosure filled with a porous medium[J].International Journal of Engineering Science, 2001, 39(1): 125–134. |
| [3] | Baytas A C, Pop I. Free convection in a square porous cavity using a thermal nonequilibrium model[J].International Journal of Thermal Sciences, 2002, 41(6): 861–870. |
| [4] | Misirlioglu A, Baytas A C, Pop I. Free convection in a wavy cavity filled with a porous medium[J].International Journal of Heat and Mass Transfer, 2005, 48(9): 1840–1850.DOI:10.1016/j.ijheatmasstransfer.2004.12.005 |
| [5] | Ma X, Zabaras N. A stabilized stochastic finite element second-order projection method for modeling natural convection in random porous media[J].Journal of Computational Physics, 2008, 227(8): 8448–8471. |
| [6] | Gao D Y, Chen Z Q, Chen L H. A thermal lattice Boltzmann model for natural convection in porous media under local thermal non-equilibrium conditions[J].International Journal of Heat and Mass Transfer, 2014, 70: 979–989.DOI:10.1016/j.ijheatmasstransfer.2013.11.050 |
| [7] | Subich C J, Lamb K G, Stastna M. Simulation of the Naviere-Stokes equations in three dimensions with a spectral collocation method[J].International Journal for Numerical Methods in Fluids, 2013, 73(1): 103–129. |
| [8] | Li B W, Sun Y S, Yu Y. Iterative and direct Chebyshev collocation spectral methods for one-dimensional radiative heat transfer[J].International Journal of Heat and Mass Transfer, 2008, 51(5): 5887–5894. |
| [9] | Zhang J K, Li B W, Chen Y Y. Hall effects on natural convection of participating MHD with thermal radiation in a cavity[J].International Journal of Heat and Mass Transfer, 2013, 66: 838–843.DOI:10.1016/j.ijheatmasstransfer.2013.07.090 |
| [10] | Li Z C, Chen S Y, Chen C S, et al. A spectral collocation method for a rotating Bose-Einstein condensation in optical lattices[J].Computer Physics Communications, 2011, 182(6): 1215–1234.DOI:10.1016/j.cpc.2011.02.002 |
| [11] | Khandelwal M K, Bera P, Chakrabarti A. Influence of periodicity of sinusoidal bottom boundary condition on natural convection in porous enclosure[J].International Journal of Heat and Mass Transfer, 2012, 55(11): 2889–2900. |
| [12] | Canuto C, Hussaini M Y, Quarteroni A, et al. Spectral methods:fundamentals in single domains[M]. Berlin: Springer, 2006. |
| [13] | Peyret R. Spectral method for incompressible viscous flow[M]. New York: Springer, 2002. |
| [14] | Manole D M, Lage J L. Numerical benchmark results for natural convection in a porous medium cavity[J].ASME Publications Heat Transfer Division, 1993, 216: 55–55. |
| [15] | Moya S L, Ramos E, Sen M. Numerical study of natural convection in a tilted rectangular porous material[J].International Journal of Heat and Mass Transfer, 1987, 30(4): 741–756.DOI:10.1016/0017-9310(87)90204-3 |
| [16] | Chen Y Y, Li B W, Zhang J K. Spectral collocation method for natural convection in a square porous cavity with local thermal equilibrium and non-equilibrium models[J].International Journal of Heat and Mass Transfer, 2016, 96: 84–96.DOI:10.1016/j.ijheatmasstransfer.2016.01.007 |
