基于高品質因數波導型縮比例船形腔的多物理場耦合仿真研究

劉景源,裴士倫,張天爵,殷治國,汪 洋

(中國原子能科學研究院 核技術綜合研究所,北京 102413)

高能強流質子加速器在核物理與粒子物理等前沿研究領域、大眾健康和先進能源等國民經濟領域,均有著十分廣泛而重要的應用[1-2]。近年來,中國原子能科學研究院提出1套可產生2 GeV、6 MW連續波質子束的強流回旋加速器組合方案,其中包含100 MeV直邊分離扇回旋加速器、800 MeV螺旋分離扇回旋加速器和2 GeV 連續波(CW)固定場交變梯度加速器(FFA),3臺加速器的高頻腔均工作在44.4 MHz,最終可實現加速器組合全流程連續束、等時性高效率加速和高平均質子束流功率輸出[3]。在分圈式等時性加速器中,當束流加速到引出階段,提高螺旋軌道圈的間距,加大留給引出裝置的安裝空間,可有效降低束流被引出裝置阻擋,提高引出效率,實現單圈引出。在回旋式等時性加速器中,在運動方向上,束流的縱向聚焦力較弱,縱向空間電荷效應引起的流強閾值與束流圈能量增益的3次方呈正比[4-5],較橫向空間電荷效應引起的流強閾值低很多。根據PSI實驗室高功率質子加速器的長年運行調試經驗,受空間電荷影響,束流流強與圈能量增益的3次方呈正比,該實驗室590 MeV回旋加速器先后使用的矩形腔體或歐米伽腔體已較難滿足更高的電磁儲能需求[6-8],經過對矩形、歐米伽形和跑道形等多種形狀的波導型高頻腔進行系統研究,提出全流線形的船形高頻腔設計方案[9]。高頻腔徑向變軌范圍約2.8 m,加速電壓可達到2 MV,諧振頻率44.4 MHz、Q值高達90 000。

為掌握船形高頻腔的加工工藝,采用ANSYS HFSS和ANSYS workbench設計工作在177.6 MHz的1∶4縮比例船形腔,開展腔體頻率調諧的仿真研究工作和腔體的樣機研制,同時利用現有的高頻功率源開展相關高頻功率實驗研究。本文將采用穩態熱分析和Fluent兩種熱仿真軟件對該腔體樣機進行多物理場耦合仿真計算,研究船形腔溫度分布、結構變形和頻率變化。

1 設計指標要求

該船形腔經過適當調整,可在不對腔體的高頻性能造成任何影響的情況下,改造為一臺引出電子束能量5～6 MeV、引出靶上平均流強16～20 mA、引出靶上功率90～100 kW的高能電子加速器,圖1所示為加速器基本結構與原理示意圖[10]。

圖1 高能電子加速器基本結構與原理示意圖

綜合考慮2 GeV CW FFA高頻腔的性能指標和1∶4縮比例船形腔樣機在電子回旋加速器工程應用中的要求,縮比例腔需滿足以下要求:1) 縮比例腔體諧振頻率調諧范圍要求能達到±170 kHz左右;2) 高功率運行時,損耗在腔體上的RF熱功率約為100 kW,冷卻水若能將這部分熱量全部帶走,冷卻水管出口相對35 ℃入口水溫的溫升要求<10 ℃;3) 在計入大氣壓力、電動缸拉力或壓力、RF熱功率、重力及各種位移約束等的情況下,機械結構變形導致的等效應力應盡可能低,以方便選用合適形態的無氧銅板材及不銹鋼材。圖2示出縮比例腔機械設計總裝結構,該腔體主要由主腔體、支撐組件及調諧電動缸組件等組成,頻率調諧采用電動缸上下拉壓船形腔體復合板外壁的方法實現。

圖2 縮比例腔機械設計總裝結構

2 多物理場耦合仿真理論

