基于臺階墊補的磁場優化算法研究

蒲雪江，任志強，萬鑫淼，林鵬太，楊玙菲，李智慧

(四川大學 原子核科學技術研究所，四川 成都 610065)

設計回旋加速器主磁鐵時，通常要留有一定的磁極裕量。在完成磁鐵加工組裝后，基于磁場的實際測量結果，通過磁鐵磁極修正來獲得理想的目標磁場和得到相應磁鐵形狀，這樣可以保證工程的穩妥實施[1]。實際磁場和目標磁場的偏差會直接影響回旋加速器的性能，故選擇合適的磁場優化(得到磁場偏差盡可能小)方法非常重要。

長期以來，國內外研究者用了很多種磁場優化方法，如活極頭削切[2-3]、磁極邊緣墊補[4-7]、磁極面優化[8-10]以及磁極張角優化[11-18]等。本文主要是以圓形軸對稱磁鐵模型為研究對象，開展基于臺階的磁場優化研究。通常而言，按照一定方法來改變磁鐵的幾何形狀就可得到目標磁場。故本文基于磁極形狀和對應磁場幅值響應之間的分析，采用調整磁極面臺階高度的方式[9]優化磁場，求得目標磁場對應的磁極臺階高度。

1 單個墊補塊對磁場的貢獻

1.1 平均場定義

本文研究的圓形軸對稱磁鐵模型如圖1所示。利用CST(EMS模塊，靜磁求解器)[19]軟件計算三維磁場后，考慮到求解精度盡可能高，選擇3次樣條插值法對中心平面上的網格點進行插值計算，然后把中心平面上相同徑向位置處圓周上插值點的幅值軸向分量求平均，作為該位置處的磁場平均值，最終得到沿徑向變化的磁場平均值曲線，稱作平均場Bz,av(r)，用數學公式表達[20]如下：

(1)

其中：Bz(r，θp)為中心平面上在半徑r、角度為θp處的磁場幅值軸向分量；N為圓周上的插值點個數。

圖1 圓形軸對稱磁鐵模型Fig.1 Circular axisymmetrical magnet model

1.2 單個墊補塊對磁場的影響

設磁鐵的整個磁極面半徑為rmax，沿徑向可均勻分為m個臺階，每個臺階的寬度為Δr=rmax/m，第i個磁極臺階的中心位置為ri=(m-0.5)Δr。在徑向位置ri處，墊補高度為Δh的磁極臺階塊，可以得到磁場響應函數為σi(r)。如在徑向200～250 mm處的臺階增加0.4 mm的高度，就可得到如圖2所示的磁場平均場響應曲線和臺階墊補塊。

圖2 高度為0.4 mm的磁極臺階墊補塊和 磁場平均場響應曲線Fig.2 Pole shimming step block with height of 0.4 mm and response curve of magnetic field average field

由圖2可看出，磁場響應曲線近似呈雙高斯分布，且越靠近臺階中心幅值越大，反之位置越遠，幅值呈現迅速減小趨勢。因為高斯分布函數與曲線幅值、寬度和中心位置均相關，考慮到在優化磁場過程中調節磁極臺階高度時，通常會盡量把臺階寬度均勻化，因此本文分析影響σi(r)的兩個因素ri和Δh。

第1步分析徑向位置對磁場響應函數的影響。由前面分析可知，第i個臺階墊補塊在r處的磁場響應函數σi(r)[1]可表示為：

σi(r)=Aifi(r)

(2)

其中：Ai為第i個磁極臺階處的單位高度墊補塊對磁場的響應函數幅值；fi(r)為磁場響應的形狀因子函數，結合磁極臺階幾何特征和圖2墊補塊的磁場曲線特點，可用一歸一化雙高斯函數來擬合形狀因子。

由于磁極臺階是軸對稱，單個臺階墊補塊的磁場響應函數可看成左右相同兩半邊墊補塊磁場響應的疊加。另外，對磁場模擬計算后，由圖3可發現：對于半邊形狀因子擬合曲線，或雙高斯函數擬合后的曲線，或此時模擬計算得到的磁場響應函數曲線，在沿遠離函數中心位置方向幅值呈減小趨勢，其中穿越f(r)=0后反向增大的這一現象可以稱作反向疊加現象。

圖3 實際形狀因子、 擬合形狀因子與拆分的形狀因子比較Fig.3 Comparison of actual shape factor, fitted shape factor and split shape factor

