邊界局域模引起鎢偏濾器靶板侵蝕和形貌變化的數值模擬研究?

黃艷孫繼忠桑超峰胡萬鵬王德真?

1)(大連工業大學信息科學與工程學院,大連 116034)2)(大連理工大學物理與光電工程學院,大連 116024)(2016年9月22日收到;2016年10月24日收到修改稿)

邊界局域模引起鎢偏濾器靶板侵蝕和形貌變化的數值模擬研究?

黃艷1)2)孫繼忠2)?桑超峰2)胡萬鵬2)王德真2)?

1)(大連工業大學信息科學與工程學院,大連 116034)2)(大連理工大學物理與光電工程學院,大連 116024)(2016年9月22日收到;2016年10月24日收到修改稿)

鎢材料在高瞬時熱流作用下的熔化、流動是國際熱核聚變實驗堆面壁材料最突出的問題.本文將熱傳導方程與Navier-Stokes方程結合,建立了二維流體動力學模型,研究在邊界局域模(ELM)強熱流轟擊下,鎢熔化層在表面張力、壓強梯度力、磁場力等作用下的流動,以及偏濾器靶板的侵蝕和形貌演化.結果表明,在ELM過程中,熔化層中的液體不斷地向邊緣區域流動,在打擊點區域形成一個熔池,在熔化層的邊緣區域形成類似山峰結構的凸起,加重了鎢偏濾器靶板的侵蝕.在空間分布為高斯形狀入射能流的作用下,鎢熔化層兩側的山峰結構是對稱的;當能流密度小于3000 MW.m?2時,表面張力對熔化層的流動起主要作用.本文在模型的數值求解中,采用交錯網格的方法進行離散,克服了液體表面追蹤的算法難點,保證了鎢偏濾器靶板侵蝕程度計算的準確性.

邊界局域模,鎢偏濾器靶板,熔化,流動

1 引 言

金屬鎢因為具有熱導率高、熔點高、濺射率低、氚滯留率低等優點,被選作國際熱核聚變實驗堆(ITER)的等離子體面壁材料,但鎢材料在瞬時事件,例如邊界局域模(ELM)、破裂、垂直位移事件等的高熱流作用下會熔化、汽化.熔化層在表面張力、磁場力、壓強梯度力等作用下會流動,導致材料形貌變化,影響材料的導熱能力,誘發電弧,甚至液體鎢會以液滴的形式飛濺出去.這些后果不但會縮短器壁的壽命,而且產生的鎢雜質(原子、塵埃、液滴)將嚴重威脅等離子體穩態運行[1].因此鎢的熔化問題(特別是在準穩態運行時ELMs熱流作用下)是以鎢作為ITER面壁材料最突出的問題[2].

目前,在多個聚變實驗裝置(例如TEXTOR和JET)中,人們已經觀察到鎢偏濾器靶板的熔化以及熔化層流動的現象[3?5].但是ITER由于加熱功率遠高于現有聚變實驗裝置,并且刮削層寬度遠小于現有等離子體,其等離子體與器壁相互作用的一系列問題(比如發生第一類ELMs和破裂時打到靶板的能量,全鎢偏濾器靶板的侵蝕和雜質的輸運等)與目前利用現有托卡馬克裝置外推預測的結果之間可能存在較大的差異[6],所以,理論上對該問題的分析研究就變得尤為重要.

Federici等[7?9]在一系列研究工作中建立了一維的數值求解模型,研究高熱流對器壁材料的侵蝕.但其工作緊緊局限于一維研究,這樣沒有辦法考慮能流密度空間分布的不均勻性.美國普渡大學的Hassanein等開發了HEIGHTS程序包,研究等離子體不穩定性事件中高能流對器壁材料的侵蝕[10?15].Hassanein等雖然在器壁材料侵蝕方面做了很多工作,但是對各種力作用于鎢熔化層帶動熔化層流動,加劇靶板侵蝕的相關問題并未開展研究.Bazylev等[16?21]建立了MEMOS模型,研究高能流對鎢材料的侵蝕.但是,Bazylev等在求解模型方程過程中,先沿靶板豎直方向求解熱傳導方程,然后求解Navier-Stokes方程,之后再沿水平方向求解熱傳導方程,而實際中熱傳導過程在豎直和水平方向是同時進行的,所以有必要對其求解過程進一步改進.另外,對能流密度小于3000MW.m?2的情況,尚未見到考慮熔化層流動時偏濾器靶板侵蝕的詳細研究.

