邊界元法計(jì)算束流位置探測(cè)器的靈敏度系數(shù)

李 勇 何 俊 隨艷峰 袁任賢 杜垚垚 張 婉 王徐建 尹 頔 于令達(dá) 劉 智 魏書(shū)軍 岳軍會(huì) 曹建社

1（中國(guó)科學(xué)院高能物理研究所 北京 100049）

2（中國(guó)科學(xué)院大學(xué) 北京 100049）

3（中國(guó)科學(xué)院上海高等研究院 上海 201204）

高能同步輻射光源（High Energy Photon Source，HEPS）是我國(guó)第一臺(tái)第四代同步輻射光源，目前正處于工程建設(shè)階段［1-2］。其能量為6 GeV，自然發(fā)射度為34 pm·rad，建成后將為各種基礎(chǔ)科學(xué)研究提供能量高達(dá)300 keV的高性能X射線［3-4］。第四代光源突出的特點(diǎn)是亮度高、發(fā)射度小，光源點(diǎn)尺寸比三代光源要小得多，而一般的光學(xué)實(shí)驗(yàn)對(duì)光源點(diǎn)的穩(wěn)定性要求為束流尺寸的10%，相應(yīng)的束流軌道穩(wěn)定性也提到了亞微米級(jí)，束流位置測(cè)量的分辨率與精度要求也提高了。束流位置測(cè)量由束流位置探測(cè)器（Beam Position Monitor，BPM）來(lái)完成，通常由束流位置探頭和電子學(xué)組成［5-6］。BPM的位置靈敏度系數(shù)是探頭信號(hào)與束流實(shí)際位置轉(zhuǎn)化過(guò)程中的關(guān)鍵參數(shù)。在BPM探頭加工完成后，可以使用標(biāo)定系統(tǒng)對(duì)其進(jìn)行標(biāo)定，對(duì)位置靈敏度系數(shù)進(jìn)行精確測(cè)量［7］。在BPM設(shè)計(jì)階段，利用電磁仿真軟件如CST等，也可以模擬標(biāo)定過(guò)程，得到設(shè)計(jì)的理論值［8］。本文則主要介紹一種基于有限元思想的計(jì)算位置靈敏度系數(shù)的方法——邊界元分析方法。

1 理論

1.1 位置靈敏度系數(shù)的計(jì)算

BPM探頭通常由腔體和四個(gè)拾取電極組成，如圖1所示。BPM輸出信號(hào)是束流位置、束流強(qiáng)度、束團(tuán)長(zhǎng)度以及其他因素（如線纜長(zhǎng)度、電子學(xué)工作頻率、環(huán)境溫度）等的函數(shù)，第i個(gè)拾取電極上的輸出信號(hào)在時(shí)域上可表示為：

式（1）中x、y為束流的實(shí)際位置；I為束流流強(qiáng)；σ是束團(tuán)長(zhǎng)度；Aelse是其他因素綜合之后的影響因子。通常用4個(gè)電極輸出信號(hào)的差和比來(lái)計(jì)算束流位置：

式（2）中歸一化的電信號(hào)U、V的值分別反映束流在X方向和Y方向的位置信息。經(jīng)過(guò)式（2）計(jì)算，流強(qiáng)的影響被消除。對(duì)于相對(duì)論粒子（β≈ 1），電子學(xué)的工作頻率與束團(tuán)長(zhǎng)度的影響也被消除了。束流位置與歸一化的電信號(hào)的關(guān)系可表示為［9-10］：

階數(shù)n可以取不同的值，以n=1為例，束流位置可表示為：

式中：A0（B0）為水平（豎直）方向的電中心的偏移，產(chǎn)生的原因是電極差異、線纜和電子學(xué)通道增益差異。它可以用準(zhǔn)直技術(shù)（Beam Base Alignment）來(lái)消除［11-12］。A2和B1則分別是兩個(gè)方向的耦合，通常也是與束流位置相關(guān)的函數(shù)，但是在管道中心時(shí)近似為0，而A1和B2則是BPM最重要的一個(gè)參數(shù)，稱之為位置靈敏度系數(shù)，通常用k表示，其中水平與豎直方向的位置與靈敏度系數(shù)滿足：

