王樹(shù)青 廖 鵬 彭 麗 楊 超 印茂偉

1 MeV介質(zhì)壁質(zhì)子直線加速器數(shù)值模擬優(yōu)化

王樹(shù)青 廖 鵬 彭 麗 楊 超 印茂偉

(西南科技大學(xué)國(guó)防科技學(xué)院 四川綿陽(yáng) 621010)

采用自主開(kāi)發(fā)的電磁粒子模擬軟件對(duì)中國(guó)工程物理研究院研制的1 MeV介質(zhì)壁質(zhì)子加速器進(jìn)行數(shù)值仿真，并對(duì)影響介質(zhì)壁加速器的主要參數(shù)進(jìn)行了系統(tǒng)分析，得到了可靠的優(yōu)化方案。模擬結(jié)果顯示:保持加速梯度20 MV/m的情況下，HGI越厚，得到的束流能量峰值越高，束流值越大，需綜合考慮加速器尺寸和所需的束流和束能，以選取合適的HGI厚度;峰值電流隨加速梯度的減小而增大，即在特定要求，可以通過(guò)損失部分加速梯度得到更大的輸出電流;輸出束流峰值隨輸入電流大小線性增長(zhǎng)，且輸入束流每增加3 mA，對(duì)應(yīng)輸出峰值增加0.12 A，增長(zhǎng)比率為1:40;輸出束流峰值隨輸入束流半徑的增加而減少，隨法蘭內(nèi)半徑的增大而增加;隨著延時(shí)間距的增大，束流先增大后減小，最佳延時(shí)為21，21.5，21.7 ns，且電流值為不考慮延時(shí)的1.73倍。

介質(zhì)壁加速器 電磁粒子模擬 質(zhì)子治療

介質(zhì)壁質(zhì)子加速器進(jìn)行放射治療的優(yōu)點(diǎn)表現(xiàn)如下:首先，加速器體積大大縮小，采用固態(tài)元器件所建立的介質(zhì)壁加速器的加速梯度理論上可以超過(guò)100 MV/m，因此輸出質(zhì)子束的加速器的長(zhǎng)度在2 m左右，考慮到束線等因素，加速器的長(zhǎng)度可以控制在5 m左右，這對(duì)于輸出相同能量的傳統(tǒng)加速器而言(約70 m)，長(zhǎng)度大大縮短。其次，目前質(zhì)子治療系統(tǒng)由于加速器系統(tǒng)龐大不能實(shí)現(xiàn)移動(dòng)調(diào)節(jié)，因此采用結(jié)構(gòu)相對(duì)簡(jiǎn)單的治療機(jī)架實(shí)現(xiàn)病灶的定位調(diào)節(jié)，定位機(jī)架的成本超過(guò)加速器主體所需費(fèi)用。而采用介質(zhì)壁質(zhì)子加速器進(jìn)行放射治療，由于加速器主體較小，可以采用通過(guò)調(diào)節(jié)加速器主體來(lái)實(shí)現(xiàn)定位調(diào)節(jié)的方案，因此成本大為降低，初步測(cè)算不到目前傳統(tǒng)質(zhì)子加速器治療系統(tǒng)成本的一半。上述優(yōu)點(diǎn)將有利于提高治療效率，降低治療成本，使質(zhì)子放射治療技術(shù)為一般患者所接受，有利于質(zhì)子放射治療技術(shù)的應(yīng)用和推廣。鑒于醫(yī)用加速器誘人的市場(chǎng)前景，國(guó)外的相關(guān)企業(yè)，如:LLNL and UC Davis，Compact Particle Acceleration Corporation等，已開(kāi)展了相關(guān)的技術(shù)和應(yīng)用研究。在國(guó)內(nèi)，介質(zhì)壁加速器的研究工作主要集中在中國(guó)工程物理研究院流體物理研究所，其主要工作包括:微波介質(zhì)陶瓷材料的制備與加工、SiC光導(dǎo)開(kāi)關(guān)的研究工作。然而，對(duì)介質(zhì)壁加速器系統(tǒng)化數(shù)值仿真的研究報(bào)道甚少，筆者對(duì)此進(jìn)行了相關(guān)研究。

1 1 MeV介質(zhì)壁質(zhì)子加速器整體布局

1 MeV介質(zhì)壁加速器由注入器(離子源、匹配段)、踢束器(Kicker)、介質(zhì)壁加速段(DWA)、漂移段和束流診斷系統(tǒng)組成，如圖1所示［2］。離子源引出的質(zhì)子束經(jīng)匹配段(LEBT)聚焦后，由Kicker將一個(gè)脈沖寬度150 ns左右的束團(tuán)導(dǎo)入介質(zhì)壁加速段，在5 cm的長(zhǎng)度內(nèi)加速到約1 MeV。

