李景春，龔學余，路興強，黃千紅

（南華大學 核科學技術(shù)學院，湖南 衡陽 421001）

電子回旋共振加熱（ECRH）是磁約束聚變裝置中等離子體加熱的重要方式。EC波不僅可很好地被等離子體吸收，且在ECRH 期間還可觀察到許多新的物理現(xiàn)象，如在ECRH 作用下m／n＝2／1 的撕裂模以及逃逸電子被有效抑制［1-2］，ECRH 作 用 后 湍 流 功 率 譜 明 顯 降 低［3］。另外，EC波能以很窄的波束傳播，因此發(fā)射器不必很靠近等離子體；近年來，隨著長脈沖、高功率的回旋管源的顯著發(fā)展［4-5］，ECRH 在等離子體加熱中越發(fā)得到重視。在ECRH 的研究中，第1步都是進行波跡計算，然后計算沿波跡方向波的衰減（等離子體對波功率的吸收），并確定功率沉積剖面［6］。起初為了計算的方便，人們在圓截面下計算等離子體的EC 波的傳播和吸收［7-8］，然而幾乎所有的現(xiàn)代托卡馬克都采用拉長的非圓截面位形，非圓截面下的磁場、溫度和密度分布與圓截面不同，它們影響著EC波在等離子體中的傳播和功率沉積，另外由于工程上的限制，人們無法得到任意平衡位形，非圓截面下的ECRH、ECCD研究有助于尋找最佳的平衡位形，以實現(xiàn)最優(yōu)的EC波功率吸收。因此，數(shù)值計算EC波在非圓截面等離子體中的波跡和功率沉積，并探究特定位形下等離子體參數(shù)對EC波功率吸收的影響顯得十分必要。

非圓截面等離子體的平衡位形一般通過數(shù)值求解平衡方程獲得，數(shù)值求解Grad-Shafranov方程可直接求解（R，Z）坐標下的兩維差分方程，也可利用變分方法，通過求矩方程給出近似解［9-10］。本文采用一種直接變分法求解平衡方程，得到托卡馬克的磁場，繼而得到溫度和密度分布，再將其插值到ray-tracing程序中，研究非圓截面下的磁場、發(fā)射角、等離子體密度對EC波在等離子體中的傳播軌跡和功率沉積的影響。

1 物理模型

射頻波在托卡馬克中的傳播由射線軌跡（波跡）方程描述。在滿足幾何光學近似下，波跡方程在環(huán)坐標（r，θ，φ）下可寫成如下形式：

式中：（r，θ，φ）中相應(yīng)的正則波矢k＝（kr，kθ，kφ）；kr是徑向波數(shù)；極向和環(huán)向波數(shù)分別表示為kθ＝m／r，kφ＝n／R，m 和n 分別為極向和環(huán)向模數(shù)，R＝R0＋rcosθ，R0為托卡馬克環(huán)的大半徑，D0（ω，k）為電子回旋波色散關(guān)系的實部。

對于電子回旋波，由于波相速度遠大于電子熱速度，應(yīng)用冷等離子體色散關(guān)系，電子回旋波的色散方程可表示為：

式中：平方根項前面的正號對應(yīng)于尋常波（O-波）的情況，波的電場矢量平行于磁場，尋常波滿足折射定律；負號對應(yīng)于非尋常波（X-波）的情況，波的電場矢量位于垂直于磁場的平面上，由縱向分量和橫向分量所組成，非尋常波不滿足折射定律。N＝kc／ω，X＝ω2p／ω2，Y＝ωce／ω，ωp＝（nee2／meε0）1／2為等離子體頻率，ωce＝eB／me為電子回旋頻率?？煽闯觯臻g磁場的分布能改變色散方程，以致改變波跡。

對于波功率沉積的求解，由于托卡馬克磁場中捕獲粒子在香蕉軌道中來回運動，對粒子分布產(chǎn)生一定影響，于是采用反彈平均的Fokker-Planck方程來描述粒子運動［11］：