為了精確測(cè)量束流位置，BPM在安裝之前，通常會(huì)在工作臺(tái)上校準(zhǔn)位置靈敏度系數(shù)。在標(biāo)定時(shí)，向金屬絲或天線上輸入射頻信號(hào)來(lái)模擬束流電場(chǎng)，改變天線的位置記錄信號(hào)大小的差和比，利用式（3）可以同時(shí)得到包括靈敏度系數(shù)在內(nèi)的多個(gè)系數(shù)［7］。實(shí)驗(yàn)上通常會(huì)對(duì)靈敏度系數(shù)的測(cè)量值和設(shè)計(jì)的理論值相比較。

1.2 邊界元法理論

使用邊界元法計(jì)算BPM的位置靈敏度系數(shù)，先要確定電極上的感應(yīng)電荷，電極一般對(duì)稱分布在束流管道上。對(duì)于相對(duì)論束流，束流的電場(chǎng)收縮于橫截面上，電極上感應(yīng)電荷的求解問(wèn)題可以簡(jiǎn)化為管道邊界上感應(yīng)電荷分布的問(wèn)題，即由束流管道包圍的封閉區(qū)域的靜電問(wèn)題［9，13-14］，其二維示意圖如圖2（a）所示，電荷ρ0表示束流，被真空室邊界Г包圍。

圖2 二維靜電問(wèn)題求解(a)和管道被分成N小段的離散邊界Г (b)示意圖Fig.2 Two-dimensional electrostatic problem (a) and discretization of the boundary Г by a finite number N of small line segments (b)

描述這類問(wèn)題的基本方程是泊松方程，標(biāo)量勢(shì)φ（）和電荷密度ρ（）滿足：

為了找到標(biāo)量勢(shì)的解，根據(jù)兩個(gè)標(biāo)量場(chǎng)Φ、Ψ的格林方程：

其中：V是一個(gè)被封閉曲面S包圍的任意三維體積。選擇：

其中：φ（）是由式（6）描述的標(biāo)量勢(shì)，而G（）是相應(yīng)的格林函數(shù)，它滿足：

將式（6）和式（9）代入到式（7）中，可以得到標(biāo)量勢(shì)的積分表示：

而φ關(guān)于表面法線n?的導(dǎo)數(shù)可以表示為：

此外，假設(shè)束流管道表面S是理想導(dǎo)體，標(biāo)量勢(shì)滿足：

它的積分形式為：

在二維情況下，格林函數(shù)可以通過(guò)直接求解方程（9）來(lái)確定，即：

假設(shè)包含電荷ρ（x，y）的二維區(qū)域用Ω表示，束流管道用Г表示，如圖2所示，靜電感應(yīng)的標(biāo)量勢(shì)的積分形式可表示為［13］：

如果將束流看作點(diǎn)電荷，則電荷分布為：

在這個(gè)條件下，可以直接對(duì)面積Ω進(jìn)行積分，有：

這個(gè)等式是所有進(jìn)一步計(jì)算的基礎(chǔ)。

式（17）可以用數(shù)值法 求解。將邊界劃分成一個(gè)個(gè)小的線段［14］，如圖2（b）所示，第i個(gè)線段的電勢(shì)φi可表示為：

假設(shè)邊界線段單元足夠小，則電荷密度σj（s）在該同一個(gè)線段單元內(nèi)近似為常數(shù)，方程（18）可寫(xiě)為：

我們將管道的邊界Г劃分成N個(gè)邊界元素Гj，方程（19）可表示為：

其中的矩陣［φi］、［Gi0］是N維列向量，而矩陣［Gij］是一個(gè)N×N矩陣。束流管道的邊界Г的電勢(shì)為零。通過(guò)矩陣求逆，方程（22）可以改寫(xiě)為：