圖1 介質(zhì)壁加速器總體布局Fig.1 The overall layout of DWA

2 電磁粒子模擬算法

網(wǎng)格生成器用于生成一個(gè)非均勻的、用于時(shí)域有限差分計(jì)算的網(wǎng)格，該網(wǎng)格由多個(gè)連續(xù)排列的Yee元胞組成。它能生成直角坐標(biāo)系和圓柱坐標(biāo)系下的三維網(wǎng)格以及直角坐標(biāo)系、圓柱坐標(biāo)系和極坐標(biāo)系下的二維網(wǎng)格。在進(jìn)行模擬之前，只需指定大致的網(wǎng)格步長(zhǎng)，并標(biāo)注那些關(guān)鍵點(diǎn)必須位于全網(wǎng)格點(diǎn)位置的區(qū)域，網(wǎng)格生成器就會(huì)自動(dòng)生成非均勻網(wǎng)格。如果某兩個(gè)區(qū)域所指定的網(wǎng)格步長(zhǎng)不同，則會(huì)在兩個(gè)區(qū)域連接處的幾個(gè)Yee元胞上采用漸變的網(wǎng)格步長(zhǎng)，以保證平滑的過(guò)渡。

時(shí)域電磁場(chǎng)求解器包括中心差分求解器和時(shí)偏求解器，用于在每個(gè)時(shí)間步計(jì)算和更新整個(gè)模擬區(qū)域內(nèi)的電磁場(chǎng)值，并處理各種電磁場(chǎng)邊界條件。在對(duì)電磁場(chǎng)邊界條件的處理上，也采用模塊化設(shè)計(jì)的思想，針對(duì)每一種邊界都提供功能實(shí)現(xiàn)部分和接口調(diào)用部分，并且盡量考慮到今后擴(kuò)展的需要，形成功能實(shí)現(xiàn)和接口調(diào)用的文件規(guī)范。另外，作為時(shí)域電磁場(chǎng)求解器的輔助部分，還提供了一些小的求解器:泊松求解器和靜磁場(chǎng)求解器用于提供電磁場(chǎng)初始條件，頻域求解器用于諧振腔本征模分析。

粒子求解器用于在每個(gè)時(shí)間步先計(jì)算帶電粒子的作用場(chǎng)，再推動(dòng)粒子運(yùn)動(dòng)，然后計(jì)算粒子引起的電荷密度場(chǎng)和電流密度場(chǎng)并處理各種帶電粒子邊界條件。在對(duì)帶電粒子邊界條件的處理上，也采用模塊化設(shè)計(jì)的思想，形成文件規(guī)范。粒子計(jì)算部分涉及多種帶電粒子，粒子既可以通過(guò)初始放置產(chǎn)生，又可以通過(guò)各種發(fā)射過(guò)程產(chǎn)生，其中后者是軟件實(shí)現(xiàn)的一個(gè)重點(diǎn)。

3 模型及參數(shù)介紹

圖2所示為介質(zhì)壁加速器剖面圖。它由3個(gè)HGI組成，中間利用金屬法蘭分隔，并饋入梯形電壓(如圖3所示)，上升沿和下降沿分別為2 ns，半高寬10 ns，頂寬8 ns，輸入電壓300 kV;入口段采用束發(fā)射，且輸入電流 30 mA，束壓40 keV，發(fā)射半徑5 mm，粒子束橫向均勻分布，法蘭內(nèi)半徑8 mm，厚度為2 mm;電介質(zhì)材料(即HGI)內(nèi)半徑17.5 mm，外半徑25 mm，HGI厚度為1.5 cm，介質(zhì)材料均勻;質(zhì)子束注入后21 ns時(shí)開(kāi)始饋入電壓;模擬總時(shí)間設(shè)為60 ns，徑向網(wǎng)格0.5 mm，角向60°，軸向1 mm，加速段前端4 cm，后端50 cm。

圖2 介質(zhì)壁加速器建模圖Fig.2 The modeling diagram of DWA

圖3 單個(gè)HGI電壓波形Fig.3 A single HGI voltage waveform

4 模擬結(jié)果及分析

圖4為加速段電場(chǎng)分布，由圖4(a)可知，HGI所在區(qū)域軸向電場(chǎng)大體是勻強(qiáng)電場(chǎng)。由圖4(b)可知，HGI所在區(qū)域徑向電場(chǎng):對(duì)于質(zhì)子而言，加速段入口處為橫向聚焦效應(yīng)，出口處呈橫向發(fā)散。