其中，μ0＝cosθ，其他各量的具體表示參見文獻［11］。該方程是一個積分微分方程，求解得出電子分布函數(shù)后，波功率沉積可由下式計算：

對于磁場位形的求解，本文利用能量原理，采用直接的變分方法，對平衡方程進行求解［12-13］。給定等離子體壓強和電流分布后，便可求解得到極向磁通ψp，以得到非圓截面下的磁場位形。

得到磁場后，通過下列各式可進一步求出等離子體中的密度和溫度分布［14］：

其中：ψpa為邊界處的極向磁通；αte、αti、αdn均可取0.5、1.0、1.5，分別對應(yīng)不同的溫度、電子密度剖面分布。

2 數(shù)值計算與結(jié)果討論

對平衡方程進行數(shù)值求解，采用與文獻［12］相同的參數(shù)，得到磁場分布如圖1所示。結(jié)果與文獻［12］一致。

計算HL-2A 上的環(huán)向、極向磁場分布，結(jié)果如圖2所示。HL-2A的參數(shù)取為：R0＝1.65m，a＝0.4m，Ip＝450kA，Bφ0＝2.42T，κ＝1.3，δ＝0.3。

為了使程序能計算托卡馬克不同磁場位形下的EC 波波跡，將等離子體每點的磁場、溫度、密度作為輸入變量，映射到（r，θ）網(wǎng)格點上，然后在計算時根據(jù)需要進行數(shù)值插值。

圖1 Z＝0時，環(huán)向、極向感應(yīng)磁場隨R 的剖面分布Fig.1 Profile of poloidal and toroidal induction components against Rfor Z＝0

圖2 Z＝0時，HL-2A 上環(huán)向、極向感應(yīng)磁場隨R 的剖面分布Fig.2 Profile of poloidal and toroidal induction components against Rfor Z＝0on HL-2A

為了保證在邊界條件μ＝-1和μ＝1上的對稱性，加入點j＝1和j＝N2，且有離散化分布函數(shù)用fi，j＝f（pi，μj）代替，方程最終的離散化形式為：

此方程實際為一個大型九對角帶狀方程組，本文采用矩陣三角分解法中的追趕法的一種改進方法——近似因子分解法加以求解。

采用HL-2A 裝置的參數(shù)，其中電子回旋波波功率和頻率分別?。篜 ＝0.5 MW，ω＝68GHz；等離子體密度?。簄e0＝2.5×1019m-3，nea＝0.5×1019m-3；等離子體溫度?。篢eo＝1.1keV，Ti0＝0.65keV，Ti（e）a＝50eV。圖3為相同中心磁場、等離子體溫度和密度下，EC波O 模分別從非圓截面、圓截面托卡馬克等離子體頂部發(fā)射時的波跡，其中對應(yīng)極向發(fā)射角θ＝185°，環(huán)向發(fā)射角與平行折射率相關(guān)，φ＝arcsin N‖，分別取N‖＝0、0.574、0.643。由圖3可看出，相比圓截面，非圓截面下EC 波波跡會明顯地向弱場側(cè)偏移。說明磁場、溫度和電子密度分布的共同改變對波跡產(chǎn)生了影響。為研究單個變量分布對波跡產(chǎn)生的影響，固定磁場分布，改變αte、αti、αdn以給定不同溫度、密度分布，且每次計算只改變單個變量分布，發(fā)現(xiàn)溫度、電子密度分布的改變并未對波跡產(chǎn)生明顯影響，可見非圓截面等離子體中主要是磁場分布影響EC波波跡。