本文在前期工作的基礎上[22,23],將熱傳導方程和流體力學的Navier-Stokes方程結合建立了二維流體動力學模型,研究在類似ITER運行參數下發生第一類ELM過程中能流密度小于3000 MW.m?2,鎢熔化層在表面張力、壓強梯度力、磁場力等的作用下流動時,偏濾器靶板的侵蝕和形貌變化,詳細地分析了ELM各參數對鎢偏濾器靶板侵蝕程度的影響,預測未來ITER發生ELM時鎢偏濾器靶板的侵蝕程度,為未來裝置的設計建造提供理論依據.另外,在模型的數值求解中,對于偏濾器靶板的溫度分布直接求解二維熱傳導方程;方程采用交錯網格的方法進行離散,成功解決了液體表面追蹤的算法難點,保證了鎢偏濾器靶板侵蝕程度計算的準確性.

2 物理模型與數值求解

未來ITER發生第一類ELMs時打到偏濾器靶板的能量預期為1—3 MJ.m?2,沉積時間為0.1—1 ms[24,25],在此范圍的能量作用下,厚度為1 cm的鎢偏濾器靶板的熔化厚度為幾十微米到上百微米[7,22,26],遠遠小于鎢靶板的厚度,這屬于淺水環流問題,因此可以只研究液體沿靶板表面的流動[27].如圖1所示,液體流動沿x軸方向(對應著實際托卡馬克裝置的極向),垂直偏濾器靶板豎直向下的方向為z軸正方向.圖1表示在空間分布為高斯形狀能流的作用下,鎢偏濾器靶板熔化,熔化層流動產生的形貌變化.模擬過程中考慮了熔化、汽化、熔化層流動和熱輻射效應.

研究中考慮鎢熔化層所受表面張力、壓強梯度力、磁場力,忽略熔化層的重力,熔化層中的鎢液體視為不可壓縮液體,有

表達式中ρ,C,v,T和k分別表示鎢的密度、比熱、流動速度、溫度和熱傳導系數,材料的熱傳導系數隨溫度改變,采用k=1/(aT+b)的表達形式[28],對于鎢材料a=1.70× 10?6m.W?1,b=6.41×10?3m.K.W?1[29];Q是體能量沉積項;p是作用于鎢液體表面的等離子體壓強和蒸汽壓強的總壓強.

圖1 熔化層及熔化層流動示意圖Fig.1.Sketch of melt layer and melt motion of the target.

液體表面為自由面,表面處既要滿足能量平衡條件,又要滿足力平衡條件:

固體和液體的交界處,使用經典斯蒂芬邊界條件[31]:

(6)式中下腳標s,l分別代表固體和液體,ks是固體部分的熱傳導率,kl是液體部分的熱傳導率,Lm表示熔化潛熱,Vm表示熔化層前端的熔化速度.計算過程中,靶板底端邊界條件取冷卻溫度350 K.

利用邊界條件(5),可以將方程(1),(2)整理成描述不可壓縮流體流動的一維St.Venant方程組的近似形式:

(8)和(9)式中h表示熔化層的厚度;ux表示水平方向流動速度對熔化層厚度所取的平均;ν=μ/ρ為鎢液體的運動黏度;kα表示表面張力系數對溫度的偏導數,是一個負值,表達式為kα=?α/?T;Vs表示液體表面的速度,包括兩部分機制,Vv為汽化速度,Vem為由于液體流動導致熔化層厚度增加或者減少而產生的液體表面速度;Jz表示z方向的電流密度;By表示環向磁場強度(本文y軸正方向設為環向磁場的反方向).(9)式右端各項的物理意義分別為:第一項為壓強梯度力,第二項表示由于液體質量增加所產生的摩擦力,第三項表示液體流動產生的黏滯力效應,第四項表示表面張力效應,第五項表示單位體積磁場力.

液體表面的汽化過程導致液體表面不斷向下移動,模型中引進了運動坐標z′:

