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

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

[ t+∆t t+∆t ( ) t+∆t t+∆t ] [ t+∆t t+∆t ( ) t+∆t t+∆t ]
T −T p − p T −T p − p
ρ n,i−1/ 2 i−1 i ρ n i−1 i ρ n,i+1/ 2 i i+1 ρ n i i+1
(ρ s s) + + (ρ s s) + −
A t i−1/ 2 d ρ d A t i+1/ 2 d ρ d
n,i−1/ 2 i−1/ 2 s,i+1/ 2 i+1/ 2
[ t+∆t t+∆t ( ) t+∆t t+∆t ] [ t+∆t t+∆t ( ) t+∆t t+∆t ]
ρ s,i−1/ 2 T i−1 −T i ρ s p i−1 − p i ρ s,i+1/ 2 T i −T i+1 ρ s p i − p i+1
(ρ s s) − − (ρ s s) −
A t i−1/ 2 d ρ d A t i+1/ 2 d ρ d
s,i−1/ 2 i−1/ 2 s,i+1/ 2 i+1/ 2
( )
ρ t+∆t −ρ V P
t
= ρ f A f + （14）
S ∆t
 
∑ ∑ ∑ 
 ρ n,f ρ s,f ρ n,f ρ s,f 
ρ f A f =    A t v t+∆t + A t v t+∆t + S n,f +    （15）
 t n,Nb n,Nb t s,Nb s,Nb t t S s,f
 A A A A 
f n,f Nb(f) s,f Nb(f) n,f s,f
观察式（14）可以发现，当使用面梯度对热机效 （1）温度场的求解。通过热力学基本关系将
应项中的温度梯度进行离散时，节点 P 上的温度受 He II 的熵形式的能量方程转换为温度方程并进行
到当前网格和相邻网格点的温度的影响，这避免了 求解，从而获得新时间步的温度场 T t+∆t 。并根据求
体心梯度引起的临近网格温度对当前节点无影响， 得的温度场更新物性。
[15]
而相隔一个网格的温度却对其产生影响的问题 。 （2）动量预测。利用上一个时间步所求出的压

1.3 面梯度离散方法的算法构建
力场，代入到动量式（7）和式（8）中，获得式（16）和
本文对 He II 控制方程的求解使用 Super-PISO
式（17）：
算法 [12-13] 。该算法作为压力隐式算子分离法 [22] 的
 
分支，其核心思路仍然是通过迭代求解速度场和压 v = 1  A v −ρ s s∇T t+∆t − ρ n ∇p +S n,i（16）
 ∑


t
t
t
∗





力场。不过在原始的算法中，热机效应项中温度梯 n,i A t  k=i±1 n,k n,k ρ 
n,i
度的离散采用体心梯度。若要使用面梯度离散方  
 ∑ 
t
t
t
法进行求解，必须在原始算法的基础上进行调整。 v = 1     A v +ρ s s∇T t+∆t − ρ s ∇p +S s,i （17）

∗


s,i
s,k s,k
t 

调整后的算法流程如图 3 所示，具体流程如下： A s,i k=i±1 ρ
式中：上标“*”为通过动量预测获得的物理量。对
温度场的求解
式（16）和式（17）进行求解获得常流体和超流体的
T t+Δt
预测速度场 v 和 v 。
∗
∗
动量预测 n,P s,P
（3）物理量在面上的离散。此步骤也是面梯度
* *
v n,P , v s,P
物理量在面上的离散 离散方法中最关键的一点。在得到预测速度场 v ∗ n,P
* * ∗
v n, f , v s, f 和 v 后，使用预测速度重新构建动量方程并将预
s,P
压力方程的求解 测速度场插值到面上，得到：
t+Δt p *  
速度场的重构 方程迭代循环 v = 1  ∑ A t v ∗ −(ρ s s) ∇ f T t+∆t − ( ) ∗     
ρ n

∗


** ** n,f A t   n,Nb n,Nb f ρ f ∇ f p +S n,f 
v , v s,Pn,P n,f Nb(f)
（18）
迭代循环
判断  ( ) 
不满足 ∑  
残差要求 v = 1     A t v ∗ +(ρ s s) ∇ f T t+∆t − ρ s ∗   
∗
满足残差要求 s,f A t   s,Nb s,Nb f ρ f ∇ f p +S s,f 
s,f Nb( f)
下一个时间步的物理场 （19）
（4）压力方程的求解。将式（18）和式（19）代入
图 3 面梯度离散方法算法流程图
到式（6）中获得压力方程：
Fig. 3 Flow chart of the algorithm for face-gradient
discretization scheme
     
( )
∑ ∑  ∑  ρ n ρ s s ρ s ρ s s 
    
  ρ n,f  t ∗  ρ s,f  t ∗  t+∆t
  A v +S n,f +  A v +S s,f − − ∇ f T 


      ·S f
 t  n,Nb n,Nb t  s,Nb s,Nb 
 A   A   A t A t 
f n,f Nb(f) s,f Nb(f) n s f （20）
( ) ( )  
∑ ρ t+∆t −ρ t
 ρ n ρ n ρ s ρ s   
∗ 

=   + ∇ f p  ·S f −

  V i
 ρA t ρA t  ∆t
f n f s f
对式（20）进行求解，获得面上的压力场 p 。 （5）速度场的重构。在步骤（2）中求得的速度
∗
f