由前面的已知條件可計算得到沿徑向的各單個臺階墊補塊的磁場響應曲線，然后擬合響應函數的幅值點，可得到任意位置處的單個墊補塊的磁場響應函數的幅值，如圖4所示。

圖4 不同徑向位置處墊補單位高度 (0.4 mm)臺階塊的磁場響應幅值Fig.4 Magnetic field response amplitude of shimming step block per unit height (0.4 mm) at different radial positions

如圖4可知，由于在徑向不同位置上墊補相同寬度Δr、相同高度Δh時，墊補塊的磁場響應函數幅值是在中心局部較小區域內增加，隨后持續減小。結合半邊磁極墊補塊和一整個墊補塊的響應幅值擬合曲線可發現：越靠近磁極中心，正向疊加效應越明顯，磁場漏磁減少，磁場幅值變化增大；越靠近磁極外邊緣，磁場反向疊加效應增大，漏磁也明顯增大，幅值變化減小。中間局部區域的貢獻幅值存在抖動現象，其原因是存在計算誤差。而且當各徑向位置處臺階高度相等時，隨半徑增大磁場減小的幅度近似呈現指數增長趨勢[21]。

第2步分析增加墊補塊的高度對磁場響應的影響。結合前一步分析的結果，采取在同一徑向位置處線性增加墊補塊高度，就可得到增加的高度和磁場響應幅值之間的關系；同理，對徑向位置因素分析后，可得到不同位置、不同高度的墊補塊對磁場的響應影響。現以臺階基準高度為95 mm、m=10、不同位置的臺階上線性增加墊補塊的高度為例分析，得到的結果如圖5所示。

由圖5可知，同一臺階處臺階墊補塊的高度在一定范圍內增加，其磁場響應函數幅值近似呈現線性增加的趨勢。

綜合上述分析，對于模擬求解得到的單個墊補塊的磁場響應函數σi(r)，可用形狀因子函數、磁場響應幅值、臺階位置3個變量，雙高斯函數擬合的方式得到的擬合函數來近似代替。

圖5 不同徑向位置、不同臺階墊補塊高度時 磁場響應幅值的變化趨勢Fig.5 Variation trend of magnetic field response amplitude at different radial positions and heights of step shimming block

2 磁場優化算法原理

本文算法主要思路是以初始平均場和目標平均場之間的差值，和單個墊補塊的磁場貢獻初值為基礎，采用調整磁極面上的徑向不同位置處的臺階高度的方法，迭代修正單個墊補塊的磁場貢獻，來減小計算平均場與目標平均場之間的偏差，最終求出目標平均場的近似磁極臺階高度。

根據磁場的對稱性，涉及到的平均場和臺階貢獻的計算區域為中心平面上沿徑向0～482.6 mm。優化算法涉及到的參數列于表1。

表1 優化算法基本參數Table 1 Basic parameter of optimization algorithm

給出1組初始平均場，可得到和目標平均場的偏差ΔB0=Bobj-B0。經過前文的分析，也可以得到單個臺階增加Δh時的磁場貢獻初值為σ0(r)。由于在迭代計算過程中墊補塊的磁場響應函數是不斷修正的，因此不需要非常精準地預測下一次臺階高度的修正量。為簡化計算，在臺階附近及其內部區域的響應幅值是主要的，其余范圍是次要的，可用一系列分段函數來近似替代σ0(r)，即若要調整第i個臺階的高度，則只保留其臺階寬度范圍內的磁場響應幅值，其余范圍均為0，如式(3)所示。

σ0,i(r)=

(3)

再根據墊補公式(式(4))可求得臺階高度修正系數k1：

(4)

在計算過程中，此方程組的解為最小二乘解，用矩陣方式表示為：

然后根據式(5)可求得新的臺階高度h1。

h1=k1Δh+h0

(5)

接下來利用CST軟件計算更新后的三維模型、MATLAB[22]分析磁場數據，可得到新的平均場B1，并判斷臺階墊補后的磁場精度：

max|Bobj-B1|≤pr

(6)

其中，pr為設定的精度。若滿足判據，則完成迭代計算，整個磁場優化過程結束，此時求得的磁極臺階高度可作為目標磁極臺階高度。否則需要進行迭代計算，其過程如下。

結合第t-1代已經得到的ht-1、Bt-1和σt-1(r)數據，由式(4)可得到kt，再利用式(5)修正臺階高度后，將更新的臺階高度代入CST軟件進行計算、分析后可求得Bt，然后利用判據(式(6))進行磁場精度判斷，若滿足判據則完成迭代求解。反之，需對臺階貢獻進行修正，表達式為：

