Page 48 - 《真空与低温》2025年第4期

李子兀等：一种新型的超流氦二流体模型梯度项离散方法 463

discretization method can effectively suppress such oscillations. This research provides a more reliable and stable discretiza-
tion method for the numerical simulation of the superfluid helium two-fluid model，which is conducive to promoting the re-
search progress in related fields.
®
Key words：superfluid helium；OpenFOAM ；flow and heat transfer in superfluid

0 引言 震荡问题。

在大型科学装置（如大型强子对撞机、NHM- 1 算法的构建
FL 等）中，通常选择使用超流体氦（He II）来冷却磁
二流体模型将 He II 视为两种完全渗透的流体，
体，以获得更高的磁场强度 [1-2] 。研究表明，对于
分别为携带熵和黏度的常流体，以及不携带熵且完
NbTi 磁体，使用 He II 代替普通液氦（He I）作为冷
全无黏的超流体。二流体模型由方程式（1）～（4）组
[3]
却剂，可以将临界电流提高 20%～40% 。使用 He II
成，其中式（1）为连续性方程，式（2）和式（3）分别为
来冷却托卡马克装置中的超导磁体由来已久 [4-6] 。
常流体和超流体的动量方程，式（4）为能量方程。
然而，由于超流氦低温系统的复杂性，以及用常流 ∂ρ
氦（He I）或超临界氦冷却磁体即可满足设计参数， ∂t +∇(ρ n v n +ρ s v s ) = 0 （1）
这种冷却方式在随后的托卡马克中并没有继续使 ∂ρ n v n ρ n
2
+∇(ρ n v n v n ) = −ρ s s∇T − ∇p+η∇ v n − F sn
用 [7-8] 。不过，随着托卡马克运行参数的不断提高， ∂t ρ
（2）
对磁体性能提出了更高的要求：即更快的响应时间
和更高的约束能力，这也对低温系统提出了更高的 ∂ρ s v s +∇(ρ s v s v s ) = ρ s s∇T − ρ s ∇p+ F sn （3）
∂t ρ
参数要求。因此，使用 He II 而不是 He I 作为冷却
∂ A GM ρ n ρ s |v n −v s | 4
剂来冷却磁体系统被认为是提高托卡马克性能的 (ρs)+∇(ρsv n ) = （4）
∂t T
方向之一。早在 1999 年，欧洲原子能共同体就开
式中：下标 s 为超流体；下标 n 为常流体；ρ 为密度；
展了 LCT 项目，成功证明了在下一代聚变装置中
v 为速度；T 为温度；p 为压力；η 为黏度；s 为熵。动
[9]
使用 He II 冷却大型超导磁体线圈的可行性 。此
量方程中的 ρ s s∇T 为热机效应项。He II 微观上的
外，在下一代聚变装置—中国聚变工程试验堆
量子湍流效应，在宏观上表现为互摩擦项。在本文
（CFETR）中，也考虑用 1.8 K 超流体氦来冷却超导 [19-21]
中的形式被设定为 ：
[10]
磁体，从而获得更高的参数 。
2
F sn = A GM ρ n ρ s |v n − v s | (v n −v s ) （5）
然而，对于 He II 的流动传热研究难度相当之
式中：A G 为互摩擦系数。
M
大，其中一个原因就在于对目前广泛被接受的 He II
1.1 方程的离散
物理模型—二流体模型 的求解极其困难。为了
[11]
®
由于 OpenFOAM 平台使用的为有限体积法对
对 二 流 体 模 型 进 行 求 解， Soulaine 等 [12-13] 提 出 了
控制方程进行求解，对控制方程在时间上使用欧拉
Super-PISO 算法，并利用该算法对 He II 的热对流、
格式离散，空间上使用中心差分格式可以得到：
迫流传热以及 He II 在多孔介质中的流动传热现象
ρ t+∆t −ρ t ∑( )
进行了研究。除此之外，该算法还成功捕捉到了 V i + ρ n,f v t+∆t +ρ s,k v t+∆t ·S f = 0 （6）
∆t n,f s,f
在 He II 圆柱热对流中由量子湍流导致的独特的四 f
[14]
涡结构 。然而，该算法将热机效应项中的温度梯 式中：f 为 f = i±1/2表示网格面上物理量的离散值，
度定义为体心梯度，这在某些情况下会引起数值求 故 ∑ 为网格面上物理量的求和。同时常流体动量
f
[15]
解的锯齿形震荡 。尤其在对 He II 进行数值模拟 方程可以离散为：
时，热机效应项的数值在整个方程组中占比较 t t+∆t ∑ t t+∆t t+∆t ρ n t+∆t
A v = A v −ρ s s∇T − ∇p +S n,i （7）
n,k n,k
n,i n,i
大 [16-17] ，温度梯度定义为体心梯度更容易引发数值 k=i±1 ρ
震荡。为了解决该问题，本文基于交错网格的思 与之相似，超流体动量方程可以离散为：
[18]
路 ，提出一种二流体模型温度梯度项的面梯度离 t t+∆t ∑ t t+∆t t+∆t ρ s t+∆t
A v = A v +ρ s s∇T − ∇p +S s,i （8）
s,i s,i
s,k s,k
®
散方法，并将该离散方法植入到 OpenFOAM ，用于 k=i±1 ρ
解决在 He II 数值模拟中由热机效应项导致的数值 式中：S n 和 S s 为常流体和超流体动量方程离散过