常正則，李 梅，韓瑞雄，孫良瑞，張祥鎮(zhèn)，姜永誠，桑民敬，葉 瑞，李少鵬，葛 銳

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

超導(dǎo)腔是現(xiàn)代超導(dǎo)加速器的核心先進(jìn)部件，可在單位長(zhǎng)度上為束流提供更高的加速電壓和高頻功率，同時(shí)可節(jié)省設(shè)備占用的空間。超導(dǎo)腔多采用純鈮材料制成，使用時(shí)需將其安裝在由液氮和液氦兩層冷屏逐級(jí)熱隔斷的低溫恒溫器中，超導(dǎo)腔本體通常采用冷氦氣和液氦沖刷浸泡的方式進(jìn)行降溫冷卻。超導(dǎo)腔的降溫過程具有以下特點(diǎn)。1) 超導(dǎo)腔是薄壁空腔結(jié)構(gòu)，較易產(chǎn)生形變，無法承受過大的熱應(yīng)力和壓力。2) 在不同的溫度區(qū)間，超導(dǎo)腔的降溫速率要求不同。120 K以上時(shí)要求降溫平緩，表面溫差小，避免熱應(yīng)力過大；120～40 K時(shí)要求降溫速率加快，避免“氫中毒”現(xiàn)象產(chǎn)生[1]；40 K以下則可直接積液氦進(jìn)行冷卻。3) 超導(dǎo)腔形狀復(fù)雜，有部分區(qū)域傳感器無法測(cè)量，且傳感器反應(yīng)時(shí)間長(zhǎng)，難以實(shí)現(xiàn)實(shí)時(shí)反饋。

國內(nèi)外對(duì)超導(dǎo)腔降溫過程的研究較少，王國平等[2]對(duì)BEPCⅡ超導(dǎo)腔垂直測(cè)試杜瓦的漏熱進(jìn)行了分析與實(shí)驗(yàn)研究，得到準(zhǔn)確的漏熱量和漏熱區(qū)域。劉以勇[3]對(duì)超導(dǎo)波蕩器進(jìn)行了冷卻測(cè)試，得到了各部分溫度與時(shí)間的對(duì)應(yīng)曲線。Gorbounov等[4]發(fā)展了一種結(jié)合數(shù)值和分析方法的Fortran計(jì)算程序，可用于計(jì)算單相氦冷卻超導(dǎo)電纜的降溫過程。Chang等[5]對(duì)超導(dǎo)腔運(yùn)行過程中的熱負(fù)荷進(jìn)行了計(jì)算，并與實(shí)驗(yàn)得到的液氦蒸發(fā)量進(jìn)行了對(duì)比。O’Rourke等[6]通過在低熱負(fù)荷真空管內(nèi)預(yù)冷卻低溫恒溫管，可實(shí)現(xiàn)AIST超導(dǎo)加速器(SCA)低溫恒溫器在從室溫(300 K)到4 K降溫區(qū)間內(nèi)液氦損失的蒸發(fā)量降為0。Wenskat等[7]對(duì)歐洲XFEL的AMTF垂直測(cè)試站的降溫過程進(jìn)行了研究，通過對(duì)貼在恒溫器外壁面的3種不同TVO傳感器和兩個(gè)安裝在頂部和底部的Cernox TM傳感器進(jìn)行監(jiān)測(cè)，降溫過程基于PLC可實(shí)現(xiàn)自動(dòng)控制。Dhakal等[8]對(duì)Nb/Cu復(fù)合超導(dǎo)射頻腔從室溫冷卻到1.6 K期間，超導(dǎo)腔的剩余磁場(chǎng)所受的影響進(jìn)行了研究。Huang等[9]對(duì)大晶粒Nb超導(dǎo)射頻腔的表面電阻在冷卻過程中受溫度梯度的影響進(jìn)行了研究。

