Page 83 - 《真空与低温》2026年第1期

80 真空与低温 第 32 卷 第 1 期

、
准则数，具体定义如下： 式中： q g q l为壁面传热热流密度； S lat,i为第 i 层气
、
、
gβ∆Tξ 3 体与壁面接触面积； V g V l为该层控制体积； c g,i c l,i
Gr = （10）
ν 2 分别为第 i 层的气相和液相节点对应比热容。
µc p
Pr = （11）
λ
式中：g 为重力加速度；β 为流体热膨胀系数；∆T 为
气相或液相节点与界面的温差；ξ 为计算节点距离
贮箱底部的高度；ν 为流体运动黏度；μ 为流体动力 T g,i+1
黏度；c p 为流体定压比热容。 T g,i T wg,i T fg,i
q g
贮箱总体积守恒，所以有： T g,i−1
q foam to wg
（12）
V tot = V l +V g ·
q gs
m e q lat
对上述方程进行展开与计算，可以得到气相区 q sl
q foam to wl
比内能计算公式：
˙ ˙ ( ′′ )
Q g − Q gs + ˙m e h −u g T l,i+1 q l
s
∆u g = ∆t （13）
ρ g V g
T l,i T wl,i T fl,i
根据热力学状态方程，气相的状态可由任意 T l,i−1
两个独立的状态参数确定。故采用气相的比内能
u g与密度 ρ g作为输入，通过 CoolProp 热物性数据
库 [13] 的状态函数得到气相区的压力。该开源数据
库提供了各种工质的高精度热物性数据，其核心状 图 2 贮箱集总参数数值模型计算示意图
态方程与 NIST REFPROP 兼容，计算精度具有良好 Fig. 2 Lumped-parameter coupled model of tank

的一致性。Python 开发环境能够方便地通过 Cool- 上下边界的温度演化关系则由热通量与能
Prop 提供的接口直接处理以比内能和密度为输入 量守恒共同约束，不再适用中心差分格式，需单独
的热物性数据，有效简化了真实流体精确物性计算 处理。
的实现过程，提升模型开发效率，为准确描述推 气相顶部边界（第 N 层，与贮箱前底接触）：
进剂贮箱气相区域在传热、界面蒸发和外部注气
∂T g,N dT g,N λ g,N S q g,N S
等多因素综合作用下的压力动态变化提供了保障。 = + （16）
∂t dx ρ g,N c pg,N V g,N ρ g,N V g,N c pg,N
2.1.2 温度场计算
气相底部边界（第 0 层，与气液界面相邻）：
为了计算贮箱内气液相区域流体温度场的瞬
( ′′ )
λ g,0 S q gs S
态变化过程，分别对气相区和液相区构建一维温度 ∂T g,0 = dT g,0 − + h −h g ˙m e
s
∂t dx ρ g,0 c pg,0 V g,0 ρ g,0 c pg,0 V g,0 ρ g,0 c pg,0 V g,0
场热传导模型。如图 2 所示，基于贮箱的轴对称特
（17）
性，以贮箱一侧的截面为例展开分析。对于气、液
液相顶部边界（第 N 层，与气液界面接触）：
相区域，网格均从底部向上编号，用于求解气相温 ( )
′
∂T l,N dT l,N λ l,N S q sl S h l −h ˙m e
s
度 T g,i、液相温度 T l,i。贮箱金属壁面与泡沫区域均 = + +
∂t dx ρ l,N c pl,N V l,N ρ l,N c pl,N V l,N ρ l,N c pl,N V l,N
采用与气液分层节点相对应的环形分层网格结构， （18）
以此准确捕捉泡沫、壁面与气液两相之间的热传 液相底部边界（第 0 层，与贮箱后底接触）：
递，其中气液相对应的壁面节点和泡沫区域节点对 dT l,0 λ l,0 S q l,0 S
∂T l,0
= + （19）
、
、
、
应温度记为 T wg,i T wl,i T fg,i T fl,i，壁面外与气液相 ∂t dx ρ l,0 c pl,0 V l,0 ρ l,0 c pl,0 V l,0
、 、
边界对应换热量记为 q foam to wg q foam to wl q g q l。 上述边界节点表达式考虑了边界节点与相邻
、
对于气相区或液相区中第 i 个节点有如下方程： 节点之间的热传导，以及气液相变带来的传热传质。
2 h 和 分别表示饱和液相比焓与实际液相比焓， ′′
′
d T g,i q g S lat,i ∂T g,i s h l h s
+ （14）
λ g,i = ρ g,i c g,i
dx 2 V g,i ∂t 和 h g分别表示饱和蒸汽比焓与实际气相比焓。
2
d T l,i q l S lat,i ∂T l,i 2.2 贮箱壁面热力学模型
+ （15）
λ l,i = ρ l,i c l,i
dx 2 V l,i ∂t 为了计算贮箱壁面区域温度场的瞬态变化过