166.6 MHz 超導(dǎo)腔加強筋結(jié)構(gòu)的多物理場耦合仿真

周健榮 常正則 張新穎 韓瑞雄 李 梅 朱柯宇張 沛 葛 銳 李少鵬

（1 中國科學(xué)院高能物理研究所 北京 100049）

（2 中國科學(xué)院高能物理研究所射頻超導(dǎo)與低溫技術(shù)研究中心 北京 100049）

（3 中國科學(xué)院大學(xué) 北京 100049）

1 引 言

高能同步輻射光源（High Energy Photon Source,HEPS）是由中國科學(xué)院高能物理研究所設(shè)計并建造的第四代同步輻射光源,加速器采用166.6 MHz 和499.8 MHz 超導(dǎo)腔對束流進行加速并為束流提供同步輻射損失的能量[1]。166.6 MHz 超導(dǎo)腔為電子束提供1.2 MV 加速電壓和180 KW 功率,初步設(shè)計采用4.5 K 飽和液氦浸泡方式進行冷卻,氦池壓力為1.23×105Pa[2]。耦合器饋入功率在超導(dǎo)腔中建立起電磁場,由于剩余電阻作用會在超導(dǎo)腔壁面產(chǎn)生功率損耗,壁面的功率損耗如果不能及時通過液氦帶走,會帶來超導(dǎo)腔局部溫度升高,一旦超過鈮的超導(dǎo)轉(zhuǎn)變溫度（9.2 K）將導(dǎo)致整個超導(dǎo)腔發(fā)生失超[3-5]。

166.6MHz 超導(dǎo)腔的幾何結(jié)構(gòu)比較復(fù)雜,其內(nèi)導(dǎo)體與小束管之間設(shè)計有加強筋以改善超導(dǎo)腔的應(yīng)力分布,通過對兩種結(jié)構(gòu)的加強筋進行仿真模擬,發(fā)現(xiàn)在加強筋與內(nèi)導(dǎo)體接觸的邊上開孔會降低加強筋的強度,使超導(dǎo)腔承受更大的峰值應(yīng)力[6-7]。功率損耗主要發(fā)生在內(nèi)導(dǎo)體壁面上,使用4.5 K 飽和液氦冷卻時會在內(nèi)導(dǎo)體與小束管之間的流體區(qū)域內(nèi)產(chǎn)生氣泡,氣泡的排出受到加強筋結(jié)構(gòu)的影響,如果氣泡聚集過多會降低液氦的冷卻效果,使超導(dǎo)腔內(nèi)導(dǎo)體溫度升高,有引發(fā)失超的風(fēng)險。在加強筋中部開孔的基礎(chǔ)上,頂部進行適當(dāng)尺寸的開孔有助于氣泡排除,但可能使強度降低和影響超導(dǎo)腔的頻率穩(wěn)定。如果采用2 K 超流氦進行冷卻,超流氦的高導(dǎo)熱性能將熱量及時帶走,不會有氣泡堆積的風(fēng)險[8]。因此,需要對166.6 MHz超導(dǎo)腔的傳熱機理進行系統(tǒng)分析,以評估優(yōu)化采用4.5 K 飽和液氦進行冷卻時的加強筋結(jié)構(gòu)設(shè)計。