先進(jìn)光源技術(shù)研發(fā)與測(cè)試平臺(tái)(PAPS)中的650 MHz雙cell超導(dǎo)腔采用波導(dǎo)型高階模耦合器設(shè)計(jì)，可應(yīng)用于高能環(huán)形正負(fù)電子對(duì)撞機(jī)(CEPC)、國際直線對(duì)撞機(jī)(ILC)和上海硬X射線自由電子激光等大型超導(dǎo)腔加速器項(xiàng)目中，因此對(duì)650 MHz雙cell超導(dǎo)腔的自動(dòng)降溫進(jìn)行研究具有重要意義。以上對(duì)超導(dǎo)腔的降溫研究中，缺乏對(duì)其傳熱機(jī)理的分析研究，難以實(shí)現(xiàn)降溫過程的自動(dòng)化。因此，本文以PAPS中的650 MHz雙cell超導(dǎo)腔為例，開展傳熱機(jī)理分析和數(shù)值仿真研究，建立完善的計(jì)算方法和計(jì)算程序，實(shí)現(xiàn)超導(dǎo)腔降溫過程的自動(dòng)化。

1 超導(dǎo)腔自動(dòng)降溫過程機(jī)理

1.1 超導(dǎo)腔結(jié)構(gòu)模型

650 MHz雙cell超導(dǎo)腔結(jié)構(gòu)示意圖如圖1所示，其冷卻結(jié)構(gòu)主要由葫蘆狀的超導(dǎo)腔本身及外部包裹超導(dǎo)腔的外圓筒氦槽構(gòu)成，超導(dǎo)腔的材質(zhì)為純鈮(Nb)，氦槽的材質(zhì)為工業(yè)純鈦(Ti-TA2)。

圖1 650 MHz雙cell超導(dǎo)腔結(jié)構(gòu)示意圖Fig.1 Scheme of 650 MHz double-cell superconducting cavity

650 MHz雙cell超導(dǎo)腔降溫過程的研究對(duì)象主要集中在超導(dǎo)腔、氦槽及其上游的輸運(yùn)管道上，忽略其他附屬部件與它們之間的導(dǎo)熱，這是由于超導(dǎo)腔內(nèi)和氦槽通常處于真空度1×10-5Pa量級(jí)，基本可忽略空氣的導(dǎo)熱和對(duì)流熱流。氦槽外部有液氮和液氦兩層冷屏隔絕，其他部件與超導(dǎo)腔之間的接觸面積較小，因此也可忽略接觸面帶來的導(dǎo)熱熱流。

1.2 自動(dòng)降溫的實(shí)現(xiàn)步驟

超導(dǎo)腔降溫過程的基本思路為分段降溫，即隨降溫時(shí)間的推移，逐步調(diào)整入口參數(shù)，避免腔體溫差過大及降溫速度過慢，最終實(shí)現(xiàn)整體平穩(wěn)、快速、均勻的降溫曲線。設(shè)每段的降溫時(shí)間為1 h，即在該降溫段內(nèi)，入口參數(shù)保持不變，超導(dǎo)腔允許的最大溫差為ΔTmax(K)，降溫速率為Trate(K/h)，初始超導(dǎo)腔和氦槽壁面溫度為Tinitial(K)，降溫目標(biāo)溫度為Ttarget(K)，則總降溫段數(shù)n為：

n=ceiling((Tinitial-Ttarget-ΔTmax)/Trate)+1

(1)

式中：ceiling為向上取整函數(shù)。

第1段冷氦氣的入口溫度為Tinitial-ΔTmax，可確保壁面上不會(huì)產(chǎn)生超過ΔTmax的溫差。入口流量根據(jù)傳熱學(xué)的相關(guān)公式進(jìn)行迭代計(jì)算得到恰當(dāng)?shù)膩砹髁髁?，確保在第1個(gè)降溫階段(1 h)后，目標(biāo)壁面溫度恰好降至Tinitial-Trate。此時(shí)調(diào)整入口溫度和入口流量，進(jìn)入下一個(gè)降溫階段。第n段的冷氦氣入口溫度為Tinitial-ΔTmax-nTrate，同樣迭代計(jì)算本段所需的流量，確保1 h后目標(biāo)壁面溫度降至Tinitial-nTrate，如此循環(huán)往復(fù)直至最后1個(gè)降溫階段。由于要將溫度穩(wěn)至目標(biāo)溫度，因此最后1段的入口溫度即為Ttarget，以目標(biāo)壁面溫度降至Ttarget+0.1 K以內(nèi)為整體降溫過程結(jié)束的標(biāo)志。降溫時(shí)間不再是1 h，而是取決于上一個(gè)降溫段結(jié)束時(shí)的目標(biāo)壁面溫度，因此最后1段降溫時(shí)間τn=(Tn-1-Ttarget)/Trate，得到降溫時(shí)間后，同樣使用迭代計(jì)算方法得到最后1段降溫過程所需的流量。

