一种捕捉间断的两相氢流动模拟方法

《International Journal of Hydrogen Energy》：A discontinuity-capturing approach to the simulation of two-phase hydrogen flows

　　

本研究聚焦于在管道网络两相液氢流动背景下Godunov型方法的开发与实现，该方法在一款自主开发的代码（Arpad）中进行验证。研究动机在于克服液氢因其在低温下的极端物性而导致的数值预测热流体动力学方面的不确定性。研究人员发现，相较于成熟的可压缩流动基准问题（包括具有非线性边界条件处理的情形），本研究开发出的Godunov型方法将平均相对误差降低至1.90%。所选算例对低温管理系统具有工业相关性。获得的数值结果为两相氢流动提供了参考解。这些结果也表明，这种Godunov型方法的变体可以作为提高热流体氢动力学准确预测、支持净零航空液氢基础设施设计过程的合适解决方案。

在航空业脱碳背景下，氢气因其能够提供大量能量输出且零碳排放而备受工业界关注。航空航天领域寻求通过使用液氢推进机队来优化空间和重量，这引发了关于液氢作为能源载体影响的深入研究。然而，液氢作为一种深冷流体，其行为与标准工程流体存在本质差异，对微小热泄漏极为敏感，且与环境条件存在显著温差，导致其在管道系统中输送时通常呈现瞬态两相条件。两相现象给燃料输送系统设计带来重大挑战：泵在接近饱和条件下运行导致净正吸入压头（NPSH）几乎为零，阀门和流量计因气泡形成引发振荡和压降，换热器则需应对流动状态变化和传热系数波动。此外，沸腾和凝结前沿引起的空隙率、压力和流量急剧变化可导致系统级流动不稳定，这在预冷（管道淬火）和低温空化等过程中尤为常见。尽管现有研究在系统级架构和设计集成方面取得进展，但对液氢内部流动热流体行为的精确表征日益重要。

当前工业流程采用基于低保真度模型（0维和1维）的系统级软件设计流体管理系统，但这些工具在两相流条件下的物理行为再现能力有限。0维模型（如Siemens Simcenter Amesim）忽略动量与能量项之间的耦合效应，无法解析两相事件中流动特性的瞬态演化；1维框架虽能捕捉动量-能量耦合，但主要应用于液体火箭推进领域，其工况与民用航空存在显著差异（质量流量40-300 kg/s对比0.05-0.1 kg/s，任务时长分钟级对比小时至天级）。现有1维模拟套件中，NASA GFSSP采用基于压缩因子的状态方程，限制了对先进热物理模型的使用；EcoSimPro/ESPSS虽提供从压力基到激波捕捉的多种有限体积方法，但其Roe方法不满足熵条件，AUSM方案对接触间断的分辨率有限，且现有验证多局限于壁温测量，缺乏对流动本身的全面比较评估。学术研究中，Luchinskiy等人采用半隐式有限体积方法，Kuniyoor等人耦合逆特征线法与Crank-Nicholson格式，但均仅验证壁温数据，未评估底层流动物理。总体而言，在低温背景下将激波捕捉策略应用于1维守恒方程求解的研究相对不足，尤其缺乏针对液氢瞬态两相现象的文献证据。

鉴于此，研究人员开展了本项旨在将Godunov型方法首次应用于液氢两相流动1维降阶模型模拟的研究。Godunov型方法通过精确或近似求解Riemann问题来计算界面通量，将流动物理直接嵌入数值格式，相较于传统中心或迎风格式能更 sharp 稳健地捕捉间断，同时保持正确的波速并避免常见数值不稳定性。该方法尤其适用于液氢流动中观察到的 sharp 且常呈间断状的流动特性变化。

研究人员采用均相平衡模型（HEM）描述两相流，该模型基于三项核心假设：（1）两相处于热力学平衡状态（P_gas,T_gas=P_liq,T_liq）；（2）两相无相对运动，采用单一混合速度；（3）混合物性质按气相体积分数（α）平均。控制方程以质量、动量和能量守恒形式给出，其中动量方程包含摩擦损失源项（采用Churchill公式计算摩阻系数以适应非定常流动），能量方程包含恒定的管壁向液氢流传热项。为封闭方程组，研究人员采用基于Helmholtz自由能的多参数状态方程描述仲氢（P-H₂）热力学性质，通过CoolProp框架实现以处理临界及近临界状态的收敛问题。