常用的波導型高頻腔有矩形波導型和圓柱形波導型,兩者在電磁場分布和特性參數計算上類似[11]。2 GeV FFA中所使用的高頻腔為矩形波導型,圖3為高頻腔的示意圖,此種類型高頻腔中的模式有兩種,分別為TEmnp模和TMmnp模,m、n及p分別對應x、y及z方向上出現電磁場極大值的個數,a、b及d分別為腔體在x、y及z方向的長度。TEmnp模或TMmnp模的諧振頻率可表示為:

圖3 矩形波導型高頻腔示意圖

(1)

若束流沿z軸穿過高頻腔并獲得加速,腔體需在其運動方向上提供加速電場。考慮到高頻腔中用基模加速,則矩形波導型高頻腔的工作模只能選TM110模,即a>b>d或b>a>d,圖4示出TM110模的電磁場分布形式。為使束流能無阻礙穿過高頻腔,在腔上沿x軸(即圓形加速器的半徑方向)開設了長條形束流通道。腔內電場Ez沿x軸和y軸均呈半正弦分布,在靠近腔體沿x軸兩端的電場Ez太低,不能用于加速,由此束流通道在沿x軸方向上的長度要小于a,既要滿足設計需求又要小于a。由式(1)可知,矩形波導型高頻腔TM110模的頻率與d無關,因此該類型的高頻腔可設計成長窄型(即d較小),這樣非常有利于增加加速器在粒子束運動方向上安裝空間,實際高頻腔與圖4在形狀上會有所差異,需借助3維計算機軟件來對其進行計算和優化[9]。

圖4 矩形波導型高頻腔中TM110模場分布

矩形波導型高頻腔在高功率運行時,在腔體內表面會產生大量熱損耗,這些熱量主要通過熱傳導和熱對流被帶走。根據運行經驗,冷卻水出口溫度相對入口的溫升應≤10 ℃,船形腔冷水流速一般<3 m/s。流速過高容易引起機械震動,流速過低無法及時帶走腔體產生的熱量,經綜合考慮后縮比例腔采用2 m/s[12-13]。整個腔體本體外部布置TU1水管,要求冷卻水能帶走的熱耗散功率約為100 kW,計算分析按照平均35 ℃水的對流換熱系數來計算。選擇10 mm×10 mm、壁厚1 mm的TU1水冷管進行計算。

傳熱仿真計算中,首先需要確定水冷管的流動狀態,流動狀態根據式(2)、(3)計算得到Re。

(2)

(3)

式中:S為非圓斷面的過流截面面積;L為過流截面上流體與水冷管壁接觸的周長;v為介質流速;ρ為介質密度;η為動力黏度系數。物性參數列于表1,可計算出Re=19 875,水冷管處于湍流狀態。由于104

表1 35 ℃水的物性參數

(4)

式中:Prw為普朗特數;κw為冷卻介質的導熱率。式(4)的使用條件是:當流體被加熱時n=0.4,而流體被冷卻時n=0.3;管長與直徑比≥60;壁面溫度和冷卻水主流溫度之差在20～30 ℃[14],對于本文計算腔體,屬于流體被加熱過程,n取0.4。計算得對流換熱系數約為8 100 W/(m2·K)。

通常矩形波導型高頻腔采用鋁或銅等金屬板加工而成,腔體由薄壁金屬板焊接成型[12]。腔體周圍需設計固定結構和支撐結構,防止矩形波導型高頻腔在受到重力、大氣壓力、調諧過程中產生的拉和壓力產生變形。2 GeV FFA中所使用的高頻腔采用Y/2態TU1母材板電子束焊加工。Y/2態母材的抗拉強度、屈服強度和伸長率分別約為220 MPa、103 MPa和65%[15-16];Y/2態電子束焊拉伸試樣的斷裂均位于焊縫處,抗拉強度分別約為235 MPa和225 MPa。其一次和二次許用應力限值分別為68.7 MPa和206 MPa。

3 多物理場耦合仿真