國內(nèi)外對超導(dǎo)材料在不同應(yīng)用場景下的失超過程以及用于冷卻超導(dǎo)材料的低溫流體的換熱特性開展了一系列研究。Sang-ho Kimden 等建立了一個通用的熱穩(wěn)定性分析算法,可以對超導(dǎo)腔的動態(tài)熱穩(wěn)定性進行模擬,配合三維有限元程序,可以對超導(dǎo)腔結(jié)構(gòu)、瞬態(tài)行為以及非線性材料特性等進行穩(wěn)定性分析,運用該算法模擬了橢球形超導(dǎo)腔的失超過程[5]。T.Junquera 等通過機理分析和數(shù)值模擬方法分析了超導(dǎo)腔的熱穩(wěn)定性,研究了缺陷大小、射頻頻率、熱導(dǎo)率、冷卻條件等相關(guān)參數(shù)對超導(dǎo)腔失超的影響,并通過射頻測量方法測量了失超過程中超導(dǎo)腔電磁場的變化,有助于更全面地了解腔壁上的熱效應(yīng)引起的失超[9]。Pashupati Dhakal 對表面經(jīng)過不同方式處理的射頻超導(dǎo)腔的熱阻進行了測量,發(fā)現(xiàn)超導(dǎo)腔內(nèi)壁上的損耗功率能否有效傳遞至液氦,鈮的熱導(dǎo)率和Kapitza 熱阻是主要的影響因素,對于高RRR 值鈮制成的超導(dǎo)腔,Kapitza 熱阻對傳熱起決定作用[10]。Chengfeng 等提出了用于模擬LNG-HTS 新型混合能源傳輸管道系統(tǒng)的失超恢復(fù)轉(zhuǎn)變模型,該模型考慮了液氮沸騰傳熱和超導(dǎo)材料電流之間的復(fù)雜影響,并使用該模型進行了電磁、熱和流動分析以揭示傳輸線失超過程中的多物理場耦合效應(yīng)[11]。Xiaoqin Zhi 等建立三維瞬態(tài)模型研究了超導(dǎo)故障電流限制器（RSFCLs）的局部失超和恢復(fù)過程,使用CFD 模擬了YBCO 帶材失超時瞬態(tài)氣泡行為、傳熱系數(shù)以及YBCO 帶材的溫度變化,并結(jié)合了高速攝像機進一步揭示了液氮冷卻YBCO 帶材失超和恢復(fù)過程中的氣泡分布[12]。Yu Song 等通過U 形結(jié)構(gòu)的液氦自然循環(huán)回路實驗研究了液氦的核態(tài)沸騰熱傳遞特性,通過控制施加的熱通量依次確定了單相對流、部分核態(tài)沸騰、完全核態(tài)沸騰和薄膜沸騰4 個過程的的沸騰曲線,實驗得到的傳熱系數(shù)與文獻吻合較好。當(dāng)熱通量減少時,恢復(fù)過程中會出現(xiàn)滯后現(xiàn)象[13]。

以上關(guān)于超導(dǎo)設(shè)備冷卻與傳熱過程的研究中,對壁面電磁功率損耗引起的液氦流動傳熱過程研究較少,沒有考慮氣泡堆積對于超導(dǎo)設(shè)備傳熱的影響。本研究采用ANSYS 軟件針對166.6 MHz 超導(dǎo)腔內(nèi)導(dǎo)體附近的流場進行了流動/傳熱/電磁多物理場耦合仿真計算,以評估、優(yōu)化加強筋頂部的排氣孔設(shè)計。該工作為高能同步輻射光源選擇4.5 K 低溫系統(tǒng)提供了理論支撐,也為復(fù)雜電磁結(jié)構(gòu)的冷卻方案評估提供了一種方法。

2 幾何模型及材料說明

166.6 MHz 超導(dǎo)腔氦池組件及加強筋位置如圖1所示。超導(dǎo)腔被浸泡在飽和液氦中,工作時需要保證氦池壓力穩(wěn)定。在模擬時重點關(guān)注小束管與內(nèi)導(dǎo)體之間的流體區(qū)域產(chǎn)生的氣泡能否順利排出以及超導(dǎo)腔的溫度分布,仿真時對原始結(jié)構(gòu)進行了簡化,選取了部分區(qū)域進行仿真計算,簡化后的計算域如圖2所示。

圖1 166.6 MHz 超導(dǎo)腔氦池組件及加強筋位置示意圖Fig.1 Schematic of 166.6 MHz superconducting cavity helium vessel and position of stifferner

圖2 簡化后的計算域Fig.2 Simplified computational domain

