本文为国际期刊 International Journal of Heat and Mass Transfer 于2004年发表的一篇原创性研究论文，题为“使用高阶时空方法（high-order time–space method）求解立方体腔内自然对流问题的基准解”。主要作者为 Shinichiro Wakashima 与 Takeo S. Saitoh，均来自日本东北大学（Tohoku University）环境研究所。本研究旨在为三维立方体腔内的自然对流问题提供一套高精度、网格无关的基准数值解。
在计算流体力学和传热学领域，验证和评估新开发的数值方法（如新算法、新离散格式）的准确性和性能至关重要。为此，研究者们需要公认的基准问题（Benchmark Problems）及其精确解，用于对比和验证。其中，由两侧壁面存在温差驱动的二维方腔自然对流（即de Vahl Davis经典问题）是应用最广泛的基准问题之一。然而，随着工程应用（如晶体生长、金属铸造、相变储热系统等）对三维效应关注的加深，对三维自然对流基准解的需求日益迫切。尽管已有一些针对三维立方腔的研究工作，但在求解精度、网格独立性以及对高瑞利数（Rayleigh number）问题的处理效率方面，仍有待改进。特别是传统的时间推进法（time-marching method）在求解高瑞利数问题时，由于非线性和边界层变薄，需要极小时同步长和大量计算时间，效率低下。本研究的核心目标，正是利用Saitoh教授于1991年提出的高效“时空方法”（Time-Space Method, TSM），为普朗特数（Prandtl number）固定为0.71（空气）的三维立方腔自然对流问题，在瑞利数为10^4、10^5和10^6时，提供一套基于四阶空间精度和三阶时间精度的网格无关基准解。这不仅为学界提供了新的、高精度的验证标准，也进一步展示了TSM方法在处理此类问题上的强大能力。
本研究是一个纯粹的数值模拟研究，不涉及物理实验，研究对象是定义好的三维物理模型及其控制方程。主要工作流程可分为以下几个步骤：
第一步：物理模型与数学方程建立 研究模拟的对象是一个充满空气的立方体空腔（边长无量纲化为1）。两个相对的垂直壁面（x=0和x=1）保持恒定但不同的温度（Th和Tc，无量纲后为1和0），形成驱动对流的温差。其余四个壁面（顶部、底部及另外两个侧面）均为绝热（adiabatic）。重力作用方向为负z轴。流体被假设为不可压缩、层流流动，并遵循Boussinesq近似。基于这些假设，控制方程被转化为涡量-矢势（vorticity-vector potential）形式。具体包括：涡量输运方程（动量方程的无旋度形式）、能量（温度）方程、矢势泊松方程，以及联系涡量、矢势和速度的定义式。所有方程均已无量纲化。
第二步：数值方法与时空方法（TSM）的应用 这是本研究方法学的核心。与传统的、沿时间轴逐步推进求解的“时间推进法”不同，TSM将时间维度视为一个类似空间的坐标。这样，一个n维的非定常问题被转化为一个(n+1)维的定常边值问题。计算域从三维物理空间（x, y, z）扩展为四维时空域（x, y, z, t）。初始条件被视为时空域在t=0处的边界条件。TSM的最大优势在于其无条件的数值稳定性，时间步长可以任意选择，仅由时间方向上的精度或分辨率需求决定，从而避免了传统方法因稳定性要求而必须使用极小时间步长的限制，计算效率可提升数十至上百倍。在本研究中，作者将TSM应用于前述三维涡量-矢势形式的方程组求解。
第三步：高阶离散格式与边界条件处理 在四维时空域内，作者采用了高阶有限差分法进行离散： * 空间离散：采用四阶精度的中心差分格式 来离散控制方程和边界条件中的一阶及二阶空间导数。 * 时间离散：采用三阶精度的后向差分格式 来离散时间导数项。 * 边界处理：对于边界点，作者采用了特殊的处理方法。例如，对于绝热边界条件，不仅要求边界上的法向温度梯度为零，还通过引入虚拟网格点，推导出满足四阶精度的边界点及其内部虚拟点的温度表达式。对于等温壁面和无滑移速度边界、涡量边界条件，也采用了类似的高阶格式处理，以确保整体精度的一致性。所有变量在网格点上共置。
第四步：迭代求解与收敛标准 在时空域中，求解过程是一个迭代过程。作者使用逐点松弛的逐次超松弛迭代法（SOR method） 来求解离散后的代数方程组。整个时空域所有网格点上的解通过迭代同时更新。收敛判据基于涡量方程的残差L2范数（每个网格点的均方根误差），当该值小于10^-6时视为收敛。达到此标准时，温度方程和矢势方程的残差通常会更小（小1-2个数量级）。作者设置了不同的加速因子（acceleration parameter）以优化不同方程的迭代收敛速度。
第五步：网格无关性研究与基准解提取 为了获得可靠的、不依赖于网格密度的基准解，研究进行了系统的网格细化研究。针对每个瑞利数（10^4, 10^5, 10^6），分别在三种均匀网格上进行计算：40×40×40×10，80×80×80×10，以及120×120×120×10（前三个是空间网格数，最后一个是时间方向的网格层数，时间步长固定为1.0）。通过比较不同网格下的关键特征量，评估了网格依赖性。最终，采用理查德森外推法（Richardson extrapolation），基于最细的两套网格（120^3和80^3）的计算结果，外推出“网格尺寸趋于零”时的基准解。在外推中，假设离散误差的阶数为4（即所用空间离散格式的理论精度）。
1. 网格无关性分析与基准解数据： 研究表明，随着网格加密，计算得到的特征量逐渐收敛。例如，对于Ra=10^4，即使在最粗的网格上，结果也已接近收敛；而对于Ra=10^5和10^6，网格依赖性更为明显，特别是在努塞尔数（Nusselt number）等积分量上。通过理查德森外推，作者最终得到了一系列网格无关的基准解特征量，包括：腔体中心的矢势y分量（/2_center）和涡量y分量（x2_center）、水平方向最大速度（u1_max）及其位置、垂直方向最大速度（u3max）及其位置、垂直中截面（x=0.5）的平均努塞尔数（Nu{1⁄2}）、等温壁面的平均努塞尔数（Nu_{mean}）以及腔体中心的层结因子（s_center）。这些精确的数据以表格形式列出，为后续研究提供了详尽的对比标准。作者还指出，根据计算结果估算的实际离散误差阶数在3.3至5.2之间，与理论上的四阶精度基本吻合。
2. 流场与温度场特性分析： 基于最细网格（120^3）的计算结果，论文详细描绘了不同瑞利数下的流动与传热特性： * 中心截面（y=0.5平面）：温度、涡量（x2）和矢势（/2）的等值线图显示，其形态与经典的二维方腔自然对流结果高度相似。随着瑞利数增加，垂直等温壁面附近的速度边界层和热边界层显著变薄，而顶部和底部绝热壁面附近的边界层变化相对不那么剧烈。这符合高瑞利数下浮升力增强、流动加剧的物理预期。 * 速度分布：水平速度（u1）和垂直速度（u3）在中心截面的分布表明，最大速度值及其位置随瑞利数变化。一个关键的三维效应体现在横截面（如x=0.5或z=0.5平面）的速度等值线图上。在Ra=10^4时，水平速度的最大值点位于中心平面（y=0.5）上。而在Ra=10^5和10^6时，水平速度出现了两个峰值，分别靠近前、后壁面（y=0和y=1）的角落。垂直速度的峰值也显示出类似的双峰结构，并且随着瑞利数增加，这些峰值向角落移动。这表明在三维腔体中，除了中心的主流涡旋外，在侧壁（绝热壁）附近形成了螺旋状的流动管（spiral flow tubes），主要的质量和热量输运通过这些结构进行。 * 局部努塞尔数分布：局部努塞尔数在等温壁面（x=0和x=1）以及中心截面（x=0.5）上的分布进一步验证了上述流动结构。对于较低的Ra=10^4，传热最强的点位于壁面中心线附近。而对于较高的Ra=10^5和10^6，等温壁面上的局部努塞尔数也呈现出微弱的双峰趋势（虽然等值线图中不明显），而中心截面上的局部努塞尔数双峰则更为清晰。这证实了热量传递在中心区域主要由对流主导，而在靠近垂直等温壁面的边界层区域则主要由传导主导。
结论： 本研究成功地利用高阶有限差分时空方法（TSM），高效地获得了三维立方腔内自然对流在Ra=10^4, 10^5, 10^6及Pr=0.71条件下的收敛解。通过网格细化研究和理查德森外推，提供了一套具有四阶精度的、网格无关的基准解。这些解揭示了三维流动结构的细节，特别是在高瑞利数下，横截面上速度与局部努塞尔数呈现双峰分布，表明存在沿壁面发展的螺旋状流动管，这是区别于二维简化模型的重要三维特征。
科学价值与应用意义： 1. 提供高精度验证标准：所提供的一系列精确的特征量数据（如表2所示）构成了一个新的、高标准的基准测试案例，可用于评估和验证任何求解不可压缩流动和自然对流问题的数值方法（如有限体积法、有限元法、格子玻尔兹曼方法等）的性能和精度。 2. 展示高效数值方法的优越性：研究有力证明了时空方法（TSM）在求解具有高非线性、薄边界层特征的自然对流问题时的高效性和稳定性。TSM能够采用大时间步长，避免了传统时间推进法的稳定性限制，为求解类似复杂瞬态问题（如相变问题、高瑞利数湍流模拟的前期研究）提供了一种极具潜力的高效算法选择。 3. 深化对三维自然对流物理机制的理解：通过详尽的流场和温度场可视化与分析，清晰地展示了三维腔内自然对流随着瑞利数演化的流动结构和传热特性，特别是揭示了角区复杂三维二次流的存在形式，对相关工程应用中的优化设计具有理论指导意义。
论文在引言部分详细回顾了二维和三维自然对流基准解研究的发展历程，引用了de Vahl Davis、Saitoh和Hirose、Hortmann等人、Fusegi等人、Janssen等人的重要工作，为读者提供了清晰的学术脉络。此外，附录中的符号说明（Nomenclature） 非常完整，方便读者查阅。文中对边界条件的高阶离散处理方式也进行了具体描述，体现了作者在保证整体计算精度一致性方面的细致考量。最后，作者指出了前人工作中可能存在的不足（如压力边界条件假设的不妥），从侧面凸显了本工作计算的完整性和严密性（未采用任何对称性假设）。