得到各降溫階段的腔前入口參數(shù)后，再根據(jù)上游流程管道的壓降、摻混等公式進(jìn)行反推，得到固定時(shí)刻兌溫閥門的對(duì)應(yīng)開度，最終通過得到的閥門開度時(shí)刻表實(shí)現(xiàn)降溫過程的自動(dòng)化。根據(jù)以上需求，本文開發(fā)了對(duì)應(yīng)的Fortran程序Autocold V1.0，其中氦氣的熱物性通過實(shí)時(shí)調(diào)用NIST[10]的Refprop v9.0代碼實(shí)現(xiàn)，對(duì)于金屬熱物性則通過實(shí)時(shí)調(diào)用Cryocomp V3.0代碼實(shí)現(xiàn)。Autocold V1.0分為腔前參數(shù)計(jì)算和來流流程計(jì)算兩部分。腔前參數(shù)計(jì)算通過輸入初始溫度、目標(biāo)溫度、最大允許溫差、降溫速率等參數(shù)對(duì)超導(dǎo)腔進(jìn)行傳熱計(jì)算，得到合理的腔前流動(dòng)參數(shù)。來流流程計(jì)算通過流程計(jì)算反推各降溫階段的兌溫閥門開度參數(shù)。

1.3 腔前參數(shù)計(jì)算

1) 流場(chǎng)分析

采用CFD方法對(duì)超導(dǎo)腔流場(chǎng)進(jìn)行穩(wěn)態(tài)分析計(jì)算，對(duì)流場(chǎng)進(jìn)行建模，兩側(cè)壁面設(shè)為絕熱無滑移壁面，湍流模型采用k-ε模型[11]。采用ICEM對(duì)流場(chǎng)進(jìn)行非結(jié)構(gòu)化網(wǎng)格劃分，經(jīng)網(wǎng)格無關(guān)解驗(yàn)證后確定網(wǎng)格數(shù)量為160萬個(gè)，求解器為Ansys CFX 14.0[12]。

入口邊界設(shè)為質(zhì)量流量入口，冷氦氣流量為4 g/s，出口壓力為101 325 Pa，計(jì)算不同入口溫度下的截面流速分布，結(jié)果如圖2所示。

由圖2可知，隨入口溫度的降低，冷氦氣流速不斷減小。超導(dǎo)腔表面流速最快的區(qū)域在靠近入口處，表面流速最慢的區(qū)域在超導(dǎo)腔兩側(cè)的凹腔中。流速在入口下半部分分布不均勻，上半部分則較為均勻。氦槽表面流速的分布情況與之類似。

入口溫度：a——260 K；b——220 K；c——180 K；d——140 K；e——100 K；f——60 K；g——40 K；h——20 K圖2 不同入口溫度下的截面流速分布Fig.2 Velocity distribution under condition of different inlet temperatures

2) 非穩(wěn)態(tài)導(dǎo)熱熱流分析

考慮超導(dǎo)腔和氦槽在厚度方向(徑向)的非穩(wěn)態(tài)導(dǎo)熱熱流，厚度方向的畢渥數(shù)Bi為：

Bi=hδ/λ

(2)

式中：h為表面對(duì)流換熱系數(shù)，W/(m2·K)；δ為厚度，m；λ為固體導(dǎo)熱系數(shù)，W/(m·K)。

不同溫度條件下兩種材料(Nb和Ti-TA2)的導(dǎo)熱系數(shù)如圖3所示?？煽闯?，在40 K溫度以上時(shí)，兩者的導(dǎo)熱系數(shù)在20 W/(m·K)以上。因此Bimax=100×0.004/20=0.02<0.1，滿足集總參數(shù)法的使用條件[13]，所以可忽略厚度方向上的溫度變化。

