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

464 真空与低温 第 31 卷 第 4 期

程中产生的已知的源项；A 为离散后速度变量的系 互摩擦项中的非线性，本文中选择将其显式处理。
数； ∆t为迭代时间步长；V i 为网格单元体积；S f 为网 在传统的计算方法中，对于式（7）和式（8）中热
格单元的各个面的面矢量；上标 t 为当前时间步， 机效应项的温度梯度定义为体心梯度。以一维均
为已知的；上标 t +∆t为下一时间步，为未知的；下 匀结构化网格为例（如图 1 所示），使用体心梯度进
标 i 为当前网格；下标 i±1为临近网格。为了处理 行离散时，获得的离散后的方程为：

T
T
T
p
t+∆t
p
p
a T t+∆t +a T t+∆t +a T i−2 +a p t+∆t +a p t+∆t +a p t+∆t = C （9）
i
i+2
i−2
i−2
i
i−2
i+2
i
i
i+2
i+2
| {z } | {z }
温度梯度项体心离散 压力梯度项体心离散

而使用面梯度离散方法的主要思想则是通过面梯
i+2 i+1 i i−1 i−2
度离散方式，使梯度项离散格式包含相邻网格的影
响，而非间隔一个网格的影响，从而避免产生体心
d
i+3/2 i+1/2 i−1/2 i−3/2 离散方法引起的数值振荡。

图 1 网格结构示意图
Fig. 1 Schematic diagram of the grid structure

由于本文的研究内容为面梯度离散方法，体心 奇数网格的解
T i+1
梯度离散格式的具体推导过程请参考附录。体心 T i−1
梯度离散后的压力泊松方程如式（9）所示，不难发 T 偶数网格的解
现体心离散方法得到的温度梯度和压力梯度项中，
T i+2
仅包含当前节点 i 以及与 i 间隔一个网格节点（即 T i T i−2
i+2 和 i-2 点）物理量的影响，而不包含节点 i 点相
邻节点的影响；相邻节点影响仅出现在其他项中。
对于一般流体（非超流氦）的流动问题中，梯度项并 位置
不占主导作用，其他项中的相邻点影响较强，此时
图 2 体心梯度导致数值震荡的原理图
体心离散方法所得到的格式数值振荡并不明显。
Fig. 2 Schematic diagram of numerical oscillations caused by
然而，在超流氦流动问题中，热机效应项的存在导
body center gradient
致梯度项的作用占主导作用 [16-17] ，此时离散方程的
数值表现主要由梯度项的离散格式决定。在梯度 1.2 热机效应项面梯度离散方法
项的主导作用下，体心离散方法间隔一个网格影响 在本节中，作者同样以一维均匀结构化网格为
的特性影响下，这个问题的离散点满足的方程组将 例，对面梯度离散方法的具体原理进行说明。对于
解耦成奇数网格点的方程组与偶数网格点的方程 式（7）和式（8）来说，其在求解区域上是处处成立的，
组。在梯度项的作用下，两组方程的解出现越来越 因此在网格面上，有如下方程：
大的差，最终发展成如图 2 所示的锯齿状数值震荡。
 
( )
∑ 
1   t t+∆t t+∆t ρ n t+∆t  
t+∆t
v =   A v −(ρ s s) ∇ f T − ∇ f p   （10）
n,f t  n,Nb n,Nb f +S n,f
A n,f  Nb(f) ρ f 
 
( )
∑ 
1   t t+∆t t+∆t ρ s t+∆t  
t+∆t
v s,f =   A s,Nb s,Nb +(ρ s s) ∇ f T − ∇ f p +S s,f   （11）
v
f
A t   Nb( f) ρ f 
s,f
式中： ∑ 为网格临近网格面上物理值的求和，以 对于离散后的连续性式（6），在一维均匀网格
Nb( f)
网格 i 为例； ∑ 为面 i±1/2以及 i±3/2上物理量的 中表示为：
Nb( f) ρ t+∆t −ρ t ∑( )
求和； ∇ f 为面上梯度，以图 1 中面“i-1/2”上的热机 V i + ρ n,f v t+∆t +ρ s,f v t+∆t S f = 0 （13）
∆t n,f s,f
效应项温度面梯度为例。面梯度可以离散为： f
将式（10）、式（11）和式（12）代入到式（13），可
T i+1 −T i
∇ i−1/ 2 T = （12）
d 以得到：