学术研究报告：高径比对GaAs熔体液桥热毛细对流失稳的影响
一、 研究作者、机构与发表信息
本研究的主要作者为周游、曾忠（通讯作者）、刘浩和张良奇。其中，周游、曾忠、张良奇来自重庆大学航空航天学院，刘浩来自重庆交通大学西南水运工程科学研究所。该研究成果以题为“高径比对GaAs熔体液桥热毛细对流失稳的影响”的论文形式，发表于《力学学报》（Chinese Journal of Theoretical and Applied Mechanics）2022年2月第54卷第2期。该研究获得了国家自然科学基金、中央高校基本科研业务费及重庆市教委科学技术研究项目的资助。
二、 学术背景与研究目的
本研究属于流体力学，特别是微重力流体物理与晶体生长流体动力学交叉领域。研究背景源于浮区法晶体生长技术，这是一种在空间微重力环境下生产高品质单晶的关键技术。在此环境中，重力驱动的浮力对流被极大抑制，由温度梯度引起的表面张力梯度驱动的热毛细对流成为熔体中物质与热量传输的主导机制。液桥模型作为浮区法的简化模型，是研究热毛细对流的经典物理模型。
热毛细对流从稳定的二维轴对称流动向三维不稳定流动转变的临界条件（失稳）及其内在机制，是保障晶体生长过程稳定、控制晶体质量的核心科学问题。前人研究表明，流体的普朗特数是决定失稳模式（静态失稳或振荡失稳）和失稳机制（热毛细机制主导或水动力学惯性机制主导）的关键参数。例如，典型低普朗特数熔体（如Pr=0.011）通常发生静态失稳，由水动力学惯性机制主导；而典型高普朗特数熔体（如Pr>1）则通常发生振荡失稳，由热毛细机制主导。然而，对于普朗特数处于中间范围的熔体，特别是GaAs熔体（Pr=0.068），其失稳行为如何，尤其是几何参数（如液桥的高径比）对其失稳模式和机制的影响，尚缺乏系统深入的研究。已有的研究多集中于高或低普朗特数熔体，对Pr=0.068这类“中间”普朗特数熔体在高径比变化下的行为关注不足。
因此，本研究的目的是：针对GaAs熔体（Pr=0.068），系统研究液桥高径比对其热毛细对流从稳态到失稳转变的临界条件、失稳模式（振荡或静态）以及内在物理机制的影响。研究旨在填补这一参数区间内，几何效应如何与流体物性耦合影响失稳行为的认知空白，为相关晶体生长过程的精确控制和优化提供理论依据。
三、 详细研究流程与方法
本研究采用数值模拟与理论分析相结合的方法，核心是基于谱元法的线性稳定性分析与扰动能量分析。整个工作流程可分为以下几个详细步骤：
建立物理与数学模型：
- 物理模型： 研究采用经典的圆柱形液桥模型。液桥被约束在两个同轴、不同温度的圆盘之间，上端为低温冷壁，下端为高温热壁。假设熔体为不可压缩牛顿流体，表面张力随温度线性变化，自由表面为绝热且形状固定（忽略变形）。核心几何参数为高径比（Aspect Ratio, As），定义为液柱高度H与半径R之比（As = H/R）。本研究系统地改变了高径比，范围从0.4到2.5。
- 控制方程与无量纲化： 基于质量、动量和能量守恒定律，建立了描述流体流动与传热的控制方程组（纳维-斯托克斯方程与能量方程）。研究采用特征尺度（如高度H、运动粘度ν等）对控制方程进行无量纲化，引入了两个关键的无量纲参数：普朗特数（Pr，表征动量扩散与热扩散的相对速率）和Marangoni数（Ma，表征热毛细驱动力的强度，是本研究中的关键控制参数，其临界值Mac即为失稳临界点）。
- 边界条件： 在固壁（冷、热圆盘）上施加无滑移速度和等温边界条件；在自由表面上，施加切向应力平衡条件（即表面张力梯度与粘性应力平衡），法向速度为零，以及绝热条件。
求解基态与线性稳定性分析：
- 基态求解： 失稳前的流动是二维轴对称且定常的，称为“基态”。研究团队采用基于时间分裂法的谱元法这一高精度数值方法，直接求解无量纲化后的定常控制方程，获得了不同高径比下、不同Ma数对应的二维轴对称基态流场和温度场。谱元法结合了有限元法在几何处理上的灵活性和谱方法的高精度，通过在每个单元内使用高阶多项式（本研究采用6阶）逼近解，确保了计算精度。
- 引入扰动与线性化： 在获得的基态解上，叠加一个三维、无限小的扰动（包括速度、压力和温度扰动）。将叠加了扰动的场代入控制方程，并忽略扰动的高阶项（即进行线性化），得到关于扰动量的线性控制方程组。
- 正则模分析与特征值问题： 假设扰动量具有沿周向传播的行波形式，即正比于 e^(λt + ikθ)，其中k是周向波数，λ是复数增长率（λ = λ_r + iλ_i）。λ_r > 0 表示扰动增长（流动不稳定），λ_r < 0 表示扰动衰减（流动稳定），λ_r = 0 即为临界稳定状态。λ_i 对应扰动的振荡频率，λ_i = 0 对应静态失稳（转向三维定常流），λ_i ≠ 0 对应振荡失稳（转向三维振荡流）。
- 求解广义特征值问题： 将正则模形式的扰动代入线性化方程，利用谱元法进行空间离散，最终将问题转化为一个大型的广义特征值问题 A x̃ = λ B x̃。研究者使用开源的ARPACK程序包（基于隐式重启Arnoldi方法）高效求解此特征值问题，从而得到对于给定高径比和Ma数下的特征值λ和特征向量x̃（即扰动模态）。
- 确定临界参数： 通过扫描不同的Ma数，找到使最大λ_r从负变为正的临界Ma数，即临界Marangoni数（Mac）。同时，记录下该临界状态对应的临界波数（kc）和临界频率（fc = λ_i / (2π)）。研究中采用了插值公式来精确确定Mac和fc。
扰动能量分析揭示失稳机制：
- 为深入理解失稳的物理根源，研究进行了扰动动能平衡分析。通过对扰动动量方程进行一系列数学操作，并积分到整个计算域，得到了扰动动能增长率的分解表达式。
- 该表达式将扰动动能的来源分解为若干项，其中关键的两部分是：
- Iv项（惯性项）： 表示从基态流场到扰动流场的动能传递率，反映了水动力学惯性机制的贡献。此项可进一步分解为Iv1至Iv5五个子项，分别对应基态速度不同分量的梯度与扰动速度的相互作用。
- M项（热毛细项）： 表示由扰动引起的热毛细力（源于扰动温度场导致的表面张力变化）沿轴向（Mz）和周向（Mφ）做功的功率，反映了热毛细机制的贡献。
- 通过计算在临界状态下各项能量贡献的占比，可以定量判断是哪种物理机制在驱动失稳过程中起主导作用。
网格依赖性验证与程序有效性检验：
- 为确保计算结果可靠，研究进行了网格无关性验证。针对Pr=0.068, As=1的案例，比较了不同网格分辨率（单元内采用5阶、6阶、7阶谱离散）下计算得到的临界参数（Mac, kc），结果表明结果收敛，最终选用6阶谱离散进行后续所有计算。
- 通过将自主开发的谱元法程序的计算结果与已发表的权威文献结果进行对比，验证了程序的正确性。例如，将As=1, Pr=0.068和Pr=0.183的计算结果（转化为临界雷诺数Rec=Mac/Pr）与Levenstam等人（2001）的结果对比，相对误差很小（Rec误差<1.8%，fc误差<1.5%）。同时，对低Pr熔体（Pr=0.009）在不同高径比下的计算结果也与Li等人（2016）的线性稳定性分析和直接数值模拟结果吻合良好。
四、 主要研究结果
研究系统考察了高径比（As从0.4到2.5）对GaAs熔体（Pr=0.068）液桥热毛细对流失稳的影响，主要结果如下：
失稳模式随高径比发生转变：
- 研究发现，GaAs熔体的失稳模式并非固定，而是强烈依赖于高径比。这区别于典型的低Pr或高Pr熔体。
- 当 0.4 ≤ As ≤ 1.18 时，热毛细对流失稳表现为振荡失稳。即当Marangoni数超过临界值Mac后，流动从二维轴对称定常态转变为三维周期性振荡对流。临界频率fc不为零。
- 当 1.20 ≤ As ≤ 2.5 时，热毛细对流失稳表现为静态失稳。即失稳后流动转变为三维定常对流。临界频率fc为零。
- 这一发现清晰地表明，对于Pr=0.068的熔体，仅凭普朗特数无法唯一确定其失稳模式，几何参数高径比起到了关键的调控作用。
临界参数随高径比的变化：
- 临界Marangoni数（Mac）整体上随高径比As的增加而下降，意味着更“瘦高”的液桥在更小的温差（或更弱的热毛细力）下就容易失稳。
- 在Mac-As关系曲线上，于As≈0.67和As≈1.18附近出现了两个局部峰值，表明在这些特定高径比下，流动相对更稳定（需要更高的Ma才能失稳）。
- 失稳的临界波数（kc）随As增加而减小，从As=0.4时的kc=6，逐渐降至As=1.18时的kc=2，在As≥1.20后进一步降至kc=1。这符合“kc * As ≈ 常数”的物理直观，即液桥周长所能容纳的涡胞数目与其高度（或说轴向尺度）成反比。
失稳机制的能量分析：
- 扰动能量分析揭示了失稳的物理机制。对于GaAs熔体（Pr=0.068），其失稳同时包含水动力学惯性机制和热毛细机制的贡献，但以前者为主导。
- 水动力学惯性机制的贡献主要体现在能量项Iv4和Iv5上，这两项代表了基态轴向速度的径向梯度（∂w0/∂r）和轴向梯度（∂w0/∂z）向扰动动能传递能量的过程。
- 热毛细机制的贡献体现在能量项Mz和Mφ上，代表了扰动温度场引起的表面张力梯度做功。
- 两种机制对总扰动动能增长的贡献比例随高径比变化：
- 在振荡失稳区域（As ≤ 1.18），Iv4和Iv5的贡献占主导（约75%-104%），Mz+Mφ的贡献占比较小（约4%-27%），但在As接近1.18时，热毛细机制的贡献显著增加。
- 在静态失稳区域（As ≥ 1.20），依然是Iv4和Iv5占主导（约71%-126%），但Mz+Mφ的贡献比例呈现先增后减再增的非单调变化。特别地，在As较大（如1.67-2.22）时，周向热毛细力做功项Mφ甚至变为负值，起到了抑制失稳的作用。
- 研究还通过可视化临界状态下的扰动温度场和速度场，直观展示了不同高径比下失稳涡结构（如温度“热极”与“冷极”的分布）和扰动流态的差异，这些流场结构与能量分析结果相互印证。例如，在静态失稳模式下，自由液面扰动速度的方向与典型低Pr熔体有所不同，部分区域出现了从冷极指向热极的反向流动，这反映了热毛细机制参与的复杂性。
五、 研究结论与意义
本研究得出以下核心结论： 1. 对于普朗特数为0.068的GaAs熔体，其液桥热毛细对流的失稳模式（振荡或静态）并非由物性参数唯一决定，而是受到液桥高径比的显著调控。在As ≤ 1.18时为振荡失稳，在As ≥ 1.20时为静态失稳。 2. 在整个研究的高径比范围内（0.4 ≤ As ≤ 2.5），GaAs熔体热毛细对流的失稳主要由水动力学惯性机制驱动，但热毛细机制也始终参与并贡献能量，两者贡献比随高径比动态变化。 3. 临界Marangoni数（Mac）总体随高径比增大而减小，但在As ≈ 0.67和1.18处存在局部极大值，表明在这些几何比例下流动稳定性相对较高。
本研究的科学价值在于： * 深化了对中间普朗特数流体热毛细失稳机理的理解： 明确了在Pr=0.068附近，失稳是惯性机制与热毛细机制共同作用的“混合模式”，且几何形状可以改变两者的竞争关系，甚至诱发失稳模式的转变。 * 揭示了几何参数对失稳行为的调控作用： 系统阐明了高径比不仅影响失稳的临界条件（Mac, kc），还能根本性地改变失稳的模式（振荡/静态），这为通过设计晶体生长系统的几何尺寸来控制流动稳定性提供了理论指导。 * 提供了完整的参数空间图谱： 研究给出了GaAs熔体在宽高径比范围内的失稳临界参数、模式及机制图谱，填补了该领域的研究空白。
其应用价值主要体现在为空间浮区法生长GaAs等化合物半导体晶体的工艺优化提供关键理论参考。通过选择合适的高径比，可以主动调控熔体对流的失稳阈值和模态，从而抑制有害的振荡或三维流动，促进生长界面稳定，最终提高晶体的均匀性和质量。
六、 研究亮点
- 重要的新发现： 首次系统揭示并报告了对于GaAs熔体（Pr=0.068），液桥热毛细对流的失稳模式会随高径比发生从“振荡失稳”到“静态失稳”的转变。这一发现挑战了仅凭普朗特数划分失稳模式的传统认知，突出了几何效应的重要性。
- 机理阐释的深度： 不仅给出了宏观的失稳临界参数和模式，更通过深入的扰动能量分析，定量揭示了失稳过程中水动力学惯性机制与热毛细机制各自贡献的比例及其随高径比的变化规律，从能量输运的角度阐明了失稳的物理本质。
- 研究方法的可靠性： 采用了高精度的谱元法进行基态求解和线性稳定性分析，并进行了严格的网格无关性验证和与经典文献结果的对比，确保了研究结论的可靠性。
- 研究对象的针对性： 聚焦于具有重要应用背景（半导体晶体生长）但研究相对较少的中间普朗特数熔体（GaAs, Pr=0.068），探讨了高径比这一关键工程参数的影响，具有明确的工程应用指向性。
七、 其他有价值的内容
研究中对不同高径比区间内失稳机制的细分讨论颇具价值。例如，指出在振荡失稳区内，主导的惯性项会从Iv4和Iv5共同主导（小As时），转变为Iv5主导（中等As时），再转变为Iv4主导（大As时）。在静态失稳区内，当As很大时，周向热毛细力项Mφ会转变为稳定项。这些细节表明，即使同一种主导机制（水动力学惯性），其内部具体的能量传递路径（即哪个速度梯度项起主要作用）也会因几何形状而异，这为更精细化地理解和预测流动不稳定性提供了更深入的视角。