其中,z=0的位置是靶板表面的初始位置.因為發生ELMs時,極向能流密度在寬為0.1 m的范圍內基本按高斯規律衰減,打擊點中心的能流密度最大[32],導致作用于鎢熔化層各點的壓強和表面張力不同.因此本文采用二維熱傳導模型研究極向ELM能流不均勻導致鎢熔化層在多種力作用下流動時,鎢偏濾器靶板的侵蝕和形貌變化.利用運動坐標(10),熱傳導方程(3)可以寫成如下形式:

(11)式中,Vs表示鎢靶板的表面速度(包括靶板表面的汽化速度以及由于液體流動導致熔化層厚度增加或者減少而產生的液體表面速度,見(8)式),ux為鎢液體在水平方向的流動速度,固體部分ux為零.模擬過程中假設作用在偏濾器靶板上的極向ELM能流空間分布為高斯分布,打擊點位置能流密度最大.靶板寬度取1 cm,靶板厚度取1 cm.

3 模擬結果

3.1 模型驗證

為了驗證物理模型的合理性,首先將本模型與已有模型進行對比.文獻[17]針對鎢材料和鈹材料分別在電子束裝置JEBIS和JUDITH中的熔化流動實驗進行了模擬研究,模擬結果與實驗結果非常符合,而且文獻[17]中所使用的模型已經被廣泛接受,因此具有一定代表性.本文針對文獻[17]中空間為高斯分布的能流作用于鎢靶板的例子進行了模擬,模擬參數與文獻中使用的參數相同,入射能流密度的峰值為1350 MW.m?2,沉積時間為1.8ms,kα= ?9.0×10?5N.m?1.K?1[17].計算結果表明,在能量沉積過程中鎢靶板汽化很弱,作用于鎢靶板的壓強梯度力很小,不會引起鎢液體流動,主要是表面張力引起鎢熔化層流動,在熔化層的邊緣處形成兩個像山峰一樣的凸起,如圖2(a)所示.能流脈沖結束后靶板兩端凸起的厚度約為8.35μm,文獻[17]中模擬的結果也是在熔化層邊緣處形成高度約為8μm的山峰,如圖2(b)所示,與本文模型得到的山峰高度基本一致.

在模型的數值求解過程中,T,h和ux三個變量均采用交錯網格的方法進行離散,優點是對速度分量無需再進行內插運算,同時有效地避免了常規算法處理液體流動這類問題易出現的數值噪聲[33].由于本文將坐標原點設在靶板的邊緣處,熔化層流動時,ELM打擊點兩側的流動速度方向不同,在求解過程中分段采用上風格式,成功地追蹤了熔化層流動過程中液體表面在各個時刻的位置,且所得的熔化層表面較光滑,說明本文采用的計算方法是成功的.

圖2 鎢靶板表面形貌(其中入射能流密度峰值為1350 MW.m?2,作用時間為1.8ms)(a)本文模擬結果;(b)內嵌圖,文獻[17]中的模擬結果Fig.2.Crater depth of tungsten under the peak power density of 1350 MW.m?2and the heat load du ration of 1.8 ms:(a)Calcu lated resu lts;(b)inset,profi le of the melt su rface in the Ref.[17].

3.2 表面形貌演化

未來ITER穩態運行時打到靶板的能流密度預期為5—20 MW.m?2[24,25],如果第一類ELMs間隙間沉積在偏濾器靶板的能流密度取10 MW.m?2,則在ELMs間隙間能流作用下,厚度為1 cm的鎢靶板表面溫度會達到1000 K左右[7,22]. 本文首先計算了能流密度峰值為2900 MW.m?2,持續時間為0.8 ms的ELM作用在初始溫度為1000 K的鎢偏濾器靶板上,考慮熔化層在表面張力和蒸汽壓強梯度力(在討論的能流密度范圍內,磁場力作用遠小于表面張力,所以本文忽略了磁場力)的作用下流動時,在熔化層兩端形成的像山峰一樣凸起的高度,偏濾器靶板的表面溫度以及熔化層厚度.計算過程中kα取?9.0×10?5N.m?1.K?1[17],同時與不考慮熔化層流動時偏濾器靶板的表面溫度和熔化層厚度進行了比較.