圖3 不同溫度條件下Nb和Ti-TA2的導(dǎo)熱系數(shù)Fig.3 Thermal conductivity factor of Nb and Ti-TA2 under condition of different temperatures

然后考慮軸向和周向的導(dǎo)熱熱流，在最極端的溫差情況下，兩側(cè)凹腔表面流速為0，僅靠導(dǎo)熱降溫，同時(shí)沖擊點(diǎn)處的對(duì)流換熱系數(shù)無限大，沖擊點(diǎn)的超導(dǎo)腔壁溫等于來流冷氦氣的溫度，為定壁溫邊界?？蓪⒀剌S向(周向類似)的導(dǎo)熱問題簡(jiǎn)化為第一類邊界條件下的半無限大平板非穩(wěn)態(tài)導(dǎo)熱問題，其解析解為：

(3)

式中：T(x,τ)為某一指定τ時(shí)刻x處的溫度；Tw為壁面溫度；T0為初始壁溫；α為熱擴(kuò)散率；η為函數(shù)自變量；erf為誤差函數(shù)，為單調(diào)升函數(shù)。

以超導(dǎo)腔允許最大溫差20 K、初始壁溫為300 K、第1輪入口冷氣溫度為280 K計(jì)算，1 h后的腔體最高溫度為287.4 K，下降12.6 K，因此設(shè)定的降溫速率Trate不應(yīng)大于該值。

3) 對(duì)流熱流分析

超導(dǎo)腔降溫過程可視為管內(nèi)強(qiáng)迫對(duì)流過程，對(duì)于管內(nèi)湍流，可使用Gnielinski公式計(jì)算平均努賽爾數(shù)Nuf[14]：

(1+(d/l)2/3)ct

f=(1.82lgRe-1.64)-2

(4)

式中：下標(biāo)f表示流體；Re為雷諾數(shù)；Pr為普朗特?cái)?shù)；l為管長(zhǎng)，m；d為管內(nèi)徑，m；ct為熱流方向上的溫度校正項(xiàng)；f為管內(nèi)湍流流動(dòng)的Darcy阻力系數(shù)，按照Filonenko公式[14]計(jì)算。

對(duì)氣體有：

ct=(Tf/Tw)0.45Tf/Tw=0.5～1.5

(5)

式中，Tf為流體溫度，K。

對(duì)于管內(nèi)層流，采用Sieder-Tate公式計(jì)算平均Nuf：

(6)

確定管內(nèi)流動(dòng)的平均Nuf后，則平均表面對(duì)流換熱系數(shù)h為：

h=Nu×λ/d

(7)

分別計(jì)算氦槽側(cè)壁面和超導(dǎo)腔側(cè)壁面的平均對(duì)流換熱系數(shù)，然后迭代求解得到冷氦氣的出口溫度。

(Tout-Tin)cpm=hqSq(Tq-Tf)+

hwSw(Tw-Tf)

(8)

式中：Tout和Tin分別為出、入口氦氣溫度；cp為比定壓熱容；m為質(zhì)量流量；hq為超導(dǎo)腔側(cè)的對(duì)流換熱系數(shù)；Sq為超導(dǎo)腔側(cè)的換熱面積；Tq為超導(dǎo)腔溫度；hw為氦槽側(cè)壁面溫度；Sw為氦槽側(cè)換熱面積。

腔前參數(shù)計(jì)算流程如圖4所示。根據(jù)計(jì)算流程，可計(jì)算得到每段降溫過程所需的腔前參數(shù)，為下一步來流流程計(jì)算提供了邊界條件。

1.4 來流流程計(jì)算

圖4 腔前參數(shù)計(jì)算流程Fig.4 Process diagram of cavity parameter