数值求解方面，研究人员在自主开发的Arpad代码中实现了基于Godunov型的有限体积框架。空间离散采用一阶分段常数重构，在每个单元界面形成局部Riemann问题；采用Rusanov近似Riemann求解器计算界面通量，该求解器以 fastest 局部波速估计给出耗散性通量格式。时间推进采用前向Euler显式方法。为处理非线性源项，研究人员采用算子分裂策略，将非齐次系统分解为对流步（齐次系统求解）和源项步（常微分方程组精确求解）。边界条件处理包括：亚声速入口/出口采用Dirichlet-Neumann组合（入口给定密度和法向速度、压力外推；出口给定静压、其余变量外推），固壁边界采用速度Dirichlet条件（法向速度为零）配合其余变量Neumann条件。特别地，针对管道全截面破裂场景，研究人员发展了基于非线性寻根方法（如二分法）的声速出口边界处理：假设等焓和等熵行为，通过求解centroid与边界间的焓平衡方程，结合CoolProp数据库确定出口热力学状态，以准确预测瞬态质量流排放。

方法验证通过七个基准算例系统开展。算例1为Sod激波管问题（空气介质），验证间断捕捉能力。网格收敛性分析显示4000单元网格GCI最低，但2000单元网格在HLL和HLLC求解器下已能提供足够 sharp 的接触间断解。计算结果准确再现激波、膨胀波和接触间断的波系结构，Rusanov求解器下密度、速度、压力和内能的L₂范数相对误差分别为1.64%、0.99%、2.00%和2.75%。算例2为Fanno流动（氮气），验证稳态摩擦模型。以沿管长静压积分为关注量的GCI分析表明1000单元网格足够准确，数值解与解析解吻合良好，温度和压力的归一化L₂范数相对误差分别为0.74%和2.33%。算例3为Rayleigh流动（氮气），验证恒定传热效应。类似GCI分析确认1000单元网格适用，数值结果与解析解在温度和压力分布上均表现一致，相对误差分别为1.90%和2.42%。算例4为Fanno-Rayleigh组合流动（氮气），同时考核摩擦与传热的耦合效应。在仅Mach数分布有解析解的情况下，数值预测的压力和温度趋势与预期相符，相对误差分别为0.07%和4.1%。算例5为管道全截面破裂（甲烷），验证非定常摩擦模型。以出口质量流率时间积分为关注量的GCI分析选用4000单元网格，数值结果与工业风险评估工具PolyPiRRaM和PHAST的对比表明，Arpad能有效预测高压甲烷管道破裂后的排放率动态，尽管初期对PolyPiRRaM略有高估、对PHAST略有低估。算例6为两相水激波管，验证相变现象和波速预测。以沸腾前沿密度为关注量的GCI分析选用2000单元网格，三种Riemann求解器的结果与Daude和Galon的参考解对比确认了对凝结波和沸腾波的正确分辨，同时验证了CoolProp用于两相流的适用性。算例7为两相仲氢激波管，为本研究首次报道的P-H₂配置。2000单元网格下HLLC求解器表现最优，其一阶离散结果与Rusanov求解器二阶离散（minmod限制器）相当，但后者出现非物理振荡（违反Godunov阶障碍定理）。由于缺乏文献参考解，结果以定性物理一致性为评判标准，显示与水流案例相似的波系结构；特别指出在接近过冷液体极限（α→0）时，Riemann求解器选择对体积分数预测具有显著影响。

研究讨论部分需要特别关注的发现包括：在CoolProp调用策略上，研究人员采用温度-密度（T-ρ）直接查询而非闪蒸计算，以避免迭代寻根过程中的数值发散问题（Brent方法在小扰动下可能产生非物理压力值）；在单相区处理上，通过显式设置质量分数Q为0或1防止CoolProp返回非物理体积分数。与NIST数据库的对比验证了CoolProp对P-H₂饱和液/气密度预测的准确性（相对误差约10^-2量级）。

研究结论可概括如下：本研究提出了基于Godunov型方法的有限体积格式用于管道两相氢流动模拟，该方法采用均相平衡模型耦合CoolProp热物性库评估仲氢性质，以非定常摩擦模型描述壁面效应、以恒定源项描述传热。方法性能通过七个基准算例评估：空气激波管确认了格式的间断分辨能力；氮气Fanno、Rayleigh及组合流动验证了稳态摩擦和热传递的实现；甲烷管道破裂考核了非定常摩擦模型，与PolyPiRRaM和PHAST的风险分析工具对比显示良好预测精度；基准算例的平均L₂范数相对误差为1.90%。水和仲氢的两相激波管验证了相变捕捉能力，其中仲氢案例首次报道，其波系结构与物理预期一致，且揭示了Riemann求解器选择在接近过冷液体极限时的关键影响。这些结果为两相仲氢流动提供了首个参考数据集。研究同时指出当前框架的局限性：均相平衡假设限制其对分离流动的适用性；恒定传热项未捕捉低温管道的实际传热变化；1维简化忽略了几何和湍流影响。未来工作将扩展至变传热系数和经验沸腾关联式，并与GFSSP、ESPSS等工具进行对比验证。