圖5為59 ns時(shí)質(zhì)子平均能量沿Z方向的演化，由圖5可知，隨著時(shí)間的推移，受到3個(gè)HGI同時(shí)作用的質(zhì)子束團(tuán)能量能達(dá)到0.935 MeV，且該束團(tuán)逐漸追趕未被加速和小于3個(gè)HGI同時(shí)作用的質(zhì)子束團(tuán)。

圖6為HGI出口5 cm處面電流大小隨時(shí)間演化，由圖6可知，在加速器出口5 cm，31.6 ns時(shí)電流值最大，達(dá)到1.2 A。

圖4 加速器加速段電場(chǎng)分布Fig.4 Electric field distribution of accelerator section in DWA

圖5 59 ns時(shí)質(zhì)子平均能量沿Z方向的演化Fig.5 Average energy of proton evolutions along the Z direction at 59 ns

圖6 HGI出口5 cm處面電流隨時(shí)間演化Fig.6 The plane current with time evolution at 5 cm of HGI export

表1為其它參數(shù)不變且保持加速梯度20 MV/m的情況下，加速器輸出束流大小和能量隨HGI厚度的變化關(guān)系。由表1可知，HGI越厚，加速電壓峰值越大，得到的束流能量峰值越高，束流值越大。這需要綜合考慮加速器尺寸和所需的束流和束能大小，以選取合適的單個(gè)HGI厚度。

表2為其它參數(shù)不變，僅通過(guò)調(diào)節(jié)HGI厚度來(lái)改變加速梯度，分析輸出束流大小隨加速梯度的變化關(guān)系。由表2可知，峰值電流隨加速梯度的減小而增大，即在特定要求，可以通過(guò)損失部分加速梯度得到更大的輸出電流。

表3為其它參數(shù)不變，改變輸入束流大小，分析輸出束流與輸入束流大小的關(guān)系。由表3可知，在距離加速段后5 cm位置處測(cè)得電流遵循如下規(guī)律:輸出束流峰值隨輸入電流大小線性增長(zhǎng)，且輸入束流每增加3 mA，對(duì)應(yīng)輸出峰值增加0.12 A，增長(zhǎng)比率為1:40。

表4為其它參數(shù)不變，改變輸入束流半徑，分析輸出束流與輸入束流半徑的關(guān)系。由表4可知，在距離加速段后5 cm位置處測(cè)得電流遵循如下規(guī)律:輸出束流峰值隨輸入束流半徑的增加而減少。

表5為其它參數(shù)不變，改變3個(gè)HGI的分隔法蘭內(nèi)半徑，研究輸出束流與法蘭內(nèi)半徑的變化。在距離加速段后5 cm位置處測(cè)得電流遵循如下規(guī)律:輸出束流峰值隨法蘭內(nèi)半徑的增大而增加。

表6為其它參數(shù)不變，根據(jù)圖6顯示的HGI軸向勻強(qiáng)電場(chǎng)的特性，將3個(gè)HGI兩端觸發(fā)電壓延時(shí)設(shè)為2:1，分析束流隨脈沖時(shí)序的關(guān)系。由表6可知，在距離加速段后5 cm位置處測(cè)得電流遵循如下規(guī)律:隨著延時(shí)間距的增大，束流先增大后減小，最佳延時(shí)為21，21.5，21.7 ns，且電流值為不考慮延時(shí)的1.73倍。

表1 改變HGI厚度束流和束能變化Table 1 The changes of the beam current and energy by changing the thickness of HGI

表2 改變HGI的加速梯度束流變化Table 2 The changes of the beam current by changing the accelerating gradients of HGI

表3 改變輸入束流大小輸出束流變化Table 3 The changes of the output beam by changing the input beam size

表4 改變HGI內(nèi)半徑束流變化Table 4 The changes of the output beam by changing the inner radius of HGI

表5 改變輸入束流半徑輸出束流變化Table 5 The changes of the output beam by changing the input beam radius

表6 改變3個(gè)HGI脈沖時(shí)序束流變化Table 6 The changes of the beam current by changing three pulse sequences of HGI

5 結(jié)論