通過(guò)求解該矩陣方程，就可以得到管道上的感應(yīng)電荷密度分布σj。

2 感應(yīng)電荷密度的計(jì)算

不同加速器管道的形狀各不相同，圖3是不同裝置的BPM橫截面示意圖。對(duì)于圓形或橢圓形的管道，可根據(jù)方位角θ來(lái)劃分邊界，首先以圖3（a）所示的管道為例，如果將管道分成N=3 600段，束流位于管道中心時(shí)，管道上的感應(yīng)電荷為均勻分布，每段上的歸一化電荷密度將為1/3 600= 2.78×10-4，每段邊界元的長(zhǎng)度為19 μm，遠(yuǎn)小于典型的CST網(wǎng)格劃分的尺寸。當(dāng)邊界元長(zhǎng)度與機(jī)械加工精度（~10 μm）、紐扣電極定位精度（~20 μm）和束流位置測(cè)量精度（~10 μm）同一量級(jí)時(shí)，截?cái)嗾`差對(duì)結(jié)果的影響幾乎可以忽略。現(xiàn)在假設(shè)束流位于（x0=y0=1 mm），管道上的電荷將不再均勻分布，這時(shí)需要用式（23）計(jì)算其電荷分布，則σ隨θ的分布如圖4（b）中黑線所示，其中紐扣電極（即拾取電極）直徑Φb=8 mm，4個(gè)紐扣電極分別位于45°、135°、225°、315°處，橫坐標(biāo)θ為各單元對(duì)應(yīng)的角度，可以看出θ= 45°處電荷密度最大；在相同條件下，對(duì)于圖3（b）中的橢圓形管道，束流位于管道中心時(shí)，最大值出現(xiàn)在θ= 90°和θ= 270°處，當(dāng)束流位于x0=y0= 1 mm處，則最大值出現(xiàn)在θ= 83°（整個(gè)管道上離束流最近）處，如圖4（b）中紅線所示。對(duì)于八邊形、跑道型或其他不規(guī)則的管道，計(jì)算的方法與圓或橢圓類似，將管道分成N段，確認(rèn)每一小段的坐標(biāo)后，用式（23）計(jì)算。注意如果是根據(jù)方位角來(lái)劃分單元，因?yàn)榫喙艿乐行木嚯x不同，劃分的線段元的長(zhǎng)度是不同的。橫截面為八邊形的北京正負(fù)電子對(duì)撞機(jī)BPM的感應(yīng)電荷分布如圖4（c）所示（假設(shè)束流分別位于y0=0，x0= 0 mm、10 mm、20 mm、30 mm處），其中P1~P4對(duì)應(yīng)的點(diǎn)如圖3（c）所示，與圓和橢圓不同的是，八邊形管道橫坐標(biāo)采用的是邊界上距離起始點(diǎn)的路程，其中起始點(diǎn)為P0（s=0）點(diǎn)。因?yàn)閳D形是上下對(duì)稱的，圖4（c）中僅顯示了管道上半部分的感應(yīng)電荷分布。當(dāng)束流位于管道中心時(shí)，兩個(gè)紐扣電極之間的中點(diǎn)M1處的電荷密度最大，而當(dāng)束流沿著水平軸向右移動(dòng)時(shí)，最大點(diǎn)逐漸往P1方向移動(dòng)。

圖3 高能光源儲(chǔ)存環(huán)BPM[15-16]圓形管道(a)、高能光源增強(qiáng)器BPM橢圓管道(b)和北京正負(fù)電子對(duì)撞機(jī)BPM[17]八邊形管道(c)的橫截面示意圖Fig.3 Cross section diagram of the HEPS storage ring BPM[15-16] on circular pipe (a), HEPS booster BPM on elliptical pipe (b), and BEPC II BPM[17] on octagonal pipe (c)

圖4 HEPS儲(chǔ)存環(huán)BPM截面圖(a)、HEPS儲(chǔ)存環(huán)和增強(qiáng)器BPM的感應(yīng)電荷分布圖(b)和北京正負(fù)電子對(duì)撞機(jī)八邊形BPM感應(yīng)電荷分布圖(c)（彩色見(jiàn)網(wǎng)絡(luò)版）Fig.4 Cross sectional diagram of the HEPS storage ring BPM (a), induced charge distributions for the storage ring and booster BPM (b), and induced charge distributions for the BEPCII octagonal pipe BPM (c) (color online)

3 位置靈敏度系數(shù)與UV Mapping圖

通過(guò)邊界元法計(jì)算得到管道上的電荷分布之后，特定位置電極上的感應(yīng)電荷Qi（i=a，b，c，d）可以通過(guò)對(duì)該電極覆蓋的邊界線段元的電荷密度求和得到。電極上的信號(hào)強(qiáng)度正比于感應(yīng)電荷的大小，用Qi替代等式（2）中的Vi（i=a，b，c，d），可得到差和比U和V，畫(huà)出Mapping圖。相應(yīng)的位置靈敏度系數(shù)可由差和比U、V對(duì)束流位置x0、y0的導(dǎo)數(shù)確定：

互為倒數(shù)的位置靈敏度S（%·mm-1）和靈敏度系數(shù)k（mm）是BPM最重要的特征值。在確定了位置靈敏度系數(shù)后，通過(guò)等式（5）可計(jì)算束流位置。

為了檢驗(yàn)該方法，首先對(duì)圓形管道（HEPS儲(chǔ)存環(huán)BPM）進(jìn)行了計(jì)算，為了進(jìn)行比較，也對(duì)四個(gè)電極A、B、C、D分別位于0°、90°、180°、270°的正交型BPM進(jìn)行計(jì)算，很多不受同步光影響的直線加速器上的BPM就是此種類型。對(duì)于正交型BPM的差和比U，V分別為：

此時(shí)稱為兩電極算法；電極A、B、C、D分別位于45°、135°、225°、315°時(shí)，對(duì)應(yīng)四電極算法，U、V值使用式（2）進(jìn)行計(jì)算。圖5是兩種方法下，束流位置的橫坐標(biāo)x0與縱坐標(biāo)y0與U、V組成Mapping圖，以BPM幾何中心為原點(diǎn)，束流位置（x0，y0）取值范圍為±7 mm×±7 mm，相鄰取值點(diǎn)間隔為1 mm，紐扣電極半徑rb= 4 mm，管道半徑Rb=11 mm。從圖5中可以看出，在相同的管道尺寸之下，兩電極算法的線性區(qū)域更大。

圖5 兩電極(a)與四電極(b)算法Mapping圖Fig.5 Mapping of a round pipe with two-electrode (a) and four-electrode (b) algorithms

對(duì)于規(guī)則的圓形管道，兩電極與四電極在管道中心的靈敏度系數(shù)ki（0，0）還可以分別由以下解析公式進(jìn)行計(jì)算［13，18-19］：

式中：Rb是束流管道半徑；α為紐扣電極所占的張角；θ為紐扣電極中心所處的方位角。

表1為不同方法得到的BPM位置靈敏度系數(shù)對(duì)比。邊界元法的位置靈敏度系數(shù)是經(jīng)過(guò)邊界元法計(jì)算U、V值后由式（2）、（24）和（25）計(jì)算得到，解析值是根據(jù)式（28）得到的計(jì)算值。由表1可知，用邊界元方法計(jì)算所得的k值與解析值十分接近，相對(duì)誤差在1%左右。表中最后一列的實(shí)驗(yàn)值是通過(guò)校準(zhǔn)天線的工作臺(tái)測(cè)量所得，與邊界元法計(jì)算值的相對(duì)誤差在2%左右。儲(chǔ)存環(huán)BPM的測(cè)量數(shù)據(jù)來(lái)自于加工的一套驗(yàn)證樣機(jī)；增強(qiáng)器BPM的數(shù)據(jù)來(lái)自已經(jīng)加工完成的50套正式BPM，k值的波動(dòng)在±0.03 mm。50個(gè)增強(qiáng)器BPM的kx、ky的測(cè)量數(shù)據(jù)也證明了邊界元法的精度在1%量級(jí)。

表1 不同方法得到的各BPM位置靈敏度系數(shù)（mm）Table 1 Position sensitivity coefficients of BPMs obtained by different methods (mm)

4 位置靈敏度系數(shù)匹配與分布

邊界元法可用于BPM設(shè)計(jì)，比如對(duì)于八邊形管道，調(diào)節(jié)兩紐扣電極之間的距離，就可以實(shí)現(xiàn)kx與ky的匹配。對(duì)于橢圓形管道，可以調(diào)節(jié)紐扣電極的方位角來(lái)實(shí)現(xiàn)kx與ky的匹配。圖6是HEPS增強(qiáng)器BPM和BEPC II八邊形BPM的靈敏度系數(shù)匹配情況（k值為束流在管道中心時(shí)的位置靈敏度系數(shù)）。圖6（a）是橢圓管道中水平/豎直方向靈敏度系數(shù)隨紐扣電極A的方位角度θA變化情況，（另三個(gè)電極的方位角為θB=180°-θA，θC=180°+θA，θD= 360°-θA），從中可以看出，θA增大時(shí)，水平方向的k值增大、豎直方向k值減小。當(dāng)θA= 49°時(shí)kx與ky比較接近，分別是11.40 mm和11.59 mm。圖6（b）是八邊形管道中水平/豎直方向靈敏度系數(shù)隨紐扣電極A與紐扣電極B之間的距離D變化情況，對(duì)于八邊形管道，紐扣電極A、B關(guān)于y軸對(duì)稱，而C、D與A、B關(guān)于x軸對(duì)稱，設(shè)計(jì)時(shí)電極A和B的位置可以在圖3（c）中P2與P3組成的直線上移動(dòng)，當(dāng)A、B之間距離減小時(shí)，水平方向的k值減小、豎直方向k值增大，當(dāng)兩紐扣電極之間距離為33 mm時(shí)，kx和ky比較接近，分別為23.79 mm和23.91 mm。

圖6 邊界元法用于橢圓管道(a)和八邊形管道靈敏度系數(shù)匹配(b)Fig.6 Sensitivity coefficient matching for an oval pipe (a) and octagonal pipe (b) by the boundary element method

上面的計(jì)算僅分析了管道中心區(qū)域的位置靈敏度系數(shù)。從圖5的Mapping圖中可以看出，管道中心處線性度比較好，而遠(yuǎn)離管道中心處呈現(xiàn)明顯的非線性。從式（5）中可知，S和k是與束流位置相關(guān)的函數(shù)，管道中心處的S（0，0）和k（0，0）是設(shè)計(jì)時(shí)最關(guān)心的參數(shù)，但是非線性區(qū)的k值也同樣是BPM的重要特征參數(shù)。對(duì)于管道半徑為11 mm的儲(chǔ)存環(huán)BPM，兩電極和四電極算法所得的k值的等高線圖如圖7所示，其中圖7（a）表示4個(gè)電極分別位于0°、90°、180°、270°時(shí)，U、V采用兩電極算法，利用kx=dx/（U（x+dx，y）-U（x，y））計(jì)算出的kx分布，其中dx=0.1 mm；圖 7（b）表示四電極分別位于 45°、135°、225°、315°時(shí)，U、V采用四電極算法計(jì)算出的kx分布。從中可以看出，兩電極算法的k值等高線是橢圓形，而四電極算法的k值等高線是雙曲線型，無(wú)論是兩電極還是四電極算法，kx的最小值在90°和270°處，同理可知ky的最小值在0°和180°處。k值等高線圖可以直觀反映BPM在管道中不同位置的靈敏度差異。

圖7 圓形管道兩電極(a)與四電極(b)算法的k值等高線Fig.7 Contour maps of k for the two-electrode (a) and four-electrode (b) algorithms

橢圓形與八邊形的管道的k值等高線圖如圖8與圖9所示，利用kx= dx/ （U（x+ dx，y）-U（x，y））計(jì)算kx，利用ky= dy/ （U（x，y+ dy）-U（x，y））計(jì)算ky，其中dx= dy= 0.1 mm。從圖8中可知四電極算法的k值等高線為雙曲線，因?yàn)楣艿赖牟粚?duì)稱，僅能保證在中心區(qū)域的kx與ky相等，在線性區(qū)之外，kx與ky受管道形狀的影響而存在差異，為了便于比較，圖8（a）與圖8（b）中的標(biāo)尺保持一致，束流位置在±10 mm×±10 mm的區(qū)域內(nèi)，對(duì)應(yīng)的ky值均小于30 mm。八邊形的k值因邊界的影響形狀變得更為復(fù)雜一些，但是無(wú)論什么形狀的管道，在管道中心，信號(hào)對(duì)束流位置變化最為靈敏（S較大，k較?。?。

圖8 橢圓形管道水平靈敏度系數(shù)kx (a)和豎直靈敏度系數(shù)ky (b)等高線圖Fig.8 Contour maps of horizontal sensitivity coefficient kx (a) and vertical sensitivity coefficient ky (b) for an oval pipe

圖9 八邊形管道水平靈敏度系數(shù)kx (a)和豎直靈敏度系數(shù)ky (b)等高線圖Fig.9 Contour maps of horizontal sensitivity coefficient kx (a) and vertical sensitivity coefficient ky (b) for an octagonal pipe

5 結(jié)語(yǔ)

使用邊界元法計(jì)算了不同橫截面形狀BPM的位置靈敏度系數(shù)，對(duì)于圓形管道，邊界元法計(jì)算的靈敏度系數(shù)與解析公式一致，相對(duì)誤差在1%左右。對(duì)于沒(méi)有解析公式的橢圓形管道、八邊形管道，邊界元法計(jì)算的值與實(shí)驗(yàn)標(biāo)定的值也基本一致，誤差在2%左右。邊界元方法是計(jì)算BPM的位置靈敏度系數(shù)的有效方法。另外還以高能光源BPM和北京正負(fù)電子對(duì)撞機(jī)BPM為例，演示了如何通過(guò)邊界元法調(diào)節(jié)紐扣電極的位置來(lái)進(jìn)行靈敏度系數(shù)匹配。最后計(jì)算了圓形、橢圓形、八邊形管道的靈敏度系數(shù)在管道內(nèi)的分布，特別是遠(yuǎn)離中心的非線性區(qū)的情況。相比于CST等電磁場(chǎng)仿真軟件，邊界元法不需要太多的計(jì)算資源，就有極高的精度，在BPM設(shè)計(jì)與相關(guān)問(wèn)題分析上具有廣泛的應(yīng)用潛力。

作者貢獻(xiàn)聲明李勇負(fù)責(zé)文獻(xiàn)調(diào)研，論文撰寫(xiě)；何俊負(fù)責(zé)程序編寫(xiě)，數(shù)據(jù)分析，論文整體內(nèi)容負(fù)責(zé)；隨艷峰負(fù)責(zé)項(xiàng)目管理與基金支持；袁任賢負(fù)責(zé)研究思路建議；杜垚垚負(fù)責(zé)BPM電子學(xué)應(yīng)用；張婉負(fù)責(zé)數(shù)據(jù)與圖片處理；王徐建負(fù)責(zé)BPM機(jī)械；尹頔負(fù)責(zé)實(shí)驗(yàn)支持；于令達(dá)負(fù)責(zé)程序調(diào)試；劉智負(fù)責(zé)軟件支持；魏書(shū)軍負(fù)責(zé)束流位置算法；岳軍會(huì)負(fù)責(zé)基金支持；曹建社負(fù)責(zé)思路建議。