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

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

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

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

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

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

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

邊界局域模,鎢偏濾器靶板,熔化,流動(dòng)

1 引 言

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

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

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

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

2 物理模型與數(shù)值求解

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

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

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

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

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

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

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

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

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

液體表面的汽化過程導(dǎo)致液體表面不斷向下移動(dòng),模型中引進(jìn)了運(yùn)動(dòng)坐標(biāo)z′:

其中,z=0的位置是靶板表面的初始位置.因?yàn)榘l(fā)生ELMs時(shí),極向能流密度在寬為0.1 m的范圍內(nèi)基本按高斯規(guī)律衰減,打擊點(diǎn)中心的能流密度最大[32],導(dǎo)致作用于鎢熔化層各點(diǎn)的壓強(qiáng)和表面張力不同.因此本文采用二維熱傳導(dǎo)模型研究極向ELM能流不均勻?qū)е骆u熔化層在多種力作用下流動(dòng)時(shí),鎢偏濾器靶板的侵蝕和形貌變化.利用運(yùn)動(dòng)坐標(biāo)(10),熱傳導(dǎo)方程(3)可以寫成如下形式:

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

3 模擬結(jié)果

3.1 模型驗(yàn)證

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

在模型的數(shù)值求解過程中,T,h和ux三個(gè)變量均采用交錯(cuò)網(wǎng)格的方法進(jìn)行離散,優(yōu)點(diǎn)是對速度分量無需再進(jìn)行內(nèi)插運(yùn)算,同時(shí)有效地避免了常規(guī)算法處理液體流動(dòng)這類問題易出現(xiàn)的數(shù)值噪聲[33].由于本文將坐標(biāo)原點(diǎn)設(shè)在靶板的邊緣處,熔化層流動(dòng)時(shí),ELM打擊點(diǎn)兩側(cè)的流動(dòng)速度方向不同,在求解過程中分段采用上風(fēng)格式,成功地追蹤了熔化層流動(dòng)過程中液體表面在各個(gè)時(shí)刻的位置,且所得的熔化層表面較光滑,說明本文采用的計(jì)算方法是成功的.

圖2 鎢靶板表面形貌(其中入射能流密度峰值為1350 MW.m?2,作用時(shí)間為1.8ms)(a)本文模擬結(jié)果;(b)內(nèi)嵌圖,文獻(xiàn)[17]中的模擬結(jié)果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穩(wěn)態(tài)運(yùn)行時(shí)打到靶板的能流密度預(yù)期為5—20 MW.m?2[24,25],如果第一類ELMs間隙間沉積在偏濾器靶板的能流密度取10 MW.m?2,則在ELMs間隙間能流作用下,厚度為1 cm的鎢靶板表面溫度會達(dá)到1000 K左右[7,22]. 本文首先計(jì)算了能流密度峰值為2900 MW.m?2,持續(xù)時(shí)間為0.8 ms的ELM作用在初始溫度為1000 K的鎢偏濾器靶板上,考慮熔化層在表面張力和蒸汽壓強(qiáng)梯度力(在討論的能流密度范圍內(nèi),磁場力作用遠(yuǎn)小于表面張力,所以本文忽略了磁場力)的作用下流動(dòng)時(shí),在熔化層兩端形成的像山峰一樣凸起的高度,偏濾器靶板的表面溫度以及熔化層厚度.計(jì)算過程中kα取?9.0×10?5N.m?1.K?1[17],同時(shí)與不考慮熔化層流動(dòng)時(shí)偏濾器靶板的表面溫度和熔化層厚度進(jìn)行了比較.

圖3給出了ELM過程中3個(gè)不同時(shí)刻靶板的表面形貌.研究表明,在能量沉積過程中鎢靶板汽化很弱,作用于鎢靶板表面的壓強(qiáng)梯度力很小,不能引起鎢液體流動(dòng).鎢靶板熔化之后,熔化層中的液體在表面張力的作用下會不斷地由打擊點(diǎn)區(qū)域向熔化層邊緣區(qū)域流動(dòng),在打擊點(diǎn)區(qū)域形成一個(gè)熔池,邊緣處形成形狀似山峰結(jié)構(gòu)的凸起,計(jì)算時(shí)間越長,在鎢液體中產(chǎn)生的熔池越深,邊緣處形成的山峰高度越高.表1列出了從t=0.6 ms到t=0.95 ms,每過0.05 ms靶板表面凸起的最大值,也列出了該時(shí)刻與0.05 ms之前時(shí)刻熔化層凸起最大值的增加量.從表1數(shù)據(jù)可以看出,ELM能流與靶板作用時(shí)間越長,熔化層兩端凸起厚度的增加速度越快.因?yàn)樵诒疚倪x取的計(jì)算參數(shù)下,靶板表面的汽化量很小,靶板表面形貌變化主要由熔化層流動(dòng)造成.ELM能量沉積時(shí)間越長,靶板熔化厚度越大,熔化層流動(dòng)速度越大,由熔化層流動(dòng)導(dǎo)致熔化層厚度增加或減少的幅值越大.