經(jīng)模擬分析得到了1 MeV介質(zhì)壁質(zhì)子加速器主要參數(shù)的優(yōu)化方案:(1)在保持加速梯度20 MV/m的情況下，HGI越厚，得到的束流能量峰值越高，束流值越大，需綜合考慮加速器尺寸和所需的束流和束能，以選取合適的HGI厚度;(2)峰值電流隨加速梯度的減小而增大;(3)輸出束流峰值隨輸入電流大小線性增長(zhǎng)，且輸入束流每增加3mA，對(duì)應(yīng)輸出峰值增加0.12 A，增長(zhǎng)比率為1:40;(4)輸出束流峰值隨輸入束流半徑的增加而減少，隨法蘭內(nèi)半徑的增大而增加;(5)隨著延時(shí)間距的增大，束流先增大后減小，最佳延時(shí)為21，21.5，21.7 ns，且電流值為不考慮延時(shí)的1.73倍。

［1］ MALCOLM B.Some future directions for repetitive pulsed power［J］.IEEE Trans.Plasma Science，2002，30: 262－266.

［2］ NUNALLY W C.Critical component requirement for compact pulse power system architectures［J］.IEEE Trans.Plasma Science，2005，33(4):1262－1267.

［3］ 施衛(wèi)，馬湘蓉，薛紅.半絕緣GaAs光電導(dǎo)開(kāi)關(guān)的瞬態(tài)熱效應(yīng)［J］.物理學(xué)報(bào)，2010，59(8):5700－5705.

［4］ SAMPAYAN S，CAPORASO G，CHEN Y J，et al.Development of a compact radiography accelerator using dielectric accelerator technology［C］.Proc of the Particle Accelerator Conference，2005.716－718.

［5］ 楊超.基于負(fù)氫離子源的全三維PIC/MCC數(shù)值模擬算法研究［D］.成都:電子科技大學(xué)，2012.

［6］ 楊超，劉大剛，周俊，等.一種新型徑向三腔同軸虛陰極振蕩器全三維粒子模擬研究［J］.物理學(xué)報(bào)，2011，60(8):084102.

［7］ 王衛(wèi)，鄧建，夏連勝，等.基于大功率激光二極管的光導(dǎo)開(kāi)關(guān)導(dǎo)通特性［J］.強(qiáng)激光與粒子束，2014，26 (4):045102.

［8］ 劉毅，夏連勝，王衛(wèi)，等.玻璃陶瓷Blumlein線的傳輸特性［J］.強(qiáng)激光與粒子束，2014，26(4):045018.

［9］ 楊超，龍繼東，王平，等.潘寧源放電的全三維電磁粒子模擬/蒙特卡羅碰撞數(shù)值算法研究［J］.物理學(xué)報(bào)，2013，62(20): 205207.

Numerical Simulation and Optimization of 1 MeV Dielectric Wall Proton Liner Accelerator

WANG Shu－qing，LIAO Peng，PENG Li，YANG Chao，YIN Mao－wei

(School of National Defense Science and Technology，Southwest University of Science and Technology，Mianyang 621010，Sichuan，China)

In this paper，1 MeV dielectric wall proton accelerator which is developed by China Academy of Engineering Physics was numerically simulated by the self－developed electromagnetic particle simulation software.A reliable optimization scheme was obtained after the system analyzing of its major parameters.The simulation results are the following:With the accelerating gradient of 20 MV/m，the thicker the HGI is，the higher beam energy peak and beam value become.Accelerator’s size，the needed beam and beam energy should be considered synthetically in order to choose a suitable thick of HGI.The peak current increases with the accelerating gradient decreasing，which means larger output current would be obtained by losing part of accelerating gradient in the specific requirements.The output beam peak grows linearly with the input current and when the input current increases by every 3 mA，the output beam peak increases by 0.12 A，in which growth rate is 1:40.The output beam peak decreases with the input beam radius increasing and increases with the flange inner radius increasing.With the rise of delay interval，the beam current increases first and then decreases.The optimal delay are 21，21.5，21.7 ns and the current is 1.73 times the value of the one without considering the delay.

Dielectric wall accelerator;Particle－in－cell simulation;Proton therapy

TL53

A

1671－8755(2015)03－0094－05

2015－03－27

國(guó)家自然科學(xué)基金重點(diǎn)基金(11035004);中物院科學(xué)技術(shù)發(fā)展基金(2013A0402018);國(guó)家自然科學(xué)基金青年基金(51407169);西南科技大學(xué)博士基金項(xiàng)目(13ZX7106);浙江省科技計(jì)劃項(xiàng)目(2013C33073)。

王樹(shù)青(1992—)，男，本科生。通信作者:楊超，男，博士，研究方向?yàn)槟M粒子理論及算法。E－mail:ychao1983@126.com