σt(r)=σt-1(r)-ΔBt(r)/kt

(7)

根據式(7)求得更新的σ(r)后，代入式(4)繼續計算，直至迭代計算到設定的最大次數為止，即得最優解。

本文算法中的三維磁場計算是利用CST軟件完成的，利用MATLAB[22]軟件編寫CST-MATLAB磁場優化算法程序，利用MATLAB與CST直接通信，然后計算出平均場結果，并給出臺階修正量，然后通過MATLAB修正臺階高度，進行自動迭代優化。本算法具有如下特點：1) 通用性，只需要掌握圓形軸對稱磁鐵的磁場響應一般規律，就可以根據理想目標磁場來求解對應的磁極形狀，為后續的磁鐵設計、磁場優化方法提供了參考；2) 便捷性，該算法程序可以自動優化計算，無需人工參與，就可以求解出磁鐵的最優形狀參數。

3 磁場優化實例

為驗證算法的有效性，利用該算法對1臺常梯度圓形軸對稱磁鐵的磁極進行優化。磁鐵的磁極半徑為482.6 mm，在0～129.6 mm內，磁場為1.6 T平滑地降低到1.5 T，在129.6～450 mm范圍內，要求磁場降落指數為0.2，設定最大迭代次數為60代。根據磁場降落指數n1[21]的定義為：

(8)

可以得到磁極半徑為129.6～450 mm內的磁場應滿足：

B=B0×(r/r0)-n1

(9)

代入r0和n1可得到隨半徑變化的平均場函數為：

B=3.971×r-0.2

(10)

初始磁鐵極面為平頭，臺階高度為120 mm，為了優化的磁場精度盡可能高，需要適當調整參數m、n。設在實際整個磁極面的中心平面上，沿徑向采樣n=50個點，m=30個臺階，求解對象為滿足目標平均場時對應的磁極臺階高度。計算結果如圖6所示。由圖6可知，墊補場和目標場曲線在0～450 mm范圍內基本重合，體現出了優化算法的有效性。整個迭代收斂過程中，收斂速度先快后慢，在20多代后出現振蕩，說明該算法迭代已經到達最優性能。在經過57次迭代計算后，在磁極半徑為0～408 mm范圍內達最大偏差為0.001 803 T。由于在實際加工過程中會存在精度限制，通常加工精度為0.1 mm，故對臺階高度保留1位小數后，可得到與實際情況更貼近的墊補場曲線。然后把墊補場曲線減去目標場曲線，就可得到經過算法墊補的平均場偏差曲線，如圖7所示。

圖6 初始場、墊補場和目標場平均場的比較Fig.6 Comparison of initial field, shimming field and target field average field

圖7 墊補場與目標場的偏差曲線Fig.7 Deviation curve of shimming field and target field

由圖7可看出，在磁極半徑為24～408 mm范圍內，磁場精度為0.002 2 T。最大偏差發生在第10個臺階158.4 mm處，為0.002 2 T。在磁極半徑小于24 mm或大于408 mm時，偏差有一小范圍的上升，說明磁場還有優化空間，算法還可改進優化。在靠近磁極邊緣的60 mm范圍內，平均場偏差迅速降低后反向增大，其主要原因是邊界的漏磁現象明顯，導致在墊補過程中，磁場貢獻不能很理想地疊加。圖8為墊補后的磁極面臺階圖。整個自動化迭代求解過程體現了算法的便捷性。

圖8 墊補后的磁極臺階高度Fig.8 Pole step height after shimming

4 結論

本文以圓形軸對稱磁鐵為模型，磁極面臺階高度為研究對象，研究了不同墊補塊對磁場變化的影響，推導出以墊補磁極面臺階高度方式的磁場優化算法，設定好參數條件、運行條件后，利用算法程序在1個實例上進行了模擬，最終得到了精度為0.002 2 T的磁場和對應的臺階高度。根據本文算法的特點可得到如下結論：1) 優化的對象為磁極面或近似為磁極面時，可以參考本文的方法；2) 對1臺回旋加速器的磁鐵，磁極臺階貢獻只需要進行1次分析計算，當需要不一樣的目標平均場時，其計算數據可參考使用，這樣減少了大量重復計算，從而提高了計算效率；3) 整個優化計算過程是程序自動化求解，無需人力參與，節約了計算時間，提高了效率。

本文優化計算過程驗證了該算法的可靠性與便捷性，為今后的磁場優化方法提供了參考。