圖5示出了在ANSYS Workbench中開展“RF-熱-結構-RF”全自洽閉環多物理場有限元分析的具體流程,使用Steady-State Thermal模塊開展固熱分析。導入到Geometry中的Solidworks機械模型包含有兩種類型的實體部件:金屬和非金屬實體部件。在HFSS中,與RF場相關的非金屬實體和金屬實體與RF場計算相關,只計算這些區域以使整個模擬計算可順利進行并效率更高。為了使計算分析更接近實際情況,RF計算完成后需調節腔體內部儲能到4 J,由此可使高頻功率損耗達到設計要求的100 kW。在固熱分析中,RF損耗將其設置為熱負載。在這種情況下,將對流系數和壁面冷卻水溫度作為邊界條件,計算得到腔體溫度分布。在靜態結構分析中,將腔體溫度分布結果作為負載導入靜態結構分析模塊,設置重力、大氣壓力、固定支撐等約束,計算得到整個腔體機械模型的變形及應力分布。為了評估結構變形對RF場的影響,選擇對應的實體設置并設置結構變形反饋功能。將計算得到的變形反饋給HFSS,然后再評估機械變形對腔體RF性能的影響。多次迭代上述過程,可獲得穩定的計算結果。但在大多數情況下,一次或兩次迭代即可。

圖5 基于Steady-State Thermal模塊的船形腔多物理場仿真流程圖

固熱分析的方法在設置邊界條件時,對水冷管各處做了溫度分布均勻、對流換熱系數相等的假設。使用Fluent開展流固熱分析時,通過設置運動參數、物性參數和計算方程,仿真過程中逐步計算得到水冷管溫度,這樣的結果更符合實際情況。圖6示出了基于Fluent的多物理場仿真分析流程,使用流體仿真模塊Fluent模塊代替Steady-State Thermal模塊,該方法較固熱分析略微復雜,對模型的要求較高,在不影響整體分析結果的基礎上對導入的模型進行一定的調整,修補不適用于Fluent的模型缺陷。綜合考慮計算機性能和計算時間等條件,在進行網格劃分和流體傳熱仿真時,不考慮船形腔的支撐組件,將HFSS計算得到的腔體內壁功率設置為Fluent的熱負載。設置冷卻水進口流速和溫度,根據式(3)得到Re,求解模型采用k-ε模型。

圖6 基于Fluent的船形腔多物理場仿真流程圖

4 計算結果及討論

表2列出HFSS仿真計算的腔體功率損耗分布,在諧振頻率177.848 MHz下腔體無載Q值約為44 157.7。

表2 功率損耗分布

4.1 溫度分布特性

圖7分別示出穩態熱分析和Fluent仿真的腔體溫度分布,穩態熱分析顯示腔體最高溫度為77 ℃,Fluent仿真計算得到腔體最高溫度為82.6 ℃。兩種方法均顯示腔體最高溫度分布出現在鼻錐兩端根部區域。對比可發現,穩態熱分析結果中腔體溫度沿y-z平面對稱分布,而Fluent仿真計算結果顯示,腔體出口附近溫度比進口附近高約4.7 ℃。腔體熱量可被冷卻水帶走,且冷卻水進出口溫升小于10 ℃,滿足設計指標要求。

a——穩態熱仿真;b——Fluent仿真

在Fluent計算結果中,腔體總流量為Q=6.713 L/s,進出口平均溫升為ΔT=3.6 ℃,計算得到腔體總功率為100.89 kW。在Fluent的計算結果中提取水冷管壁面溫度Tw、熱流密度q及管內主流溫度Tf,水冷管壁面溫度根據牛頓冷卻公式(式(5))計算得到腔體管道對流換熱系數為7 532～8 573 W/(m2·K),與式(2)～(4)計算的理論值基本一致。

q=h(Tw-Tf)

(5)

4.2 腔體變形及頻率調節