圖3給出了ELM過程中3個不同時刻靶板的表面形貌.研究表明,在能量沉積過程中鎢靶板汽化很弱,作用于鎢靶板表面的壓強梯度力很小,不能引起鎢液體流動.鎢靶板熔化之后,熔化層中的液體在表面張力的作用下會不斷地由打擊點區域向熔化層邊緣區域流動,在打擊點區域形成一個熔池,邊緣處形成形狀似山峰結構的凸起,計算時間越長,在鎢液體中產生的熔池越深,邊緣處形成的山峰高度越高.表1列出了從t=0.6 ms到t=0.95 ms,每過0.05 ms靶板表面凸起的最大值,也列出了該時刻與0.05 ms之前時刻熔化層凸起最大值的增加量.從表1數據可以看出,ELM能流與靶板作用時間越長,熔化層兩端凸起厚度的增加速度越快.因為在本文選取的計算參數下,靶板表面的汽化量很小,靶板表面形貌變化主要由熔化層流動造成.ELM能量沉積時間越長,靶板熔化厚度越大,熔化層流動速度越大,由熔化層流動導致熔化層厚度增加或減少的幅值越大.

圖3 (網刊彩色)靶板表面形貌演化,其中,能流密度峰值為2900 MW.m?2,kα取?9.0×10?5N.m?1.K?1,初始溫度為1000 KFig.3.(color on line)Temporal evolution of the tungsten melt surface.Incident heat,2900 MW.m?2,kα,?9.0 × 10?5N.m?1.K?1,and initial sample temperatu re,1000 K.

表1 各時刻熔化層兩端凸起高度的最大值以及0.05 ms時間內熔化層凸起高度的變化值Tab le 1.Themaximumheight of the lumpon themelt layer edge,and increment of the lumpin 0.05 ms.

圖4分別給出了考慮熔化層流動和不考慮熔化層流動,熔化區域靶板表面在三個不同時刻的溫度分布.從圖4可以看出,t=0.6 ms時,熔化層表面形貌變化較小,考慮熔化層流動和不考慮熔化層流動靶板表面溫度相差不大.t=0.7 ms時刻,考慮熔化層流動和不考慮熔化層流動靶板表面溫度開始不同.t=0.8 ms時刻,考慮熔化層流動時,x=0.46—0.54 cm之間區域靶板表面溫度明顯低于不考慮熔化層流動時靶板表面的溫度分布,打擊點處溫度相差了約110 K.在x=0.39—0.46 cm之間的區域和x=0.54—0.61 cm之間的區域,考慮熔化層流動時靶板表面溫度高于不考慮熔化層流動時靶板的表面溫度,這是由于鎢液體在流動的過程中將熱量帶到了熔化層的邊緣區域.另外,從圖4的結果還可以看出,熔化層的流動也降低了靶板表面各點的溫度梯度.

圖4 (網刊彩色)偏濾器靶板的表面溫度的演化,其中,能流密度峰值為2900 MW.m?2,kα取?9.0×10?5N.m?1.K?1,初始溫度為1000 KFig.4.(color on line)Temporalevolu tion of the surface temperatu re with the incident heat of 2900 MW.m?2,kαof?9.0×10?5N.m?1.K?1,and the initial sample temperatu re of 1000 K.

圖5分別給出了考慮熔化層流動和不考慮熔化層流動,不同時刻偏濾器靶板的熔化厚度.從圖5可以看出:由于熱量隨著時間不斷向下傳導,計算時間越長,靶板的熔化厚度越大;而且,不考慮熔化層流動計算所得的熔化層前沿的位置要小于考慮熔化層流動熔化層前沿的位置,與文獻[34]的規律一致,計算時間越長,二者相差越大.在t=0.6 ms,不考慮熔化層流動時,熔化層前沿位于z=77.90μm;考慮熔化層流動時,熔化層前沿位于z=79.89μm處,二者相差了約1.99μm;到t=0.7 ms時,兩種情況下熔化層前沿的位置相差了約2.24μm;到t=0.8 ms時,此參數約為2.62μm.

以上結果說明,在第一類ELM熱流與靶板的相互作用過程中,即使在低ELM能量作用下,作用在熔化層上的力也會導致熔化層流動,使得材料形貌發生變化,增加材料的熔化厚度.