本文也研究了非圓截面下，等離子體中心溫度和電子密度對EC波O 模波跡和功率沉積的影響，計算發(fā)現(xiàn)：中心溫度對EC 波波跡影響不大，而中心電子密度對波跡影響明顯。圖4a、b分別為θ＝185°、φ＝15°時，EC 波在等離子體中心的電子密度ne0＝2.0×1019、2.5×1019、3.0×1019m-3時的波跡和功率沉積。由圖可知，隨著中心電子密度的增大，波跡穿過等離子體中心后，會向等離子體弱場側(cè)偏移。這是因為中心密度增大會導致等離子體折射率變化率增大，如圖4c所示。在中心密度大的等離子體中，EC波穿過等離子體中心后，因為前后兩點的相對折射率大，波跡會發(fā)生更大角度偏折，從而偏向弱場側(cè)。從頂部發(fā)射時，波的功率沉積較低，且隨中心密度增大而降低，這一結(jié)果與通常的理論分析有所區(qū)別，對于O1模，等離子體密度的增大使其左旋極化逐漸轉(zhuǎn)變?yōu)槠叫袠O化，而O1的吸收主要來自平行極化，波的吸收應(yīng)隨等離子體密度增大而增大。這里的反常數(shù)值結(jié)果可能是由于中心電子密度增大，折射損失率增大，從而導致波的功率沉積隨中心密度增大而降低。

本文還計算了EC波O 模從中平面弱場側(cè)發(fā)射時的波跡和功率沉積，結(jié)果如圖5 所示。圖5a和b分別為EC 波在不同平行折射率下波跡的小截面圖和中平面俯視圖，對應(yīng)極向發(fā)射角θ＝185°，平行折射率分別取N‖＝0、0.574、0.707、0.766、0.866。由圖可知，當N‖＝0時，電子回旋波在等離子體中為近直線傳播；當N‖≠0時，波跡線向弱場側(cè)彎曲，隨著平行折射率的增大，彎曲幅度變大，甚至折回弱場側(cè)穿出，這種情形下波未到達共振層便傳出，對波在等離子體中的能量沉積很不利，因此是實驗上應(yīng)避免的。圖6a為θ＝185°，φ＝14°、18°、22°時波功率沉積的徑向分布，圖6b為θ＝180°、185°、195°時，EC波的單次吸收（即波饋入等離子體后第1次經(jīng)過共振層時所損失的能量與總饋入能量之比）隨環(huán)向發(fā)射角的變化。計算表明，波功率沉積隨極向發(fā)射角的增大而降低，環(huán)向發(fā)射角取10°左右時，波功率沉積達到最大，對于θ＝180°，環(huán)向發(fā)射角在14°以內(nèi)，單次吸收在90%以上。為獲得較高的加熱效率，應(yīng)盡可能保持低的極向和環(huán)向發(fā)射角。

圖3 EC波O 模從頂部發(fā)射時，分別在非圓截面、圓截面托卡馬克等離子體的傳播軌跡Fig.3 Ray trajectories for O-mode of ECW launched from top in non-circular and circular Tokamak plasma

圖4 不同中心電子密度下的EC波波跡（a）、波功率沉積（b）和等離子體折射率（c）Fig.4 Ray trajectories（a）for O-mode of ECW，power deposition（b）and plasma refractive index（c）with different central electron densities

圖5 EC波O 模從弱場側(cè)（LFS）發(fā)射時波跡的小截面圖（a）和中平面俯視圖（b）Fig.5 Ray trajectories for O-mode launched from LFS In poloidal section（a）and equatorial plane（b）

圖6 極向發(fā)射角為185°時波的功率沉積分布（a）和極向發(fā)射角為180°、185°和190°時EC波的單次吸收隨環(huán)向發(fā)射角的變化（b）Fig.6 Profile of power deposition of EC wave at poloidal injection angle of 185°（a）and EC wave single absorption versus toroidal injection angle at poloidal injection angle of 180°，185°and 190°（b）

3 結(jié)論