650 MHz雙cell超導(dǎo)腔上游降溫管路的2 500 W@4.5 K氦制冷機(jī)中，有300 K和40 K兩個(gè)兌溫閥門V1和V2，兌溫后經(jīng)一系列管路和彎頭引至主分配閥箱，到達(dá)超導(dǎo)腔及氦槽位置。采用文獻(xiàn)[15]的方法計(jì)算水平圓形直管流動(dòng)摩擦阻力Δp，其中Δp=Δpf+Δpj，Δpf為直管摩擦阻力損失，Δpj為局部摩擦阻力損失。

直管摩擦阻力損失為：

(9)

式中：ρ為流體的密度；u為流體的速度。

局部摩擦阻力損失的計(jì)算采用阻力系數(shù)法，將局部能量損失表示為流體動(dòng)能因子u2/2的1個(gè)倍數(shù)，即：

(10)

式中，ξ為局部阻力系數(shù)。管件與閥門的阻力系數(shù)需由實(shí)驗(yàn)確定，或查閱相關(guān)水力學(xué)手冊(cè)[15]。

管線靜態(tài)漏熱及壓力損失均會(huì)轉(zhuǎn)化為內(nèi)能，引起溫度變化，若出口溫度T2已知，則可求出入口溫度T1，即：

(11)

(12)

式中：Q為熱量，W；m為流體質(zhì)量流量，kg/s。

對(duì)于豎直圓形直管阻力，計(jì)算方法與水平管類似，不同之處在于需考慮流動(dòng)方向，相應(yīng)地增加重力項(xiàng)，本工況中氣流密度較小，且豎直管長(zhǎng)度較短，因此可忽略重力項(xiàng)的影響，按水平管計(jì)算即可。閥門開度與流量對(duì)應(yīng)關(guān)系取決于閥門自身的特性曲線，可查閱相關(guān)閥門的產(chǎn)品手冊(cè)[16]。根據(jù)計(jì)算流程，可反推得到來流摻混所需的兌溫路流量的配比關(guān)系及兌溫閥門開度等參數(shù)。

2 超導(dǎo)腔自動(dòng)降溫的計(jì)算結(jié)果

2.1 腔前參數(shù)計(jì)算結(jié)果

首先計(jì)算650 MHz雙cell超導(dǎo)腔冷卻過程中最常見的一種工況(即工況1)。設(shè)超導(dǎo)腔允許的最大溫差ΔTmax=20 K，設(shè)定的降溫速率Trate=10 K/h，初始超導(dǎo)腔/氦槽壁面溫度Tinitial=300 K，降溫目標(biāo)溫度Ttarget=40 K，目標(biāo)壁面為氦槽壁面?？偨禍貢r(shí)間為26 h，入口冷氦氣的參數(shù)隨降溫段數(shù)的變化如圖5所示?？煽闯?，第1段降溫過程流量略多，這是由于初始超導(dǎo)腔和氦槽的平均溫度相同，之后的降溫過程中，超導(dǎo)腔壁面的平均溫度始終低于氦槽壁面的平均溫度。入口流量在降溫前14 h時(shí)變化不大，之后則大幅下降。這是由于金屬熱容在160 K以下會(huì)大幅下降，而氦氣的熱容卻變化不大，因此此時(shí)降低同等溫度所需的冷氦氣流量會(huì)大幅下降。最后1段降溫過程出現(xiàn)流量增加的現(xiàn)象，這是由于前n-1個(gè)降溫段均是采用過余溫差進(jìn)行降溫，而最后1段需以Ttarget的入口溫度將目標(biāo)壁面降溫至Ttarget，其溫差較小，導(dǎo)致的對(duì)流換熱量也較小，因此需更大的流量。

圖5 工況1下入口冷氦氣參數(shù)隨降溫段數(shù)的變化Fig.5 Inlet cold helium parameter vs cooling section number at case 1