圖5 (網刊彩色)靶板熔化厚度的演化,其中,能流密度峰值為2900 MW.m?2,kα取?9.0× 10?5N.m?1.K?1,初始溫度為1000 KFig.5.(color on line)Temporal evolu tion of the melt front with the incident heat of 2900 MW.m?2,kαof?9.0×10?5N.m?1.K?1,and the initial sample temperatu re of 1000 K.

3.3 kα幅值的影響

靶板中雜質濃度等因素將會影響熔化層的表面張力[17],因此有必要分析不同幅值的表面張力對鎢偏濾器靶板侵蝕程度的影響.圖6給出了初始溫度為1000 K的鎢偏濾器靶板,在入射能流密度峰值為2000 MW.m?2,沉積時間為0.8ms的ELM能量作用下,表面張力系數對溫度的偏導數分別取四個不同值時(kα值取自文獻[17])靶板表面的形貌分布.從圖6可以看出,對于kα的四個不同幅值,表面張力都引起了鎢液體層流動,增加了鎢靶板的侵蝕程度.kα越大,由于鎢液體流動,在熔化層邊緣形成的山峰高度越高.因為kα越大,鎢液體的表面張力梯度越大,鎢液體流動的速度越快,所以兩端凸起的厚度越大.

圖6 (網刊彩色)靶板表面的形貌分布,其中,ELM能流密度峰值為2000 MW.m?2,沉積時間為0.8 ms,靶板初始溫度為1000 KFig.6. (color on line)Calcu lated crater depth of tungsten under ELMwith a peak heat flux of 2000 MW.m?2.The heat load duration is 0.8 ms,the initial sample temperature is 1000 K.

3.4 能流密度幅值的影響

未來ITER的第一類ELMs能量范圍預期為1—3 MJ.m?2,沉積時間為0.1—1 ms[24,25],本文計算了能流密度峰值分別為2000,2500,2900 MW.m?2的ELM(如圖7(a)所示)作用在初始溫度為1000 K的偏濾器靶板上,在持續時間為0.8 ms,kα取?9.0× 10?5N.m?1.K?1時,偏濾器靶板的表面形貌.計算結果表明:入射等離子體能流密度越大,由于鎢液體流動在熔化層邊緣處形成的凸起的厚度越大,這是由于入射能流密度越大,能流密度梯度越大,靶板熔化后,靶板表面溫度梯度越大,在kα取相同值時,鎢液體的表面張力梯度越大,鎢熔化層的流動速度越快;而且,能流密度越大靶板越先熔化、流動,所以在鎢熔化層兩端形成的山峰高度越高.

圖6和圖7(b)分別說明kα和ELMs能流密度會影響靶板表面形貌.進一步分析,圖6(紅色虛線曲線)和圖7(b)(紅色點劃線曲線)計算所得的熔化層兩側山峰的高度都約為5μm,兩幅圖的計算時間相同,靶板的初始溫度相同,但圖6的kα值為?2.5×10?4N.m?1.K?1,入射能流密度為2000 MW.m?2,圖7(b)的kα值為?9.0×10?5N.m?1.K?1,入射能流密度為2500MW.m?2.圖6的kα值約為圖7(b)kα值的2.78倍,而圖7(b)的能流密度僅約為圖6能流密度的1.25倍.可見,由鎢中雜質濃度變化等因素導致的鎢熔化層表面張力改變對靶板表面形貌的演化只起次要作用,而ELMs本身的特征參數變化對靶板表面形貌的改變起主要作用.

圖7 (網刊彩色)能流密度峰值分別為2000,2500和2900 MW.m?2時的(a)能流密度分布曲線和(b)偏濾器靶板表面形貌;ELM沉積時間為0.8 ms,kα為?9.0×10?5N.m?1.K?1,靶板初始溫度為1000 KFig.7. (color on line)(a)Incident heat loads and(b)temporal evolution of the tungsten melt surface calculated with peak heat fl uxes of 2000,2500 and 2900 MW.m?2.The heat load duration is 0.8 ms,kαis?9.0×10?5N.m?1.K?1,and the initial sample temperatu re is 1000 K.

3.5 討 論