簡化后的計算區(qū)域包括超導(dǎo)腔內(nèi)導(dǎo)體、小束管、加強筋、部分流體區(qū)域。立方體區(qū)域的流體代替了超導(dǎo)腔氦池除去內(nèi)導(dǎo)體與小束管之間區(qū)域的其余區(qū)域。由于研究的重點為內(nèi)導(dǎo)體與小束管之間的區(qū)域,假定外部區(qū)域?qū)?nèi)部的流動不產(chǎn)生影響,僅通過加強筋的圓孔進行質(zhì)量交換及保證圓孔處壓力與實際保持一致。立方體區(qū)域的高度與氦池直徑保持一致,厚度與超導(dǎo)腔距氦池端部的距離保持一致,保留了內(nèi)導(dǎo)體附近流場的特征。

原始加強筋截面為略帶一些角度的圓弧,此處簡化為豎直圓環(huán),圓環(huán)中部開有8 個半徑為6 mm 的圓孔。前期對頂部未開孔的加強筋進行了結(jié)構(gòu)分析,其強度滿足要求[7],需要評估加強筋結(jié)構(gòu)對流動傳熱的影響。為了對比同一加強筋結(jié)構(gòu)時頂部開孔大小對流場、換熱效果的影響,對加強筋頂部不開孔以及開10 mm、20 mm 方孔3 種工況的流動傳熱進行了仿真計算和對比分析,加強筋頂部不開孔、開10 mm、20 mm方孔的結(jié)構(gòu)如圖3 所示,重力方向為負(fù)Z方向。為了描述方便,約定沿Y方向,兩個圓孔稱為一排,沿負(fù)Z方向,每排分別命名為第一、二、三、四排圓孔。

圖3 加強筋結(jié)構(gòu)Fig.3 Stiffener structure

超導(dǎo)腔材料為RRR＞300 的鈮,加強筋材料為RRR=253 的鈮。鈮的熱導(dǎo)率與RRR值之間的關(guān)系如圖4 所示[14-15]。

圖4 不同RRR 值對應(yīng)的鈮熱導(dǎo)率隨溫度的變化Fig.4 Variation of niobium thermal conductivity corresponding to different RRR values with temperature

使用Fluent meshing 進行了網(wǎng)格劃分,選取加強筋頂部開10 mm 方孔的模型進行驗證,以換熱達到穩(wěn)定時EM 面最高溫度作為指標(biāo),對比了5 種網(wǎng)格計算所得的結(jié)果,圖5 為不同網(wǎng)格數(shù)量對應(yīng)的EM 面最高溫度,最終選取網(wǎng)格數(shù)量為810 萬的網(wǎng)格。對其余兩種模型劃分網(wǎng)格時按照與加強筋頂部開10 mm 方孔模型相同的網(wǎng)格劃分方式進行。

圖5 不同網(wǎng)格對應(yīng)的EM 面最高溫度Fig.5 Maximum temperature of EM surface under different meshes

3 數(shù)值模型

3.1 控制方程

由于此處壓力變化較小,因此可將飽和液氦及氦氣視為不可壓縮流體,重點關(guān)注氣泡的生成和流動過程,采用VOF 兩相流模型用于捕捉氣泡的流動。VOF 模型通過計算控制體內(nèi)每個流體的體積分?jǐn)?shù)捕捉氣液作用的界面,對于計算液體中大氣泡的流動具有優(yōu)勢。

VOF 模型每個控制體中流體的物性通過每相的體積分?jǐn)?shù)進行計算,式（1）,式（2）為密度與定壓比熱的計算公式。

式中:α為體積分?jǐn)?shù);ρ為密度,kg/m3;cp為定壓比熱,J/（kg·K）;下標(biāo)l 為液相,v 為氣相。體積分?jǐn)?shù)滿足式（3）:

連續(xù)性方程、動量方程和能量方程為式（4）—式（6）:

式中:為速度,m/s;μ為粘性系數(shù),Pa·s;H為焓值,J;keff為等效熱導(dǎo)率,W/（m·K）,為各相質(zhì)量傳遞源項,這里代表液相轉(zhuǎn)變?yōu)闅庀嗟膬糍|(zhì)量,kg;γ為潛熱,J/kg。采用Lee 模型對質(zhì)量傳輸源項進行計算,為式（7）,式（8）:

式中:re為液相的蒸發(fā)頻率,1/s;rc為氣相的冷凝頻率,1/s;Tsat為飽和溫度,K。

代表了相與相之間相互作用的體積力項,N;如果忽略相界面的剪切力,則代表了表面張力σ,N/m。采用CSF 模型計算,為式（9）:

kv為氣泡的表面曲率,

式中:為氣泡表面法線,為單位法線。

3.2 邊界條件及計算方法

將超導(dǎo)腔內(nèi)導(dǎo)體以及小束管上發(fā)生電磁損耗的面命名為EM 面,流體立方體區(qū)域上表面設(shè)置為壓力出口,如圖2 所示。超導(dǎo)腔壁面損耗與超導(dǎo)腔的Q0值相關(guān),Q0值提高可以減少壁面損耗。出于安全考慮,仿真時采用的Q0值為設(shè)計值,電磁損耗較高,圖6為通過HFSS 仿真計算得到的超導(dǎo)腔表面電磁損耗。隨著超導(dǎo)腔Q0值的提升,壁面損耗會降低,進而減少氣泡的生成。從圖6 中可以看出,超導(dǎo)腔的表面的電磁損耗主要發(fā)生在內(nèi)導(dǎo)體上,通過積分求得EM 面上的總電磁損耗為19.34 W,該值為保守估計,為超導(dǎo)腔設(shè)計的最大熱損耗。將HFSS 仿真計算所得的表面電磁損耗導(dǎo)入Fluent 中作為熱源,熱量先通過超導(dǎo)腔壁面,再傳入流體進行換熱。

圖6 超導(dǎo)腔表面電磁損耗分布Fig.6 Electromagntic loss distribution on surface of superconducting cavity

表1 給出了邊界條件及初始條件設(shè)置情況,考慮了重力,初始時設(shè)置固體與流體為相同的溫度。

表1 邊界條件及初始條件設(shè)置Table 1 Values of boundary conditions and initial conditions settings

使用Fluent 進行計算時,選取SSTκ-ω湍流模型,壓力求解器,使用SIMPLE 模型對壓力速度進行耦合求解。采用了二階迎風(fēng)格式對動量、能量方程進行求解,使用Compressive 對體積分?jǐn)?shù)進行求解。經(jīng)過時間無關(guān)解的驗證,計算開始階段（0—5 s）將時間步長設(shè)置為0.001 s,之后將時間步調(diào)整為0.01 s,在每個時間步內(nèi),收斂判定標(biāo)準(zhǔn)為1×10-6。

4 結(jié)果與討論

多相流計算時必須設(shè)置為瞬態(tài)計算,但實際上更關(guān)心流動穩(wěn)定后的工況。因此選取流動形態(tài)和溫度場基本穩(wěn)定后的結(jié)果為最終計算結(jié)果。其中,加強筋頂部不開孔工況計算總時間為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

從氣相體積分?jǐn)?shù)、速度、溫度分布等方面對3 種工況的仿真結(jié)果進行對比分析。其中,加強筋頂部不開孔工況選取60 s 時刻的數(shù)據(jù),頂部開10 mm、20 mm方孔兩種工況選取45 s 時刻的數(shù)據(jù)。

4.1 氣相體積分?jǐn)?shù)分布對比

選取中心截面、截面A、B 以及流體域（負(fù)Z方向）的氣相體積分?jǐn)?shù)進行分析。截面A、B 位置如下圖8 所示。截面A 緊貼加強筋表面,截面B 距截面A 20 mm。

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

通過圖9 和圖10 可知加強筋頂部開孔與不開孔時均會在流體域的上部聚集氣泡,但聚集的厚度相差較大。當(dāng)換熱達到穩(wěn)定時,在加強筋附近,加強筋頂部開10 mm、20 mm 方孔對應(yīng)的氣泡聚集厚度均不超過5 mm。加強筋不開孔時氣泡將聚集至加強筋第一排兩個圓孔頂部,距離內(nèi)導(dǎo)體垂直距離約30 mm。