按照不同溫區(qū)所需降溫速率不同，計(jì)算另1組典型降溫工況(即工況2)。在300～120 K的降溫范圍內(nèi)設(shè)定降溫速率Trate=10 K/h，超導(dǎo)腔允許的最大溫差ΔTmax=20 K，為了使超導(dǎo)腔快速通過氫中毒溫區(qū)，減少對(duì)超導(dǎo)腔品質(zhì)的損害，在120 K處穩(wěn)定片刻，隨后在120～40 K降溫范圍內(nèi)進(jìn)行快速降溫，設(shè)定降溫速率Trate=30 K/h，超導(dǎo)腔允許的最大溫差ΔTmax=40 K，目標(biāo)壁面為氦槽壁面。入口冷氦氣的參數(shù)隨降溫段數(shù)的變化如圖6所示(總降溫時(shí)間為20.67 h)。

由圖6可知，前16段的入口流量基本不變，這是由于120 K以上時(shí)，超導(dǎo)腔和氦槽金屬壁面的熱物性變化不大(主要是熱容)，氦氣單位質(zhì)量的比熱容變化也不大，因此每段所需的入口流量較為穩(wěn)定。第17段入口流量大幅增加的原因和上一個(gè)工況中末段流量增加的原因相同。第18段流量較大有兩個(gè)原因：1) 降溫速率加大，單位時(shí)間內(nèi)需吸收更多熱量導(dǎo)致的流量上升；2) 與上一個(gè)工況中第1段流量增加的原因相同。第19段的流量較第18段小，這是由于金屬材料熱容隨溫度的降低而減小導(dǎo)致的，但比前17段高，這是由于降溫速率加大，單位時(shí)間內(nèi)金屬材料所需帶走的熱量較之前略有增加引起的。第20段流量大幅升高的原因則是末段溫差較小。

圖6 工況2下入口冷氦氣參數(shù)隨降溫段數(shù)的變化Fig.6 Inlet cold helium parameter vs cooling section number at case 2

2.2 來流流程計(jì)算結(jié)果

圖7 超導(dǎo)腔入口溫度與流量的配比關(guān)系Fig.7 Relationship between superconducting cavity inlet temperature and flow rate

根據(jù)流動(dòng)條件選擇來流管道尺寸為φ26.9 mm×1.65 mm。根據(jù)流程計(jì)算反推得到超導(dǎo)腔入口溫度與流量的配比關(guān)系，如圖7所示。

若V1和V2兩個(gè)兌溫閥門均采用等百分比的理想閥門，其開度與流量的關(guān)系如圖8所示。

圖8 理想等百分比閥門開度與流量的關(guān)系Fig.8 Opening percent of ideal equal percentage valve vs flow rate

將流程計(jì)算得到的流量配比結(jié)果，根據(jù)閥門開度與流量的對(duì)應(yīng)關(guān)系，即可反推出閥門開度與腔前參數(shù)一一對(duì)應(yīng)的關(guān)系，最終建立閥門開度的時(shí)刻表，從而可實(shí)現(xiàn)降溫過程的自動(dòng)化。V1和V2閥門開度與降溫時(shí)間的關(guān)系如圖9所示。

圖9 V1和V2閥門開度與降溫時(shí)間的關(guān)系Fig.9 Valve opening percent of V1 and V2 vs cooling time

3 結(jié)論

本文以PAPS中的650 MHz雙cell超導(dǎo)腔的降溫過程為主要研究對(duì)象，采用準(zhǔn)則關(guān)系式法及計(jì)算流體力學(xué)方法構(gòu)建了整體換熱模型并編程計(jì)算，得到了降溫過程中超導(dǎo)腔所需的冷氦氣入口腔前參數(shù)。根據(jù)腔前參數(shù)反推流程，考慮整體布局和局部阻力的損失，計(jì)算了兌溫閥門開度與腔前參數(shù)的對(duì)應(yīng)關(guān)系。研究結(jié)果表明，通過正向的流動(dòng)傳熱計(jì)算可實(shí)現(xiàn)降溫過程的自動(dòng)化，因此可用于難以通過反饋控制的場(chǎng)合。在300 K初始溫度、40 K目標(biāo)溫度、降溫速率10 K/h條件下，總降溫時(shí)間為26 h。第1段所需流量最大，為11.8 kg/h，倒數(shù)第2段所需流量最小，為3.7 kg/h。整體降溫速率保持不變，由于金屬熱容隨溫度降低而減小，降溫過程中所需的總流量基本呈下降趨勢(shì)。