為了更好地分析不同物理量對鎢偏濾器靶板侵蝕程度的影響,本文進一步比較了各種參數下偏濾器靶板的粗糙度,粗糙度定義采用(12)式的形式[35,36]:

(12)式中h(i)表示偏濾器靶板表面每點的高度,h表示偏濾器靶板表面的平均水平高度.圖8(a)—(c)分別給出了偏濾器靶板表面粗糙度隨入射能流密度、kα值和ELM持續時間的變化.圖8(a)的計算結果表明,ELM的持續時間為0.8ms,入射能流密度約為1650MW.m?2時,靶板表面粗糙度約為零.因為此時靶板剛開始熔化,鎢液體在熔化層兩側形成的山峰狀凸起非常小,約為10?4μm.隨著入射能流密度的增加,靶板表面粗糙度線性增加.圖8(b)和圖8(c)的結果表明,偏濾器靶板表面的粗糙度分別與液體表面張力系數對溫度的偏導數kα和ELM持續時間成正比.kα越大,ELM作用時間越長,鎢熔化層流動導致熔化層邊緣處山峰的高度越高,偏濾器靶板表面形貌變化越大.

圖8 鎢偏濾器靶板表面粗糙度與(a)入射能流密度、(b)kα值和(c)ELM持續時間的關系Fig.8.Surface roughnessvarieswith(a)incident heat,(b)kαand(c)ELMdu ration.

4 結 論

本文通過將熱傳導方程與Navier-Stokes方程結合,建立了二維流體動力學模型,分別從ELM本身參數和靶板角度,詳細地研究了在類似ITER運行參數下發生第一類ELM過程中,偏濾器靶板的侵蝕和形貌演化,得到如下結論.

1)ELM能流與靶板相互作用的過程中,鎢熔化層會發生流動,作用時間為0.8 ms,kα為?9.0 × 10?5N.m?1.K?1, 入射能流密度 大于2000 MW.m?2時,偏濾器靶板表面形貌發生明顯變化.入射能流密度若為2000 MW.m?2,靶板表面粗糙度約為1.1μm,熔化層兩端凸起的厚度約為2.1μm;入射能流密度若為2500 MW.m?2,靶板表面粗糙度約為3.4μm,熔化層兩端凸起的厚度約為5.4μm;入射能流密度若為2900 MW.m?2,靶板表面粗糙度約為5.8μm,熔化層兩端凸起的厚度約為8.4μm.ELM持續時間越長,熔化層兩端凸起厚度的增加速度越快;考慮熔化層流動時,熔化層邊緣區域的溫度高于不考慮熔化層流動時的溫度,中心區域溫度相對降低,熔化層厚度相對增加.

2)當入射能流密度相同時,表面張力系數對溫度的偏導數kα越大,表面張力梯度越大,偏濾器靶板表面形貌變化越大.當表面張力系數對溫度的偏導數kα相同時,入射能流密度越大,偏濾器靶板表面溫度梯度越大,偏濾器靶板形貌變化越大.

ITER每次放電過程中,發生第一類ELMs的次數不少于103[32],偏濾器靶板表面形貌會發生顯著變化,這將嚴重影響材料的性能,甚至誘發電弧.下一步工作,我們將研究在多個ELMs熱流作用下偏濾器靶板的侵蝕.

[1]Xu W,Wan BN,Xie J K2003Acta Phys.Sin.52 1970(in Chinese)[徐偉,萬寶年,謝紀康 2003物理學報 52 1970]

[2]Coenen J W,Arnoux G,Bazy lev B,Matthews G F,Autricque A,Balboa I,C lever M,Dejarnac R,Coff ey I,Corre Y,Devaux S,Frassinetti L,Gau thier E,Horacek J,Jachmich S,KommM,KnaupM,Krieger K,Marsen S,Meigs A,Mertens P,Pitts R A,Puetterich T,Rack M,StampM,SergienkoG,Tamain P,Thompson V,JET-EFDAContributors 2015Nucl.Fusion55 023010

[3]SergienkoG,Bazy lev B,Huber A,Kreter A,Litnovsky A,Rubel M,Philipps V,Pospieszczyk A,Mertens P,SammU,Schweer B,Schmitz O,Tokar M,The TEXTOR Team2007J.Nucl.Mater.363 96

