166.6 MHz 超导腔加强筋结构的多物理场耦合仿真

2022-03-30周健荣常正则张新颖韩瑞雄朱柯宇张李少鹏

低温工程 2022年1期

关键词：加强筋圆孔气相

周健荣常正则张新颖韩瑞雄李梅朱柯宇张沛葛锐李少鹏

（1 中国科学院高能物理研究所北京 100049）

（2 中国科学院高能物理研究所射频超导与低温技术研究中心北京 100049）

（3 中国科学院大学北京 100049）

1 引言

高能同步辐射光源（High Energy Photon Source,HEPS）是由中国科学院高能物理研究所设计并建造的第四代同步辐射光源,加速器采用166.6 MHz 和499.8 MHz 超导腔对束流进行加速并为束流提供同步辐射损失的能量[1]。166.6 MHz 超导腔为电子束提供1.2 MV 加速电压和180 KW 功率,初步设计采用4.5 K 饱和液氦浸泡方式进行冷却,氦池压力为1.23×105Pa[2]。耦合器馈入功率在超导腔中建立起电磁场,由于剩余电阻作用会在超导腔壁面产生功率损耗,壁面的功率损耗如果不能及时通过液氦带走,会带来超导腔局部温度升高,一旦超过铌的超导转变温度（9.2 K）将导致整个超导腔发生失超[3-5]。

166.6MHz 超导腔的几何结构比较复杂,其内导体与小束管之间设计有加强筋以改善超导腔的应力分布,通过对两种结构的加强筋进行仿真模拟,发现在加强筋与内导体接触的边上开孔会降低加强筋的强度,使超导腔承受更大的峰值应力[6-7]。功率损耗主要发生在内导体壁面上,使用4.5 K 饱和液氦冷却时会在内导体与小束管之间的流体区域内产生气泡,气泡的排出受到加强筋结构的影响,如果气泡聚集过多会降低液氦的冷却效果,使超导腔内导体温度升高,有引发失超的风险。在加强筋中部开孔的基础上,顶部进行适当尺寸的开孔有助于气泡排除,但可能使强度降低和影响超导腔的频率稳定。如果采用2 K 超流氦进行冷却,超流氦的高导热性能将热量及时带走,不会有气泡堆积的风险[8]。因此,需要对166.6 MHz超导腔的传热机理进行系统分析,以评估优化采用4.5 K 饱和液氦进行冷却时的加强筋结构设计。

国内外对超导材料在不同应用场景下的失超过程以及用于冷却超导材料的低温流体的换热特性开展了一系列研究。Sang-ho Kimden 等建立了一个通用的热稳定性分析算法,可以对超导腔的动态热稳定性进行模拟,配合三维有限元程序,可以对超导腔结构、瞬态行为以及非线性材料特性等进行稳定性分析,运用该算法模拟了椭球形超导腔的失超过程[5]。T.Junquera 等通过机理分析和数值模拟方法分析了超导腔的热稳定性,研究了缺陷大小、射频频率、热导率、冷却条件等相关参数对超导腔失超的影响,并通过射频测量方法测量了失超过程中超导腔电磁场的变化,有助于更全面地了解腔壁上的热效应引起的失超[9]。Pashupati Dhakal 对表面经过不同方式处理的射频超导腔的热阻进行了测量,发现超导腔内壁上的损耗功率能否有效传递至液氦,铌的热导率和Kapitza 热阻是主要的影响因素,对于高RRR 值铌制成的超导腔,Kapitza 热阻对传热起决定作用[10]。Chengfeng 等提出了用于模拟LNG-HTS 新型混合能源传输管道系统的失超恢复转变模型,该模型考虑了液氮沸腾传热和超导材料电流之间的复杂影响,并使用该模型进行了电磁、热和流动分析以揭示传输线失超过程中的多物理场耦合效应[11]。Xiaoqin Zhi 等建立三维瞬态模型研究了超导故障电流限制器（RSFCLs）的局部失超和恢复过程,使用CFD 模拟了YBCO 带材失超时瞬态气泡行为、传热系数以及YBCO 带材的温度变化,并结合了高速摄像机进一步揭示了液氮冷却YBCO 带材失超和恢复过程中的气泡分布[12]。Yu Song 等通过U 形结构的液氦自然循环回路实验研究了液氦的核态沸腾热传递特性,通过控制施加的热通量依次确定了单相对流、部分核态沸腾、完全核态沸腾和薄膜沸腾4 个过程的的沸腾曲线,实验得到的传热系数与文献吻合较好。当热通量减少时,恢复过程中会出现滞后现象[13]。

