中性束注入(Neutral Beam Injection)加熱是磁約束核聚變裝置上主要的輔助加熱手段。強流離子束中性化是中性束注入的關鍵步驟，它直接決定了中性束的注入功率，進而影響對等離子體的加熱效果。本文分析了EAST 中性束注入器束中性化區的氣流特性，發展了基于直接模擬蒙特卡洛(DSMC)方法的氣流模擬程序，并用于EAST 中性束注入器束中性化區的模擬分析和優化設計。

　　中性束注入器(Neutral Beam Injector ，NBI)是將高能中性粒子束注入到聚變裝置內用于加熱等離子體、驅動等離子體電流的裝置，也是EAST 托卡馬克實現其科學目標所必須的輔助加熱系統。EAST-NBI 工作原理如圖1 所示。EAST-NBI 的離子源利用燈絲發射的電子和離子源供氣系統提供的工作氣體，在弧室內通過弧電場的作用形成均勻高密度源等離子體，源等離子體中的離子在離子源束引出系統靜電場的作用下被引出并加速形成高能離子束。與此同時，離子源內未被電離的工作氣體擴散進入中性化室，與中性化室的補充充氣一起形成高能粒子束中性化所需要的電荷轉換靶。在中性化室內，高能離子通過與電荷轉換靶的氣體靶分子碰撞，從而實現束流的中性化。由于離子束的中性化效率不能達到100%，EAST-NBI采用磁偏轉模式利用高能離子通過偏轉磁場時所受到的洛侖茲力的作用，將其偏轉180° 而脫離原來的傳輸方向，進入離子吞食器。經過偏轉系統分離出離子后的高能中性粒子束流經漂移管道注入到EAST 等離子體中。可見，真空技術網(http://smsksx.com/)認為離子束中性化過程對EAST-NBI 的加熱效果十分關鍵。

圖1 EAST-NBI 的示意圖

1、EAST-NBI 束中性化區

　　EAST-NBI 的離子束從電極系統引出加速后即開始中性化過程，因此束中性化區由加速器空腔、插板閥空腔、中性化室以及連接管道構成。根據各個部件的結構特點，將其簡化為如圖2 所示的二維幾何模型。特別地，由于裝配需要，真空室前端蓋板上的束通道開口要大于中性化室的外徑，這使得束中性化區在該連接部分形成與真空室相通的狹縫。束中性化區的機械結構確定以后，氣體靶厚就由進氣量來控制。如圖2 所示，EAST-NBI 束中性化區的氣源來自兩部分，一部分是來自等離子體發生器的前端進氣，另一部分是中性化室的補充氣體。

　　離子束的中性化效率與氣體含量(定義為氣體靶厚)密切相關，中性效率隨著靶厚的增加而增加，但增長率逐漸變小，且最終趨于飽和。在實驗運行時氣體靶厚需要較為準確地控制在最佳靶厚附近。對于EAST-NBI，如果進氣量不足，中性化效率不足比較明顯;相反地，如果進氣量過大的話，不但中性化效率的變化不大，多余的靶氣體還會增加真空系統負荷和真空室氣體壓強，進而導致再電離損失的增大。此外，束與氣體靶碰撞會發生一定的散射，靶厚的增加使得散射次數變多，束因此會更為發散。因此，為了實現束中性化區的合理設計和優化，需要一套研究氣體流動狀態的較為準確的計算方法。

圖2 EAST-NBI 束中性化區內氣體流動狀態

2、束中性化區氣流模擬程序

　　由EAST-NBI 束中性化區內氣流的稀薄程度和系統的特征長度可得，該氣流屬于過渡流區域和分子流區域。為了滿足跨流域模擬的需要，采用直接模擬蒙特卡洛(DSMC)方法發展NBI束中性化區氣流模擬程序。本文參考標準DSMC程序，其流程如圖3 所示，并根據束中性化區的氣流特點進行修改。

2.1、程序初始化

　　準備模擬所需的常數(例如，玻爾茲曼常數)，設定氣體和流場基本屬性，賦予各變量初始值等。其中氣體和流場基本屬性主要包括：坐標系(例如，直角坐標系或柱坐標系)、計算區域和邊界條件(包括固壁邊界、入口邊界、出口邊界)、劃分的網格和子網格、模擬分子基數、模擬時間步長、入口的進氣量、進氣溫度、進氣速度、氘分子半徑、氘分子質量、氣體粘滯- 溫度冪指數、氣體分子碰撞模型、氣體的自由度。

圖3 標準DSMC 程序流程圖

2.2、氣源模擬分子注入

　　根據各入口的進氣量Qin 和進氣溫度Tin，計算單位步長Δt 內進入中性化室的模擬分子個數Nin。然后，根據進氣溫度的Maxwell 速度分布函數，運用MC 方法確定入射模擬分子的各個速度分量的大小。再根據入口邊界條件中對進入分子的空間分布和時間分布假設，利用蒙特卡洛方法確定每個入射模擬分子在入口平面上的入射點坐標xin 和入射時序tin(<Δt)。最終，一般假設分子時間分布均勻，則由勻速直線運動方程求得每個入射模擬分子的位置坐標x=xin+u·tin=xin+u·Δt。

2.3、模擬分子位置推進

　　先由勻速直線運動方程求出每個模擬分子以各自速度在Δt 內運動的距離，從而確定模擬分子新的位置坐標x=xo+u·Δt。模擬分子在經歷遷移運動后有可能與邊界發生碰撞。因此需要將每個模擬分子新的坐標與所有的邊界坐標進行比較，判斷是否與邊界發生碰撞。如果發生碰撞的邊界是真空，則將模擬分子作為逸出處理，不再跟蹤計算。如果發生相互作用的邊界是壁面，則先用模擬分子的直線運動方程和邊界的面方程聯立求出碰撞點xc，再根據碰撞前模擬分子運動距離和速度求出碰撞前的運動耗時tc。然后由對應的分子壁面碰撞模型來確定碰撞后的分子的各個速度分量ur。最終，經過碰撞的模擬分子新的位置坐標為x=xc+ur·(Δt-tc)，實現了模擬分子位置的推進。

2.4、分子排序

　　根據模擬分子新的位置坐標調整模擬分子所屬的網格編號，統計各個網格內的分子總數，并按照網格的順序(在每個網格中又按照子網格的順序)對模擬分子進行排序，最終分子編號將被置于一組交叉引用數組中。這是為了之后選取“碰撞對”做準備，據此計算出每個網格所需要模擬的碰撞次數，并保證只有同一網格內的分子才會發生碰撞。

2.5、分子間碰撞

　　本文采用NTC 方法進行碰撞對的取樣，先計算出Δt 內每個網格中需要模擬的碰撞次數Ncol。然后利用蒙特卡洛方法在分子編號的交叉引數組中，隨機選取同一個子網格(或網格)中的兩個模擬分子作為“碰撞對”。再計算兩個分子的相對速度cr 和碰撞截面σ，本文采用可變徑硬球(VHS)模型，其碰撞截面隨著相對速度增加而減小。接著判斷“碰撞對”是否被選中發碰撞。如果判斷“碰撞對”發生碰撞，則依據VHS 模型的碰撞法則計算分子碰撞后的速度和內能。如果判斷沒有碰撞則繼續選取“碰撞對”，依次完成每個網格碰撞次數。

2.6、流動性質采樣

　　包括采樣總數以及每個網格累加的模擬分子數的總和、各速度分量的總和、速度平方和的總和。每個Δt 的模擬的時間、模擬分子數、碰撞總次數、推進模擬分子次數等參數實時輸出顯示，可以初步判斷模型是否存在錯誤。每經過jnis 個Δt 進行一次流動性質采樣，從而得到速度、溫度、密度等宏觀量隨時間的變化規律。重復nloop個Δt 直到統計誤差滿足要求，特別是非定常流動的模擬，需要獲得足夠多的總計樣本。最后求出截止時刻的各物理量的系綜平均。

2.7、模擬結果的輸出

　　將每個網格累加的模擬分子數的總和、各速度分量的總和、速度平方和的總和分別除以采樣總數，得到每個網格相關變量的統計平均值。最后，由這些變量統計平均值計算得到流場的性質并輸出到文本。流場性質包括每個網格的坐標、氣體密度、平動溫度、旋轉溫度、速度分量。

3、模擬結果分析

3.1、EAST-NBI 束中性化區氣流特性

　　在模擬計算中，工作氣體D2 被分別從弧室和中性化室進入EAST-NBI 的束中性化區，為方便比較，兩個進氣量均設為15 Torr·l/s。計算區域中網格的大小設為Δx=0.02 m、Δy=0.01 m，每個網格又劃分為2×2 個子網格。每個模擬分子代表的真實分子數為1014。設時間步長Δt=5×10-6 s，每5 個Δt 進行一次采樣，總的采樣次數為40000，總的模擬時間為1 s，并且氣流在4000 次采樣后達到穩定。

　　本文僅進行未引出束時的氣流模擬，因此氣體溫度和壁面溫度均設定為室溫。進入分子在時間和空間上分布均勻， 初始速度遵循Maxwell 分布。為了模擬分子間碰撞，故采用VHS模型，而分子與壁面碰撞采用漫反射模型。特別地，設模擬分子與弧室側的邊界碰撞時，77%的模擬分子在碰撞點反射，而23%從邊界其他點反射。圖4 給出了兩種進氣模式計算得到的速度矢量分布。如圖4(a)所示，氘氣穿過整個電極系統的大截面，以層流的狀態緩慢地進入束中性化區。但是氣體進入溝槽區域后隨即迅速向兩側擴散，同時其軸向速度停止增加反而急劇減小。由于溝槽區域的橫截面較大，壁面對氣流的粘滯作用將減弱，引起質量流速率的增加。與此同時，溝槽的廣大空間使得其中的氣體接近熱平衡狀態，這在氣體密度分布中也能體現出來(如圖5(a)所示)。此外在溝槽區域內的小狹縫與計算區域外的真空環境相連接，從而在較高的壓力梯度下將氣體推出。氣體繼續向前運動，進入相對十分收縮的中性化室桶體，導致在桶體入口處現了堵塞。從圖5(a)可以看出，兩個對稱的滯止點出現在中性化室入口附近的壁面上，這是源于亞聲速氣流掠過粘滯平面的結果。

圖4 兩種進氣模式的速度矢量分布

圖5 兩種進氣模式的氣體密度分布

　　而在中性化室進氣模式中，氘氣以高壓縮比從側面的進氣口噴出并在對面的壁面上形成滯止點。不過由于實際的進氣口在上面，進氣口距離對面的壁面較遠，很難形成滯止點。兩個宏觀量的分布圖中最明顯的是，在束中性化區前端的大部分區域內，氣體速度矢量幾乎為零，而且氣體密度大小十分接近。其原因是電極系統側的邊界條件采用了完全漫反射模型，反射氣流和來流相互碰撞而使得氣流趨于各向同性。流場的另一個特征是在進氣口和狹縫之間具有兩個剪切區，把束中性化區內的氣體向外擠壓，特別是在靠近進氣口一側還出現了一個環形渦流。兩種進氣模式有一個相同的流場特征，徑向壓強梯度在束中性化區的中部形成，并推動氣體向出口流動。

　　EAST-NBI 束中性化區氣流特性的定量分析結果如圖6 所示，兩種進氣模式的區別很明顯。Kn 數(0.2 ~ 2.4)表明了束中性化區的氣流屬于過渡流和自由分子流領域，與定性分析相符。在一定體積內，氣體溫度會隨著氣體分子的碰撞頻率增加而上升，所以兩條溫度曲線都是在進氣口附近的區域較高。而在出口附近的區域各個分子趨于同向，碰撞頻率相應地減小，溫度也相應地驟減。考慮到模型在出口邊界采用真空邊界條件，而在實際情況中，出口并不是完全真空甚至還有返流，因此出口附近的溫度會高于模擬結果。氣體壓強曲線可以跟日后的測量值進行比較，在真空邊界條件下模擬得到的出口壓強(< 0.03 Pa)滿足EAST-NBI 的工程要求。

圖6 EAST-NBI 束中性化區內各個中心線宏觀量的分布：(a)Knudsen 數;(b)氣體溫度;(c)氣體壓強;(d)軸向速度

3.2、EAST-NBI 束中性化區的優化設計

　　在實驗運行中，總是希望用最小的進氣量來獲得要求的氣體靶厚，換句話說，就是盡可能地提高氣體利用率，從而降低前文多次提及的剩余氣體帶來的一系列副作用。在實際裝置中，如幾何模型所示的狹縫在短邊的寬度約為5 cm，在長邊的寬度約為2 cm。由于這些狹縫直接與真空室的真空環境相連接，如同束中性化區內氣體的又一個出口，很大程度上減少了其中的氣體密度。為此，在未來EAST-NBI 的設計裝配中，將借鑒國外其他NBI 裝置的設計把狹縫盡可能的縮小，甚至是完全封閉起來。

　　和上述的模擬條件相同，工作氣體D2 被分別從弧室和中性化室進入EAST-NBI 的束中性化區，為方便比較，兩個進氣量均設為15 Torr·l/s。圖7 給出了從狹縫流出氣體所占比例與氣體靶厚隨狹縫寬度的變化情況。在兩種進氣模式中，隨著狹縫一點點增大，中性化區域的靶厚下降得十分迅速。當狹縫增大到4 cm 時，靶厚已經降到了沒有狹縫時的一半，而且由于狹縫距離氣源更近，此時從狹縫流出的氣體已占主導，超過所有逸出氣體分子的50%。根據模擬結果，在以后的裝配中至少要控制狹縫的寬度小于1 cm。

圖7 EAST-NBI 中性化室的安裝縫隙對氣體靶厚的影響

　　在EAST-NBI 綜合測試平臺上進行中性束引出實驗時，對中性化室補氣會引起離子源電極打火次數增多，而且弧室內氣體放電的穩定性變差。中性化室補充進氣的目的是增加靶厚，提高中性化效率。但由中性化室進入的氣體會向離子源方向擴散，并通過加速器并進入弧室中。這些氣體會增加沉積在電極表面的高能粒子數量，甚至影響弧室中的氣體放電過程。為了解決這一問題，需要在氣體靶厚不變(或減少很小)的前提下，盡可能的減少流向電極系統的氣體。由于EAST-NBI 真空室內各部件的尺寸和位置很難改變，改變中性化室進氣的方向是一個較為簡單可行的解決方案。為此，在EAST-NBI 束中性化區模型中模擬了不同進氣角度對氣體靶厚和離子源運行的影響。中性化室進氣角度的定義如圖8 所示。在計算中進氣角度范圍從45°到135°，進氣量均設為15 Torr·l/s。靠近弧室側的氣體數密度用于反映流向弧室的氣流大小，在圖8 中用圓點表示。目前裝置的進氣角度為90°，為了方便比較，以進氣角度為90°時的氣體靶厚和弧室側氣體數密度作歸一化處理。可以預見，隨著進氣角度增加，靠近弧室側的氣體數密度也將變大， 并且在60° 和120°之間的增大趨勢最明顯。圖8 還顯示了不同進氣角度對應的氣體靶厚(用方點表示)。可以看出，不同進氣角度的靶厚變化范圍要小于弧室側氣體數密度的變化范圍。這意味著，流向弧室的氣體量可以通過減小進氣角度實現，而同時氣體靶厚的減小量更小。而且當進氣角度小于60°時，兩個量的減小趨勢都變緩。因此提議將中性化室的進氣角度從90°改為60°，弧室側的氣體數密度和束中性化區氣體靶厚分別減小90%和97%。

圖8 EAST-NBI 中性化室進氣角度對離子源運行和靶厚的影響

4、結論

　　離子束中性化過程直接決定了中性束注入器的性能參數，而氣體靶作為中性化媒介最為重要，考慮到束中性化區的氣流屬于過渡流領域，本文將該流域的標準方法———直接模擬蒙特卡洛(DSMC)方法引入NBI 的相關研究，建立了NBI 束中性化區氣流模擬程序。利用該程序EAST-NBI束中性化區內的氣流特性，計算結果顯示：由于考慮了氣體分子間碰撞，DSMC 方法的結果較符合氣流的實際情況，并且根據氣流的速度、溫度、壓強等宏觀性質的細節分布，分析了氣體靶的形成機制。此外，還將該程序對EAST-NBI 束中性化區的優化設計進行模擬計算，結果顯示：中性化室和真空室前端蓋板的束通道之間的狹縫寬度必須控制小于1 cm，才能有效的提高氣體利用率;將中性化室的進氣角度從90°改為60°，能夠在氣體靶厚不變(或減少很小)的前提下，減少流向電極系統的氣體，增加離子源運行的穩定性。