[4]SergienkoG,Bazylev B,Hirai T2007Phy.Scr.T128 81

[5]Coenen J W,Arnoux G,Bazy lev B,Matthews G F,Jachmich S,Balboa I,C lever M,Dejarnac R,Coff ey I,Corre Y,Devaux S,Frassinetti L,Gau thier E,Horacek J,KnaupM,KommM,Krieger K,Marsen S,Meigs A,Mertens Ph,Pitts R A,Puetterich T,Rack M,StampM,SergienkoG,Tamain P,Thompson V,JET-EFDAContributors 2015J.Nucl.Mater.463 78

[6]Federici G,And rewP,Barabaschi P,Brooks J,Doerner R,Geier A,Herrmann A,Janeschitz G,Krieger K,Kukushkin A,Loarte A,Neu R,Saibene G,Shimada M,Strohmayer G,Sugihara M2003J.Nucl.Mater.313 11

[7]FedericiG,Loarte A,Strohmayer G 2003P lasma Phys.Con trol.Fusion45 1523

[8]Raff ray AR,Federici G 1997J.Nucl.Mater.244 85

[9]Federici G,Raff ray AR 1997J.Nucl.Mater.244 101

[10]Hassanein A,Konkashbaev I2000Fusion Eng.Des.51 681

[11]Sizyuk V,Hassanein A2015Phy.Plasmas22 013301

[12]Litunovsky V N,Kuznetsov V E,Lyub lin BV,Ovchinnikov IB,Titov V A,Hassanein A2000Fusion Eng.Des.49 249

[13]Shi Y,Miloshevsky G,Hassanein A2011Fusion Eng.Des.86 155

[14]Hassanein A,Konkashbaev I2003J.Nucl.Mater.313 664

[15]GencoF,Hassanein A2014Laser Part.Beams32 217

[16]W urz H,Bazy lev B,Landman I,Pestchanyi S,G ross S 2001Fusion Eng.Des.56 397

[17]Bazylev B,Wuerz H2002J.Nucl.Mater.307 69

[18]Coenen J W,Bazylev B,Brezinsek S 2011J.Nucl.Mater.415 S78

[19]Bazylev BN,Janeschitz G,Landman IS,Pestchanyi S E 2005Fusion Eng.Des.75 407

[20]Bazylev BN,Janeschitz G,Landman IS,Loarte A,Pestchanyi S E 2007J.Nucl.Mater.363 1011

[21]Igitkhanov Y,Bazy lev B2014IEEE Trans.P lasma Sci.42 2284

[22]Huang Y,Sun J Z,Sang C F,D ing F,W ang D Z 2014Acta Phys.Sin.63 035204(in Chinese)[黃艷,孫繼忠,桑超峰,丁芳,王德真2014物理學報63 035204]

[23]Huang Y,Sun J Z,Hu W P,Sang C F,W ang D Z 2016Fusion Eng.Des.102 28

[24]Miloshevsky G V,Hassanein 2010Nucl.Fusion50 115005

[25]Loarte A2003P lasma Phys.Con trol.Fusion45 1549

[26]Hassanein A,Sizyuk T,Konkashbaev I2009J.Nucl.Mater.390 777

[27]Jiang C B,Zhang Y L,Ding Z P 2007Computational F luid Mechanics(the fi rst edition)(Beijing: China Power Press)p211(in Chinese)[江春波,張永良,丁則平2007計算流體力學(第一版)(北京:中國電力出版社)第211頁]

[28]CarslawHW,Jaeger J C 1959Conduction of Heat in Solids(2nd Ed.)(Ox ford:C larendon Press)pp89–91

[29]Beh risch R 2010J.Synch.Investig.4 549

[30]Semak V V,Damk roger B,Kempka S 1999J.Phys.D:Appl.Phys.32 1819

[31]Udaykumar HS,Shyy W 1995In t.J.Heat Mass Transfer38 2057

[32]Bazylev BN,Janeschitz G,Landman IS,Pestchanyi S E 2005J.Nucl.Mater.337 766

[33]Udaykumar HS,Shyy W,RaoMM1996In t.J.Numer.Methods Fluids22 691

[34]W urz H,Pestchanyi S,Bazy lev B,Landman I,Kappler F 2001J.Nucl.Mater.290 1138