以上关于超导设备冷却与传热过程的研究中,对壁面电磁功率损耗引起的液氦流动传热过程研究较少,没有考虑气泡堆积对于超导设备传热的影响。本研究采用ANSYS 软件针对166.6 MHz 超导腔内导体附近的流场进行了流动/传热/电磁多物理场耦合仿真计算,以评估、优化加强筋顶部的排气孔设计。该工作为高能同步辐射光源选择4.5 K 低温系统提供了理论支撑,也为复杂电磁结构的冷却方案评估提供了一种方法。

2 几何模型及材料说明

166.6 MHz 超导腔氦池组件及加强筋位置如图1所示。超导腔被浸泡在饱和液氦中,工作时需要保证氦池压力稳定。在模拟时重点关注小束管与内导体之间的流体区域产生的气泡能否顺利排出以及超导腔的温度分布,仿真时对原始结构进行了简化,选取了部分区域进行仿真计算,简化后的计算域如图2所示。

图1 166.6 MHz 超导腔氦池组件及加强筋位置示意图Fig.1 Schematic of 166.6 MHz superconducting cavity helium vessel and position of stifferner

图2 简化后的计算域Fig.2 Simplified computational domain

简化后的计算区域包括超导腔内导体、小束管、加强筋、部分流体区域。立方体区域的流体代替了超导腔氦池除去内导体与小束管之间区域的其余区域。由于研究的重点为内导体与小束管之间的区域,假定外部区域对内部的流动不产生影响,仅通过加强筋的圆孔进行质量交换及保证圆孔处压力与实际保持一致。立方体区域的高度与氦池直径保持一致,厚度与超导腔距氦池端部的距离保持一致,保留了内导体附近流场的特征。

原始加强筋截面为略带一些角度的圆弧,此处简化为竖直圆环,圆环中部开有8 个半径为6 mm 的圆孔。前期对顶部未开孔的加强筋进行了结构分析,其强度满足要求[7],需要评估加强筋结构对流动传热的影响。为了对比同一加强筋结构时顶部开孔大小对流场、换热效果的影响,对加强筋顶部不开孔以及开10 mm、20 mm 方孔3 种工况的流动传热进行了仿真计算和对比分析,加强筋顶部不开孔、开10 mm、20 mm方孔的结构如图3 所示,重力方向为负Z方向。为了描述方便,约定沿Y方向,两个圆孔称为一排,沿负Z方向,每排分别命名为第一、二、三、四排圆孔。

图3 加强筋结构Fig.3 Stiffener structure

超导腔材料为RRR＞300 的铌,加强筋材料为RRR=253 的铌。铌的热导率与RRR值之间的关系如图4 所示[14-15]。

图4 不同RRR 值对应的铌热导率随温度的变化Fig.4 Variation of niobium thermal conductivity corresponding to different RRR values with temperature

使用Fluent meshing 进行了网格划分,选取加强筋顶部开10 mm 方孔的模型进行验证,以换热达到稳定时EM 面最高温度作为指标,对比了5 种网格计算所得的结果,图5 为不同网格数量对应的EM 面最高温度,最终选取网格数量为810 万的网格。对其余两种模型划分网格时按照与加强筋顶部开10 mm 方孔模型相同的网格划分方式进行。

图5 不同网格对应的EM 面最高温度Fig.5 Maximum temperature of EM surface under different meshes

3 数值模型

3.1 控制方程

由于此处压力变化较小,因此可将饱和液氦及氦气视为不可压缩流体,重点关注气泡的生成和流动过程,采用VOF 两相流模型用于捕捉气泡的流动。VOF 模型通过计算控制体内每个流体的体积分数捕捉气液作用的界面,对于计算液体中大气泡的流动具有优势。

VOF 模型每个控制体中流体的物性通过每相的体积分数进行计算,式（1）,式（2）为密度与定压比热的计算公式。

式中:α为体积分数;ρ为密度,kg/m3;cp为定压比热,J/（kg·K）;下标l 为液相,v 为气相。体积分数满足式（3）:

连续性方程、动量方程和能量方程为式（4）—式（6）:

式中:为速度,m/s;μ为粘性系数,Pa·s;H为焓值,J;keff为等效热导率,W/（m·K）,为各相质量传递源项,这里代表液相转变为气相的净质量,kg;γ为潜热,J/kg。采用Lee 模型对质量传输源项进行计算,为式（7）,式（8）:

式中:re为液相的蒸发频率,1/s;rc为气相的冷凝频率,1/s;Tsat为饱和温度,K。

代表了相与相之间相互作用的体积力项,N;如果忽略相界面的剪切力,则代表了表面张力σ,N/m。采用CSF 模型计算,为式（9）:

kv为气泡的表面曲率,

式中:为气泡表面法线,为单位法线。

3.2 边界条件及计算方法

将超导腔内导体以及小束管上发生电磁损耗的面命名为EM 面,流体立方体区域上表面设置为压力出口,如图2 所示。超导腔壁面损耗与超导腔的Q0值相关,Q0值提高可以减少壁面损耗。出于安全考虑,仿真时采用的Q0值为设计值,电磁损耗较高,图6为通过HFSS 仿真计算得到的超导腔表面电磁损耗。随着超导腔Q0值的提升,壁面损耗会降低,进而减少气泡的生成。从图6 中可以看出,超导腔的表面的电磁损耗主要发生在内导体上,通过积分求得EM 面上的总电磁损耗为19.34 W,该值为保守估计,为超导腔设计的最大热损耗。将HFSS 仿真计算所得的表面电磁损耗导入Fluent 中作为热源,热量先通过超导腔壁面,再传入流体进行换热。

图6 超导腔表面电磁损耗分布Fig.6 Electromagntic loss distribution on surface of superconducting cavity

表1 给出了边界条件及初始条件设置情况,考虑了重力,初始时设置固体与流体为相同的温度。

表1 边界条件及初始条件设置Table 1 Values of boundary conditions and initial conditions settings

使用Fluent 进行计算时,选取SSTκ-ω湍流模型,压力求解器,使用SIMPLE 模型对压力速度进行耦合求解。采用了二阶迎风格式对动量、能量方程进行求解,使用Compressive 对体积分数进行求解。经过时间无关解的验证,计算开始阶段（0—5 s）将时间步长设置为0.001 s,之后将时间步调整为0.01 s,在每个时间步内,收敛判定标准为1×10-6。

4 结果与讨论

多相流计算时必须设置为瞬态计算,但实际上更关心流动稳定后的工况。因此选取流动形态和温度场基本稳定后的结果为最终计算结果。其中,加强筋顶部不开孔工况计算总时间为62 s,加强筋顶部开10 mm、20 mm 方孔工况计算总时间为48 s,图7 给出了不同工况下EM 面最高温度、平均温度随时间的变化,下文中,图a 代表加强筋顶部不开孔工况,图b 代表加强筋顶部开10 mm 方孔工况,图c 代表加强筋顶部开20 mm 方孔工况。

图7 EM 面最高温度、平均温度随时间变化Fig.7 Variation of the maximum temperature and average temperature of EM surface with time

从气相体积分数、速度、温度分布等方面对3 种工况的仿真结果进行对比分析。其中,加强筋顶部不开孔工况选取60 s 时刻的数据,顶部开10 mm、20 mm方孔两种工况选取45 s 时刻的数据。

4.1 气相体积分数分布对比

选取中心截面、截面A、B 以及流体域（负Z方向）的气相体积分数进行分析。截面A、B 位置如下图8 所示。截面A 紧贴加强筋表面,截面B 距截面A 20 mm。

图8 截面A、B 位置示意图Fig.8 Schematic diagram of location of section A and B

通过图9 和图10 可知加强筋顶部开孔与不开孔时均会在流体域的上部聚集气泡,但聚集的厚度相差较大。当换热达到稳定时,在加强筋附近,加强筋顶部开10 mm、20 mm 方孔对应的气泡聚集厚度均不超过5 mm。加强筋不开孔时气泡将聚集至加强筋第一排两个圆孔顶部,距离内导体垂直距离约30 mm。

图9 中心截面及流体域(负Z 方向)气相体积分数Fig.9 Volume fraction of central section and fluid domain (negative Z direction)

图10 截面A、B 气相体积分数Fig.10 Volume fraction of section A and B

通过截面A 处气相体积分数对比可知,加强筋顶部开孔时,仍会有部分气相从加强筋的第一、二排圆孔排出,但所占的比例较少。顶部开20 mm 方孔较开10 mm 方孔相比,气相主要集中在顶部开孔及第一排两个圆孔排出。不开孔时,气相先从第一排两个圆孔排出,随着气泡变多,第二排两个圆孔排出的气相逐渐增多。对比截面B 处气相体积分数可知,不开孔时流体下部区域中气相占据的比例较高,上部纯气相占据的区域较其他两种工况显著增多。过多的气相会导致传热性能变差,抬升内导体温度。