為了獲得±170 kHz的頻率調諧范圍,分別在船形腔上部和下部各設置3個電動缸。利用電動缸對船形腔體銅壁施加拉力或壓力后,會使其內部真空部分的體積發生變化。因為拉力或壓力均施加在船形腔內RF磁場較強、電場較弱的區域附近,所以拉腔體會使其頻率降低、壓腔體會使其頻率升高。圖8、9分別示出了利用調諧裝置將船形腔頻率調變-164 kHz和+176 kHz時得到的機械變形及應力分布,最高應力出現在電動缸復合板位置。當腔體頻率調變-164 kHz時,主腔體最大變形為2.64 mm,最高應力為187.22 MPa(二次應力),此時船形腔上部和下部的電動缸對應的總拉力均為30 000 N。當船形腔頻率調變+176 kHz時,主腔體最大變形為2.96 mm,最高應力為187.16 MPa(二次應力),此時船形腔上部和下部電動缸對應的總壓力均為14 000 N。最高應力低于TU1二次應力3[P]=R態204 MPa。

a——機械變形分布;b——應力分布

a——機械變形分布;b——應力分布

在將船形腔諧振頻率調高或調低幾乎相同的量時,總壓力僅約為總拉力的1/2,這是因為船形腔正常工作時內部為真空狀態,主腔體銅壁外表面受到1個大氣壓的壓力作用。無論是利用電動缸拉還是壓腔體,所造成的變形都會在拉力或壓力施加點附近造成較高的應力,實際加工腔體時的焊縫應盡最大可能遠離這些區域。

圖10、11分別示出將電動缸設置為拉力30 000 N和壓力14 000 N,基于Fluent的多物理場有限元仿真得到的溫度分布腔體機械變形及應力分布,最大變形和最高應力都分布在上下電動缸復合板區域。電動缸總拉力為30 000 N時,船形腔主腔體最大形變為2.61 mm,最高應力為187.31 MPa。電動缸總壓力為14 000 N時,船形腔主腔體最大形變為2.50 mm,最高應力為187.18 MPa,腔體調諧范圍為329 kHz。

a——機械變形分布;b——應力分布

對比圖8和圖11可得:相同參數下,兩種仿真方式得到的腔體頻率變化范圍相差11 kHz;基于Fluent的多物理場仿真計算得到的變形小于基于穩態熱分析的結果;前者計算得到的腔體受到的最大應力大于后者的結果。表3、4分別給出最大變形和最高應力。造成這些差異的原因除了計算誤差和溫度分布條件不一致外,還有靜態結構分析約束條件的差異。受計算機性能限制,基于Fluent的仿真流程需要減少對支撐組件的靜力學計算,只能以腔體固定板端面作為固定面。而約束條件的改變,造成原本可以被腔體、固定板和支撐組件共同分攤的應力,只能由腔體來承擔。

表3 腔體最大變形

表4 腔體最高應力

a——機械變形分布;b——應力分布

5 縮比例船形腔冷測

圖12示出腔體在15 ℃的環境溫度下,通過分子泵將腔體真空度提高到8.1×10-8mbar,電機驅動位置從-2.5 mm變化到2.5 mm測得的腔體頻率變化曲線。腔體頻率變化范圍為177.561～177.935 MHz,腔體加工滿足設計需求。在大氣狀態下,利用2個耦合度極小的小耦合環,通過S21掃頻曲線對無載品質因數進行測量,測得Qs≈42 314,約為模擬計算值44 157.7的95.8%。

圖12 腔體冷測頻率曲線

6 結論

本文利用ANSYS workbench開展了基于穩態和Fluent熱分析的船形腔多物理耦合仿真計算。結果表明,腔體最高溫度分布在鼻錐兩端根部區域,主腔體表面溫度分布均勻,水冷管布置合理,可帶走腔體耗散能量。兩種熱分析結果對比顯示,Fluent仿真結果更接近實際情況。靜態結構分析表明應力主要集中在上下電動缸復合板區域,最高應力低于TU1二次應力的許用要求。基于穩態熱分析的多物理耦合仿真和實際冷測均顯示腔體頻率變化范圍滿足設計需求。受計算機算力影響,基于Fluent的多物理耦合仿真對模型優化,無法完全模擬腔體的受力約束情況,導致與穩態熱分析相比,腔體頻率變化范圍的結果相差11 kHz。