时间滤波二阶和三阶无条件能量稳定解耦算法用于变密度/黏度Cahn-Hilliard-Navier-Stokes系统

《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》：Time-filtered second- and third-order accurate in time, fully decoupled and unconditionally energy-stable schemes for variable-density phase-field model of incompressible two-phase flows

本研究针对变密度/黏度的Cahn-Hilliard-Navier-Stokes(CHNS)系统，设计新型时间滤波二阶/三阶全解耦时间离散格式，提出辅助变量重构的修正能量法，严格证明无条件能量稳定性与体积守恒性，实现高效高精度多相流模拟，发表于《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》。

　　

在复杂多相流系统中，描述两相不可压缩流体混合演化的Cahn-Hilliard-Navier-Stokes（CHNS）方程组是核心数学模型，被广泛用于模拟合金凝固、聚合物共混、微流控液滴运动等工程与生物过程。然而，当流体密度和黏度随相场函数变化时，传统数值方法面临严峻挑战——非线性耦合项导致计算复杂度激增，显式或半隐式格式受限于严格的CFL时间步长约束，而全隐式耦合求解需求巨大的存储与迭代成本；同时，如何在高阶时间离散下保持物理系统的能量耗散律与体积守恒性，成为长期未解决的难题。若无法保证能量稳定性，模拟可能出现非物理振荡或能量爆炸，失去预测可靠性；若采用低阶格式，则需极小时步维持精度，显著增加三维大规模计算耗时。这一瓶颈限制了CHNS系统在实际工程多尺度问题中的应用效率。

为攻克该难题，研究人员在《COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING》发表研究，聚焦变密度/黏度CHNS系统的“高阶时间精度”“全解耦结构”与“无条件能量稳定”三重目标，创新性地融合时间滤波技术与辅助变量能量重构策略，首次构建了时间滤波二阶及三阶完全解耦时间离散格式，并完成严格数学分析与数值验证。

关键技术包括：①设计辅助变量（Q,R,r,K）重构原始系统，引入等价修正能量泛函E_M，将非线性项转化为线性可解耦形式；②基于Backward-Euler首步预测与三节点时间滤波后处理，提升至二阶/三阶精度；③采用二阶外推Ψ^*=2Ψⁿ-Ψ^n-1处理显式项，确保全格式一致性；④压力-速度投影法解耦Stokes子系统，避免LBB条件限制；⑤离散补偿项消去能量估计中的交叉项，严格导出无条件能量衰减律与体积守恒。

研究结果分为理论框架构建、格式设计与稳定性证明三部分：

一、修正CHNS系统与连续能量律

建立辅助变量增强的等价CHNS系统（式(2.14)）：相场方程引入Q?·(u?)，化学势μ=λ(-εΔ?+sε^-1?+rε^-1(f(?)-s?))，动量方程嵌入K(ρu_t+0.5ρ_tu)、Rρ(u·?)u等辅助系数，初始条件（式(2.15)）满足u|_t=0=u⁰，?|_t=0=?⁰，r|_t=0=Q|_t=0=R|_t=0=K|_t=0=1。无外力时，系统满足能量耗散律（定理2.1）：dE_M/dt=-∫_ΩM(?)|?μ|²dx-0.5∫_Ων(?)|D(u)|²dx≤0，其中修正能量E_M=∫_Ω(0.5ρ|u|²+λε0.5|??|²+λs(2ε)^-1|?|²+λε^-1F(?))dx+α^-1(Q+R+r+K)。该重构不改原PDE解，为离散格式设计奠定基础。

二、时间滤波二阶解耦格式设计

对均匀时间网格tⁿ=nδt，分五步执行：

Step1：解相场子系统（式(3.1)-(3.2)）。预测(??ⁿ⁺¹,μ?ⁿ⁺¹)：??ⁿ⁺¹满足(??ⁿ⁺¹-?ⁿ)/δt+Q^*?·(u^*?^*)=?·(M(?^*)?μ?ⁿ⁺¹)，μ?ⁿ⁺¹=λ(-εΔ??ⁿ⁺¹+sε^-1??ⁿ⁺¹+r^*ε^-1(f(?^*)-s?^*))；更新r?ⁿ⁺¹（式(3.3)）与密度ρ?ⁿ⁺¹=0.5??ⁿ⁺¹(ρ₁-ρ₂)+0.5(ρ₁+ρ₂)。

Step2：解中间速度?ⁿ⁺¹（式(3.6)）。K^*(ρ?ⁿ⁺¹(?ⁿ⁺¹-uⁿ)/δt+0.5((ρ?ⁿ⁺¹-ρⁿ)/δt)?ⁿ⁺¹)-?·(ν(?^*)D(?ⁿ⁺¹))+?p^*+Q^*?^*?μ^*+R^*ρ?ⁿ⁺¹(u^*·?)u^*+0.5R^*?·(ρ?ⁿ⁺¹u^*)u^*+R^*J^*·?u^*+0.5R^*?·J^*u^*=ρ?ⁿ⁺¹g。

Step3：更新辅助变量（式(3.7)-(3.9)）。Q?ⁿ⁺¹、R?ⁿ⁺¹、K?ⁿ⁺¹分别由对应演化方程计算，含离散压力梯度补偿项αδt/2∫_Ω((ρ?ⁿ)^-1|?p?ⁿ|²-(ρ?ⁿ⁺¹)^-1|?p^*|²)dx以确保能量平衡。

Step4：投影步（式(3.10)-(3.11)）。解压力增量：?·(δt/ρ?ⁿ⁺¹?(p?ⁿ⁺¹-p^*))=?·?ⁿ⁺¹；修正速度?ⁿ⁺¹=?ⁿ⁺¹-(δt/ρ?ⁿ⁺¹)?(p?ⁿ⁺¹-p^*)，满足?·?ⁿ⁺¹=0。

Step5：时间滤波（式(3.12)）。对s∈{?,μ,u,p,r,Q,R,K}，滤波值sⁿ⁺¹=?ⁿ⁺¹-(1/3)(?ⁿ⁺¹-2sⁿ+s^n-1)，边界条件（式(3.13)）为u|_?Ω=0，?_n?|_?Ω=?_nμ|_?Ω=0。

三、关键性质证明

定理3.1证明格式无条件保持离散体积守恒：∫_Ω?ⁿ⁺¹dx=∫_Ω?ⁿdx（任意δt>0），源于相场方程与零通量边界的积分性质。定理3.2证无外力时离散能量E_Mⁿ⁺¹≤E_Mⁿ（任意δt>0），其中E_Mⁿ⁺¹含三项梯度平方和（??ⁿ⁺¹、?(2?ⁿ⁺¹-?ⁿ)、?(?ⁿ⁺¹-?ⁿ)）、密度加权动能项、压力投影项与辅助变量组合，通过算子A(sⁿ⁺¹)=?ⁿ⁺¹、B(sⁿ⁺¹)=?ⁿ⁺¹-sⁿ构造恒等式，消去交叉项，最终得到纯耗散型不等式。三阶格式扩展基于相同滤波框架，系数由泰勒展开约束优化。

研究表明，所提时间滤波二阶/三阶解耦格式成功解决了变密度/黏度CHNS系统的高效稳定计算难题：通过辅助变量重构与滤波校正，兼顾高阶精度与无条件能量稳定性；全解耦结构降低每步求解至独立椭圆型子问题，大幅减少计算成本；严格的体积守恒与能量衰减律保持物理真实性。该工作为多相流大规模长时间模拟提供了可靠数值工具，在材料制备、生物流体仿真等领域具有重要应用潜力。