圖9 中心截面及流體域(負(fù)Z 方向)氣相體積分?jǐn)?shù)Fig.9 Volume fraction of central section and fluid domain (negative Z direction)

圖10 截面A、B 氣相體積分?jǐn)?shù)Fig.10 Volume fraction of section A and B

通過截面A 處氣相體積分?jǐn)?shù)對比可知,加強筋頂部開孔時,仍會有部分氣相從加強筋的第一、二排圓孔排出,但所占的比例較少。頂部開20 mm 方孔較開10 mm 方孔相比,氣相主要集中在頂部開孔及第一排兩個圓孔排出。不開孔時,氣相先從第一排兩個圓孔排出,隨著氣泡變多,第二排兩個圓孔排出的氣相逐漸增多。對比截面B 處氣相體積分?jǐn)?shù)可知,不開孔時流體下部區(qū)域中氣相占據(jù)的比例較高,上部純氣相占據(jù)的區(qū)域較其他兩種工況顯著增多。過多的氣相會導(dǎo)致傳熱性能變差,抬升內(nèi)導(dǎo)體溫度。

4.2 速度分布對比

中心截面處流場速度分布如圖9 所示。

通過圖11 可知,不開孔時加強筋右側(cè)下部流體速度變化較大,加強筋頂部開10 mm、20 mm 方孔時,加強筋右側(cè)上部流體的速度變化較大。不開孔時,氣泡主要通過加強筋上部4 個圓孔流出,在靠近上部壁面附近的流體速度相對較小。加強筋頂部開孔時,產(chǎn)生的氣泡先在頂部聚集之后從頂部開孔及加強筋第一、二排圓孔流出,所以在流體上部區(qū)域速度較大。當(dāng)氣泡從加強筋排出后,頂部開10 mm 方孔工況對應(yīng)的流體速度分布比較集中,且最大速度均高于其余兩種工況對應(yīng)的最大速度。

圖11 中心截面及流體域(負(fù)Z 方向)速度分布Fig.11 Speed distribution of central section and fluid domain (negative Z direction)

4.3 溫度分布對比

圖12 為中心截面以及內(nèi)導(dǎo)體（負(fù)Z方向）溫度分布。

圖12 中心截面及內(nèi)導(dǎo)體(負(fù)Z 方向)溫度分布Fig.12 Temperature distribution of central section and inner conductor (negative Z direction)

從圖7a 可知,加強筋頂部不開孔時EM 面最高溫度在5 s、30 s 附近出現(xiàn)了較大的變化,5 s 之前,最高溫度幾乎成線性增長,5 s 至30 s 之間,最高溫度依然增長但是增長速率逐步降低,30 s 之后,最高溫度穩(wěn)定在6.25 K 附近,38 s 之后溫度又逐漸升高,48 s之后穩(wěn)定在6.41 K。可能的原因是在5 s 之后氣泡逐漸增多,由于排出不暢,流體加強了對流,溫度增長逐漸變緩。30 s 之后,積聚的氣泡逐漸到達加強筋第一排圓孔的頂部,氣泡逐漸排出,最高溫度不變約8 s 后又逐步升高,最后趨于穩(wěn)定。在此期間平均溫度緩慢升高,最后穩(wěn)定在5.15 K。

從圖7b 可知,加強筋頂部開10 mm 方孔時EM面最高溫度在5 s、10 s 附近出現(xiàn)了較大的變化,在5 s之前,最高溫度幾乎成線性增長,5 s 至7.5 s 之間,最高溫度增長變緩,7.5 s 至10 s 之間增長變快?？赡艿脑蚴? s 之后氣泡逐漸增多,5 s 至7.5 s 之間對流較強,換熱能力提高,但是氣泡沒有有效排出,隨著氣泡的積累導(dǎo)致傳熱變差,溫度增長迅速。10 s 之后,積聚的氣泡通過頂部方孔以及第一排圓孔逐漸排出,最高溫度逐漸趨于穩(wěn)定,最后穩(wěn)定在5.62 K。從平均溫度看,EM 面溫度平均溫度在30 s 之后基本保持在4.93 K,波動很小。