4.2 速度分布对比

中心截面处流场速度分布如图9 所示。

通过图11 可知,不开孔时加强筋右侧下部流体速度变化较大,加强筋顶部开10 mm、20 mm 方孔时,加强筋右侧上部流体的速度变化较大。不开孔时,气泡主要通过加强筋上部4 个圆孔流出,在靠近上部壁面附近的流体速度相对较小。加强筋顶部开孔时,产生的气泡先在顶部聚集之后从顶部开孔及加强筋第一、二排圆孔流出,所以在流体上部区域速度较大。当气泡从加强筋排出后,顶部开10 mm 方孔工况对应的流体速度分布比较集中,且最大速度均高于其余两种工况对应的最大速度。

图11 中心截面及流体域(负Z 方向)速度分布Fig.11 Speed distribution of central section and fluid domain (negative Z direction)

4.3 温度分布对比

图12 为中心截面以及内导体（负Z方向）温度分布。

图12 中心截面及内导体(负Z 方向)温度分布Fig.12 Temperature distribution of central section and inner conductor (negative Z direction)

从图7a 可知,加强筋顶部不开孔时EM 面最高温度在5 s、30 s 附近出现了较大的变化,5 s 之前,最高温度几乎成线性增长,5 s 至30 s 之间,最高温度依然增长但是增长速率逐步降低,30 s 之后,最高温度稳定在6.25 K 附近,38 s 之后温度又逐渐升高,48 s之后稳定在6.41 K。可能的原因是在5 s 之后气泡逐渐增多,由于排出不畅,流体加强了对流,温度增长逐渐变缓。30 s 之后,积聚的气泡逐渐到达加强筋第一排圆孔的顶部,气泡逐渐排出,最高温度不变约8 s 后又逐步升高,最后趋于稳定。在此期间平均温度缓慢升高,最后稳定在5.15 K。

从图7b 可知,加强筋顶部开10 mm 方孔时EM面最高温度在5 s、10 s 附近出现了较大的变化,在5 s之前,最高温度几乎成线性增长,5 s 至7.5 s 之间,最高温度增长变缓,7.5 s 至10 s 之间增长变快。可能的原因是5 s 之后气泡逐渐增多,5 s 至7.5 s 之间对流较强,换热能力提高,但是气泡没有有效排出,随着气泡的积累导致传热变差,温度增长迅速。10 s 之后,积聚的气泡通过顶部方孔以及第一排圆孔逐渐排出,最高温度逐渐趋于稳定,最后稳定在5.62 K。从平均温度看,EM 面温度平均温度在30 s 之后基本保持在4.93 K,波动很小。

从图7c 可知,加强筋顶部开20 mm 方孔时EM面最高温度在5 s、14 s 附近出现了较大的变化,在5 s之前,最高温度几乎成线性增长,5 s 至14 s 之间,最高温度增长变缓,14 s 之后逐渐稳定,换热达到稳定时,EM 面最高温度为5.59 K,平均温度为4.91 K。与7b 对比可知,开20 mm 方孔时,EM 面最高温度在前5 s 变化趋势相同,在5 s 至14 s 之间的变化趋势相差较大,换热稳定时,二者EM 面最高温度、平均温度相差较小,相对误差分别为0.5%和0.4%。

从图7a,7b,7c 对比可知,3 种工况下最高温度开始阶段前几乎线性增长,固体导热起主导作用,最高温度大约在5.5 K 后发生较大变化。加强筋顶部开10 mm、20 mm 方孔工况下,EM 面最高温度在15 s之后逐步趋于稳定,稳定在5.62 K。加强筋顶部不开孔时,EM 面最高温度在30 s 之后逐渐稳定,最后稳定在6.41 K。

通过图12 可知3 种工况下内导体最高温度的位置均位于加强筋右侧上部,内导体上部温度梯度较大,下部温度分布相对较均匀。加强筋顶部开10 mm、20 mm 方孔对应的内导体温度分布比较均匀,温度梯度较小,上部最大温差约0.8 K。加强筋顶部不开孔时内导体温度梯度较大,上部最高温差约1.5 K。在3 种工况下,内导体最高温度均处于超导转变温度9.2 K 以下。

综上可知,加强筋顶部不开孔时可以保证超导腔温度处于超导转变温区内,满足流动传热要求。其余两种工况可以改善流动传热效果,但会降低加强筋的强度。为了保证结构强度,最终选择不额外开孔的方案。