本文采用直接變分方法求解平衡方程，開發(fā)出一種等離子體平衡方程求解程序。并將此程序與求解波跡的程序進行耦合，數(shù)值計算了電子回旋波在非圓截面托卡馬克等離子體中的傳播和功率沉積。數(shù)值計算結(jié)果表明：當EC波從頂部發(fā)射時，相比于圓截面下的波跡，受磁場分布影響，非圓截面的波跡會向弱場側(cè)偏移；增大等離子體中心電子密度，波跡也會偏向弱場側(cè)，相應(yīng)的波功率沉積有所降低；當EC 波O模從弱場側(cè)發(fā)射時，N‖＝0 時，電子回旋波在等離子體中為近直線傳播；N‖≠0 時，波跡線向弱場側(cè)彎曲，隨著平行折射率的增大，彎曲幅度變大，甚至折回弱場側(cè)穿出，此時波功率沉積隨極向發(fā)射角的增大而降低，環(huán)向發(fā)射角約取10°時，波功率沉積達到最大。為獲得較高的加熱效率，應(yīng)盡可能保持低的極向和環(huán)向發(fā)射角。

［1］ DUAN X R，DING X T，DONG J Q，et al.Overview of experimental results on HL-2A［J］.Nucl Fusion，2009，49：104012.

［2］ YANG J W，ZHANG Y P，LI X，et al.Suppression of runaway electrons during electron cyclotron resonance heating on HL-2A Tokamak［J］.J Plasma Phys，2010，76：75-85.

［3］ YAN L W，DUAN X R，DING X T，et al.Overview of experimental results on the HL-2A Tokamak［J］.Nucl Fusion，2011，51：094016.

［4］ ROBINSON D C，COX M，LLOYD B，et al.Progress with heating and current drive technologies［J］.Fusion Eng Des，1999，46（2）：355-370.

［5］ THUMM M.MW gyrotron development for fusion plasma applications［J］.Plasma Phys Con-trolled Fusion，2003，45（12A）：A143-A161.

［6］ SHI Bingren，SHI Yingtian.Analysis and calculation of ray tracing and wave power deposit profiles of electron cyclotron resonance heating in the HL-2ATokamak［J］.Nuclear Fusion and Plasma Physics，2007，27（1）：1-6.

［7］ 彭曉煒，龔學余，劉文艷.HL-2A Tokamak裝置中的電子回旋波電流驅(qū)動［J］.中山大學學報，2010，19（3）：52-54.PENG Xiaowei，GONG Xueyu，LIU Wenyan.Electron cyclotron current drive in HL-2ATokamak［J］.Acta Scientiarum Naturalium Universitatis Sunyatseni，2010，19（3）：52-54（in Chinese）.

［8］ 宋紹棟，劉永，張錦華，等.HL-2A 裝置ECRH波跡和功率沉積計算［J］.核聚變與等離子體物理，2007，27（2）：123-126.SONG Shaodong，LIU Yong，ZHANG Jinhua，et al.Calculation of ECRH ray tracing and power deposition on HL-2A Tokamak［J］.Nuclear Fusion and Plasma Physics，2007，27（2）：123-126（in Chinese）.

［9］ LAO L L，HIRSHMAN S P，WIELAND R M.Variational moment solutions to the Grad-Shafranov equation［J］.Phys Fluids，1981，24（8）：1 431-1 441.

［10］GUAZZOTTO L，F(xiàn)REIDBERG J P.A family of analytic equilibrium solutions for the Grad-Shafranov equation［J］.Physics of Plasmas，2007，14：112508.

［11］FARINA D.Notes on the computation of the EC power absorption and current drive in the linear and quasilinear regimes，F(xiàn)P 99／14［R］.Milano，Italy：Istituto di Fisica del Plasma，2000.

［12］LUDWIG G O.Direct variational solutions to the Grad-Schluter-Shafanov equation［J］.Plasma Phys Control Fusion，1995，37：633-646.

［13］LUDWIG G O.Direct variational solutions of the Tokamak equilibrium problem［J］.Plasma Phys Control Fusion，1997，39：2 021-2 037.

［14］段文學，吳斌.EAST 低混雜波功率沉積和電流驅(qū)動剖面控制的數(shù)值模擬［J］.計算物理，2011，28（3）：438-444.DUAN Wenxue，WU Bin.Simulation of lower hybrid wave power deposition and current drive control in EAST［J］.Chin J Comput Phys，2011，28（3）：438-444（in Chinese）.