[35]E lsholz F,Scholl E,Scharfenorth C,Seewald G,Eich ler HJ,Rosen feld A2005J.Appl.Phys.98 103516

[36]E lsholz F,Scholl E,Rosen feld A2004Appl.Phys.Lett.84 4167

PACS:52.55.Fa,52.55.Rk,52.40.Hf,52.65.KjDOI:10.7498/aps.66.035201

N umerical study of thermal erosion and topograph ical change of d ivertor target plates induced by type-Iedge-localized modes?

Huang Yan1)2)Sun Ji-Zhong2)?Sang Chao-Feng2)Hu Wan-Peng2)Wang De-Zhen2)?

1)(School of Information Science and Engineering,Dalian Polytechnic University,Dalian 116034,China)2)(School of Physics and Optoelectronic Engineering,Dalian University of Technology,Dalian 116024,China)(Received 22 September 2016;revised manuscript received 24 October 2016)

The high-Z material tungsten(W)is a promising candidate of the plasma facing components(PFCs)for the future tokamak reactors due toits high melting point(3683 K),lowtritiumretention and lowsputtering yield.However,there are stillmany problems aboutW PFCs.One of themis thematerialmelting under off-normal transient heat fluxes—it is one of themost outstanding open questions associated with the use ofW divertor targets in international thermonuclear experimental reactor(ITER).This requires us urgently tounderstand the W melting behavior under high power flux deposition condition.In this paper,a two-dimensional(2D)fl uid dynamic model is employed by solving the liquid hydrodynamic Navier-Stokes equation together with the 2D heat conduction equation for studying the erosion of the divertor tungsten targets and its resu lting topographicalmodification during a type I-like edge-localized mode(ELM)in ITER with a Gaussian power density profi le heat load.In the present model,major interaction forces,including surface tension,pressure gradient and magnetic force responsible for melt layer motion,are taken intoaccount.The simulation results are fi rst benchmarked with the calcu lated resu lts by other code tovalidate the presentmodeland code.Simulations are carried out in a wide range of fusion plasma performance parameters,and the results indicate that the lifetime ofW plate is determined mainly by the evolution of themelt layer.As a consequence of themelt layermotion,melted tungsten is flushed tothe periphery,a rather deeperosion dent appears,and at the dent edges twohumps of tungsten formduring the ELM.The humps at both edges are almost at the same height.Calculated results showthe topographicalmodification becomes noticeab le when the W plate is exposed toa heat flux of 2000 MW.m?2for 0.8 ms(in the simu lation,the parameter kα= ?α/?Tis taken tobe ?9.0× 10?5N.m?1.K?1,where α is the surface tension coeffi cient and Tis the temperature).The values of the humps are both about 2.1μm,and the surface roughness is about 1.1μm.The longer the duration of the ELM,the more rapid ly the humps rise.Themelt flowmay account for the higher surface temperature at the pool periphery,and for the largermelt thickness.It is found that when the energy flux is under 3000 MW.m?2the surface tension is a major driving force for themotion ofmelt layer.Under the same heat flux,the bigger the kαused in the simu lation,themore severe the surface topography of the target becomes;while at the same kα,the higher the heat flux,themore severe the surface topography of the target becomes.In addition,a modified numericalmethod algorithmfor solving the governing equations is proposed.

edge-localized modes,tungsten divertor plate,melting,flowing

10.7498/aps.66.035201

?國家磁約束聚變科學項目(批準號:2013GB109001,2013GB107003)、國家自然科學基金(批準號:11275042,11575039)和遼寧省基本科研項目(批準號:2016J027)資助的課題.

?通信作者.E-mail:jsun@d lu t.edu.cn

?通信作者.E-mail:wangdez@d lu t.edu.cn

*Project supported by the National Magnetic Con finement Fusion Science Programof China(G rant Nos.2013GB109001,2013GB107003),the National Natural Science Foundation of China(G rant Nos.11275042,11575039),and the Scientifi c Research Foundation of the Liaoning Province,China(G rant No.2016J027).

?Corresponding author.E-mail:jsun@d lut.edu.cn

?Corresponding au thor.E-mail:wangdez@d lut.edu.cn