從圖7c 可知,加強筋頂部開20 mm 方孔時EM面最高溫度在5 s、14 s 附近出現(xiàn)了較大的變化,在5 s之前,最高溫度幾乎成線性增長,5 s 至14 s 之間,最高溫度增長變緩,14 s 之后逐漸穩(wěn)定,換熱達到穩(wěn)定時,EM 面最高溫度為5.59 K,平均溫度為4.91 K。與7b 對比可知,開20 mm 方孔時,EM 面最高溫度在前5 s 變化趨勢相同,在5 s 至14 s 之間的變化趨勢相差較大,換熱穩(wěn)定時,二者EM 面最高溫度、平均溫度相差較小,相對誤差分別為0.5%和0.4%。

從圖7a,7b,7c 對比可知,3 種工況下最高溫度開始階段前幾乎線性增長,固體導(dǎo)熱起主導(dǎo)作用,最高溫度大約在5.5 K 后發(fā)生較大變化。加強筋頂部開10 mm、20 mm 方孔工況下,EM 面最高溫度在15 s之后逐步趨于穩(wěn)定,穩(wěn)定在5.62 K。加強筋頂部不開孔時,EM 面最高溫度在30 s 之后逐漸穩(wěn)定,最后穩(wěn)定在6.41 K。

通過圖12 可知3 種工況下內(nèi)導(dǎo)體最高溫度的位置均位于加強筋右側(cè)上部,內(nèi)導(dǎo)體上部溫度梯度較大,下部溫度分布相對較均勻。加強筋頂部開10 mm、20 mm 方孔對應(yīng)的內(nèi)導(dǎo)體溫度分布比較均勻,溫度梯度較小,上部最大溫差約0.8 K。加強筋頂部不開孔時內(nèi)導(dǎo)體溫度梯度較大,上部最高溫差約1.5 K。在3 種工況下,內(nèi)導(dǎo)體最高溫度均處于超導(dǎo)轉(zhuǎn)變溫度9.2 K 以下。

綜上可知,加強筋頂部不開孔時可以保證超導(dǎo)腔溫度處于超導(dǎo)轉(zhuǎn)變溫區(qū)內(nèi),滿足流動傳熱要求。其余兩種工況可以改善流動傳熱效果,但會降低加強筋的強度。為了保證結(jié)構(gòu)強度,最終選擇不額外開孔的方案。

5 結(jié) 論

針對166.6 MHz 超導(dǎo)腔內(nèi)導(dǎo)體附近的流場進行了流動/傳熱/電磁多物理場耦合仿真計算,評估了3種加強筋結(jié)構(gòu)對應(yīng)的內(nèi)導(dǎo)體附近流場的氣相體積分?jǐn)?shù)、速度分布以及溫度分布,結(jié)果表明:

（1）加強筋頂部不開孔會帶來更多的氣泡堆積,抬升內(nèi)導(dǎo)體溫度,EM 面最高溫度為6.41 K,平均溫度為5.15 K。開10 mm、20 mm 方孔時,二者內(nèi)導(dǎo)體EM 面的最高溫度、平均溫度比較接近,分別為5.92 K、4.91 K。與不開孔工況相比最高溫度減少了14%,平均溫度減少了6%;

（2）加強筋頂部不開孔時,內(nèi)導(dǎo)體EM 面溫度梯度較大,開10 mm、20 mm 方孔時內(nèi)導(dǎo)體溫度分布比較均勻,與不開孔時相比最大溫差減少了47%;

（3）加強筋頂部不開孔時,氣相占了較多的體積,但內(nèi)導(dǎo)體溫度上升不太大,可以依靠固體導(dǎo)熱將熱量傳遞到被液氦浸泡的固體部分,內(nèi)導(dǎo)體溫度處于超導(dǎo)溫區(qū),不會發(fā)生失超。因此為了保證結(jié)構(gòu)強度,最終選擇加強筋頂部不額外開孔的方案。