圖3 (網(wǎng)刊彩色)靶板表面形貌演化,其中,能流密度峰值為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 各時(shí)刻熔化層兩端凸起高度的最大值以及0.05 ms時(shí)間內(nèi)熔化層凸起高度的變化值Tab le 1.Themaximumheight of the lumpon themelt layer edge,and increment of the lumpin 0.05 ms.

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

圖4 (網(wǎng)刊彩色)偏濾器靶板的表面溫度的演化,其中,能流密度峰值為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分別給出了考慮熔化層流動(dòng)和不考慮熔化層流動(dòng),不同時(shí)刻偏濾器靶板的熔化厚度.從圖5可以看出:由于熱量隨著時(shí)間不斷向下傳導(dǎo),計(jì)算時(shí)間越長,靶板的熔化厚度越大;而且,不考慮熔化層流動(dòng)計(jì)算所得的熔化層前沿的位置要小于考慮熔化層流動(dòng)熔化層前沿的位置,與文獻(xiàn)[34]的規(guī)律一致,計(jì)算時(shí)間越長,二者相差越大.在t=0.6 ms,不考慮熔化層流動(dòng)時(shí),熔化層前沿位于z=77.90μm;考慮熔化層流動(dòng)時(shí),熔化層前沿位于z=79.89μm處,二者相差了約1.99μm;到t=0.7 ms時(shí),兩種情況下熔化層前沿的位置相差了約2.24μm;到t=0.8 ms時(shí),此參數(shù)約為2.62μm.

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

圖5 (網(wǎng)刊彩色)靶板熔化厚度的演化,其中,能流密度峰值為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α幅值的影響

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

圖6 (網(wǎng)刊彩色)靶板表面的形貌分布,其中,ELM能流密度峰值為2000 MW.m?2,沉積時(shí)間為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能量范圍預(yù)期為1—3 MJ.m?2,沉積時(shí)間為0.1—1 ms[24,25],本文計(jì)算了能流密度峰值分別為2000,2500,2900 MW.m?2的ELM(如圖7(a)所示)作用在初始溫度為1000 K的偏濾器靶板上,在持續(xù)時(shí)間為0.8 ms,kα取?9.0× 10?5N.m?1.K?1時(shí),偏濾器靶板的表面形貌.計(jì)算結(jié)果表明:入射等離子體能流密度越大,由于鎢液體流動(dòng)在熔化層邊緣處形成的凸起的厚度越大,這是由于入射能流密度越大,能流密度梯度越大,靶板熔化后,靶板表面溫度梯度越大,在kα取相同值時(shí),鎢液體的表面張力梯度越大,鎢熔化層的流動(dòng)速度越快;而且,能流密度越大靶板越先熔化、流動(dòng),所以在鎢熔化層兩端形成的山峰高度越高.

圖6和圖7(b)分別說明kα和ELMs能流密度會影響靶板表面形貌.進(jìn)一步分析,圖6(紅色虛線曲線)和圖7(b)(紅色點(diǎn)劃線曲線)計(jì)算所得的熔化層兩側(cè)山峰的高度都約為5μm,兩幅圖的計(jì)算時(shí)間相同,靶板的初始溫度相同,但圖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倍.可見,由鎢中雜質(zhì)濃度變化等因素導(dǎo)致的鎢熔化層表面張力改變對靶板表面形貌的演化只起次要作用,而ELMs本身的特征參數(shù)變化對靶板表面形貌的改變起主要作用.

圖7 (網(wǎng)刊彩色)能流密度峰值分別為2000,2500和2900 MW.m?2時(shí)的(a)能流密度分布曲線和(b)偏濾器靶板表面形貌;ELM沉積時(shí)間為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 討 論

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

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

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

4 結(jié) 論

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

1)ELM能流與靶板相互作用的過程中,鎢熔化層會發(fā)生流動(dòng),作用時(shí)間為0.8 ms,kα為?9.0 × 10?5N.m?1.K?1, 入射能流密度 大于2000 MW.m?2時(shí),偏濾器靶板表面形貌發(fā)生明顯變化.入射能流密度若為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持續(xù)時(shí)間越長,熔化層兩端凸起厚度的增加速度越快;考慮熔化層流動(dòng)時(shí),熔化層邊緣區(qū)域的溫度高于不考慮熔化層流動(dòng)時(shí)的溫度,中心區(qū)域溫度相對降低,熔化層厚度相對增加.

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

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

[1]Xu W,Wan BN,Xie J K2003Acta Phys.Sin.52 1970(in Chinese)[徐偉,萬寶年,謝紀(jì)康 2003物理學(xué)報(bào) 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物理學(xué)報(bào)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計(jì)算流體力學(xué)(第一版)(北京:中國電力出版社)第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

?國家磁約束聚變科學(xué)項(xiàng)目(批準(zhǔn)號:2013GB109001,2013GB107003)、國家自然科學(xué)基金(批準(zhǔn)號:11275042,11575039)和遼寧省基本科研項(xiàng)目(批準(zhǔn)號: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