常正則,李 梅,韓瑞雄,孫良瑞,張祥鎮(zhèn),姜永誠(chéng),桑民敬,葉 瑞,李少鵬,葛 銳
(中國(guó)科學(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)腔的降溫過(guò)程具有以下特點(diǎn)。1) 超導(dǎo)腔是薄壁空腔結(jié)構(gòu),較易產(chǎn)生形變,無(wú)法承受過(guò)大的熱應(yīng)力和壓力。2) 在不同的溫度區(qū)間,超導(dǎo)腔的降溫速率要求不同。120 K以上時(shí)要求降溫平緩,表面溫差小,避免熱應(yīng)力過(guò)大;120~40 K時(shí)要求降溫速率加快,避免“氫中毒”現(xiàn)象產(chǎn)生[1];40 K以下則可直接積液氦進(jìn)行冷卻。3) 超導(dǎo)腔形狀復(fù)雜,有部分區(qū)域傳感器無(wú)法測(cè)量,且傳感器反應(yīng)時(shí)間長(zhǎng),難以實(shí)現(xiàn)實(shí)時(shí)反饋。
國(guó)內(nèi)外對(duì)超導(dǎo)腔降溫過(guò)程的研究較少,王國(guó)平等[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)曲線(xiàn)。Gorbounov等[4]發(fā)展了一種結(jié)合數(shù)值和分析方法的Fortran計(jì)算程序,可用于計(jì)算單相氦冷卻超導(dǎo)電纜的降溫過(guò)程。Chang等[5]對(duì)超導(dǎo)腔運(yùn)行過(guò)程中的熱負(fù)荷進(jìn)行了計(jì)算,并與實(shí)驗(yàn)得到的液氦蒸發(fā)量進(jìn)行了對(duì)比。O’Rourke等[6]通過(guò)在低熱負(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è)試站的降溫過(guò)程進(jìn)行了研究,通過(guò)對(duì)貼在恒溫器外壁面的3種不同TVO傳感器和兩個(gè)安裝在頂部和底部的Cernox TM傳感器進(jìn)行監(jiān)測(cè),降溫過(guò)程基于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)射頻腔的表面電阻在冷卻過(guò)程中受溫度梯度的影響進(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)、國(guó)際直線(xiàn)對(duì)撞機(jī)(ILC)和上海硬X射線(xiàn)自由電子激光等大型超導(dǎo)腔加速器項(xiàng)目中,因此對(duì)650 MHz雙cell超導(dǎo)腔的自動(dòng)降溫進(jìn)行研究具有重要意義。以上對(duì)超導(dǎo)腔的降溫研究中,缺乏對(duì)其傳熱機(jī)理的分析研究,難以實(shí)現(xiàn)降溫過(guò)程的自動(dòng)化。因此,本文以PAPS中的650 MHz雙cell超導(dǎo)腔為例,開(kāi)展傳熱機(jī)理分析和數(shù)值仿真研究,建立完善的計(jì)算方法和計(jì)算程序,實(shí)現(xiàn)超導(dǎo)腔降溫過(guò)程的自動(dòng)化。
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)腔降溫過(guò)程的研究對(duì)象主要集中在超導(dǎo)腔、氦槽及其上游的輸運(yùn)管道上,忽略其他附屬部件與它們之間的導(dǎo)熱,這是由于超導(dǎo)腔內(nèi)和氦槽通常處于真空度1×10-5Pa量級(jí),基本可忽略空氣的導(dǎo)熱和對(duì)流熱流。氦槽外部有液氮和液氦兩層冷屏隔絕,其他部件與超導(dǎo)腔之間的接觸面積較小,因此也可忽略接觸面帶來(lái)的導(dǎo)熱熱流。
超導(dǎo)腔降溫過(guò)程的基本思路為分段降溫,即隨降溫時(shí)間的推移,逐步調(diào)整入口參數(shù),避免腔體溫差過(guò)大及降溫速度過(guò)慢,最終實(shí)現(xiàn)整體平穩(wěn)、快速、均勻的降溫曲線(xiàn)。設(shè)每段的降溫時(shí)間為1 h,即在該降溫段內(nèi),入口參數(shù)保持不變,超導(dǎo)腔允許的最大溫差為ΔTmax(K),降溫速率為T(mén)rate(K/h),初始超導(dǎo)腔和氦槽壁面溫度為T(mén)initial(K),降溫目標(biāo)溫度為T(mén)target(K),則總降溫段數(shù)n為:
n=ceiling((Tinitial-Ttarget-ΔTmax)/Trate)+1
(1)
式中:ceiling為向上取整函數(shù)。
第1段冷氦氣的入口溫度為T(mén)initial-ΔTmax,可確保壁面上不會(huì)產(chǎn)生超過(guò)ΔTmax的溫差。入口流量根據(jù)傳熱學(xué)的相關(guān)公式進(jìn)行迭代計(jì)算得到恰當(dāng)?shù)膩?lái)流流量,確保在第1個(gè)降溫階段(1 h)后,目標(biāo)壁面溫度恰好降至Tinitial-Trate。此時(shí)調(diào)整入口溫度和入口流量,進(jìn)入下一個(gè)降溫階段。第n段的冷氦氣入口溫度為T(mén)initial-ΔTmax-nTrate,同樣迭代計(jì)算本段所需的流量,確保1 h后目標(biāo)壁面溫度降至Tinitial-nTrate,如此循環(huán)往復(fù)直至最后1個(gè)降溫階段。由于要將溫度穩(wěn)至目標(biāo)溫度,因此最后1段的入口溫度即為T(mén)target,以目標(biāo)壁面溫度降至Ttarget+0.1 K以?xún)?nèi)為整體降溫過(guò)程結(jié)束的標(biāo)志。降溫時(shí)間不再是1 h,而是取決于上一個(gè)降溫段結(jié)束時(shí)的目標(biāo)壁面溫度,因此最后1段降溫時(shí)間τn=(Tn-1-Ttarget)/Trate,得到降溫時(shí)間后,同樣使用迭代計(jì)算方法得到最后1段降溫過(guò)程所需的流量。
得到各降溫階段的腔前入口參數(shù)后,再根據(jù)上游流程管道的壓降、摻混等公式進(jìn)行反推,得到固定時(shí)刻兌溫閥門(mén)的對(duì)應(yīng)開(kāi)度,最終通過(guò)得到的閥門(mén)開(kāi)度時(shí)刻表實(shí)現(xiàn)降溫過(guò)程的自動(dòng)化。根據(jù)以上需求,本文開(kāi)發(fā)了對(duì)應(yīng)的Fortran程序Autocold V1.0,其中氦氣的熱物性通過(guò)實(shí)時(shí)調(diào)用NIST[10]的Refprop v9.0代碼實(shí)現(xiàn),對(duì)于金屬熱物性則通過(guò)實(shí)時(shí)調(diào)用Cryocomp V3.0代碼實(shí)現(xiàn)。Autocold V1.0分為腔前參數(shù)計(jì)算和來(lái)流流程計(jì)算兩部分。腔前參數(shù)計(jì)算通過(guò)輸入初始溫度、目標(biāo)溫度、最大允許溫差、降溫速率等參數(shù)對(duì)超導(dǎo)腔進(jìn)行傳熱計(jì)算,得到合理的腔前流動(dòng)參數(shù)。來(lái)流流程計(jì)算通過(guò)流程計(jì)算反推各降溫階段的兌溫閥門(mén)開(kāi)度參數(shù)。
1) 流場(chǎng)分析
采用CFD方法對(duì)超導(dǎo)腔流場(chǎng)進(jìn)行穩(wěn)態(tài)分析計(jì)算,對(duì)流場(chǎng)進(jìn)行建模,兩側(cè)壁面設(shè)為絕熱無(wú)滑移壁面,湍流模型采用k-ε模型[11]。采用ICEM對(duì)流場(chǎng)進(jìn)行非結(jié)構(gòu)化網(wǎng)格劃分,經(jīng)網(wǎng)格無(wú)關(guān)解驗(yàn)證后確定網(wǎng)格數(shù)量為160萬(wàn)個(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è)的凹腔中。流速在入口下半部分分布不均勻,上半部分則較為均勻。氦槽表面流速的分布情況與之類(lèi)似。
入口溫度: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,滿(mǎn)足集總參數(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ù)無(wú)限大,沖擊點(diǎn)的超導(dǎo)腔壁溫等于來(lái)流冷氦氣的溫度,為定壁溫邊界??蓪⒀剌S向(周向類(lèi)似)的導(dǎo)熱問(wèn)題簡(jiǎn)化為第一類(lèi)邊界條件下的半無(wú)限大平板非穩(wěn)態(tài)導(dǎo)熱問(wèn)題,其解析解為:
(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)腔降溫過(guò)程可視為管內(nèi)強(qiáng)迫對(duì)流過(guò)程,對(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ì)算得到每段降溫過(guò)程所需的腔前參數(shù),為下一步來(lái)流流程計(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è)兌溫閥門(mén)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ù)。管件與閥門(mén)的阻力系數(shù)需由實(shí)驗(yàn)確定,或查閱相關(guān)水力學(xué)手冊(cè)[15]。
管線(xiàn)靜態(tài)漏熱及壓力損失均會(huì)轉(zhuǎn)化為內(nèi)能,引起溫度變化,若出口溫度T2已知,則可求出入口溫度T1,即:
(11)
(12)
式中:Q為熱量,W;m為流體質(zhì)量流量,kg/s。
對(duì)于豎直圓形直管阻力,計(jì)算方法與水平管類(lèi)似,不同之處在于需考慮流動(dòng)方向,相應(yīng)地增加重力項(xiàng),本工況中氣流密度較小,且豎直管長(zhǎng)度較短,因此可忽略重力項(xiàng)的影響,按水平管計(jì)算即可。閥門(mén)開(kāi)度與流量對(duì)應(yīng)關(guān)系取決于閥門(mén)自身的特性曲線(xiàn),可查閱相關(guān)閥門(mén)的產(chǎn)品手冊(cè)[16]。根據(jù)計(jì)算流程,可反推得到來(lái)流摻混所需的兌溫路流量的配比關(guān)系及兌溫閥門(mén)開(kāi)度等參數(shù)。
首先計(jì)算650 MHz雙cell超導(dǎo)腔冷卻過(guò)程中最常見(jiàn)的一種工況(即工況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段降溫過(guò)程流量略多,這是由于初始超導(dǎo)腔和氦槽的平均溫度相同,之后的降溫過(guò)程中,超導(dǎo)腔壁面的平均溫度始終低于氦槽壁面的平均溫度。入口流量在降溫前14 h時(shí)變化不大,之后則大幅下降。這是由于金屬熱容在160 K以下會(huì)大幅下降,而氦氣的熱容卻變化不大,因此此時(shí)降低同等溫度所需的冷氦氣流量會(huì)大幅下降。最后1段降溫過(guò)程出現(xiàn)流量增加的現(xiàn)象,這是由于前n-1個(gè)降溫段均是采用過(guò)余溫差進(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)腔快速通過(guò)氫中毒溫區(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
圖7 超導(dǎo)腔入口溫度與流量的配比關(guān)系Fig.7 Relationship between superconducting cavity inlet temperature and flow rate
根據(jù)流動(dòng)條件選擇來(lái)流管道尺寸為φ26.9 mm×1.65 mm。根據(jù)流程計(jì)算反推得到超導(dǎo)腔入口溫度與流量的配比關(guān)系,如圖7所示。
若V1和V2兩個(gè)兌溫閥門(mén)均采用等百分比的理想閥門(mén),其開(kāi)度與流量的關(guān)系如圖8所示。
圖8 理想等百分比閥門(mén)開(kāi)度與流量的關(guān)系Fig.8 Opening percent of ideal equal percentage valve vs flow rate
將流程計(jì)算得到的流量配比結(jié)果,根據(jù)閥門(mén)開(kāi)度與流量的對(duì)應(yīng)關(guān)系,即可反推出閥門(mén)開(kāi)度與腔前參數(shù)一一對(duì)應(yīng)的關(guān)系,最終建立閥門(mén)開(kāi)度的時(shí)刻表,從而可實(shí)現(xiàn)降溫過(guò)程的自動(dòng)化。V1和V2閥門(mén)開(kāi)度與降溫時(shí)間的關(guān)系如圖9所示。
圖9 V1和V2閥門(mén)開(kāi)度與降溫時(shí)間的關(guān)系Fig.9 Valve opening percent of V1 and V2 vs cooling time
本文以PAPS中的650 MHz雙cell超導(dǎo)腔的降溫過(guò)程為主要研究對(duì)象,采用準(zhǔn)則關(guān)系式法及計(jì)算流體力學(xué)方法構(gòu)建了整體換熱模型并編程計(jì)算,得到了降溫過(guò)程中超導(dǎo)腔所需的冷氦氣入口腔前參數(shù)。根據(jù)腔前參數(shù)反推流程,考慮整體布局和局部阻力的損失,計(jì)算了兌溫閥門(mén)開(kāi)度與腔前參數(shù)的對(duì)應(yīng)關(guān)系。研究結(jié)果表明,通過(guò)正向的流動(dòng)傳熱計(jì)算可實(shí)現(xiàn)降溫過(guò)程的自動(dòng)化,因此可用于難以通過(guò)反饋控制的場(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。整體降溫速率保持不變,由于金屬熱容隨溫度降低而減小,降溫過(guò)程中所需的總流量基本呈下降趨勢(shì)。