朱昶勝,李玉杰,馬芳蘭,馮 力,雷 鵬

(1. 蘭州理工大學(xué) 計算機與通信學(xué)院,甘肅 蘭州 730050; 2. 蘭州理工大學(xué) 省部共建有色金屬先進加工與再利用國家重點實驗室,甘肅 蘭州 730050; 3. 蘭州理工大學(xué) 網(wǎng)絡(luò)與信息中心,甘肅 蘭州 730050)

目前,以Fe-C合金為基礎(chǔ)的鋼鐵材料被廣泛應(yīng)用于工業(yè)領(lǐng)域,在整個國民經(jīng)濟中起著舉足輕重的作用.為獲得所需性能,熱處理在加工過程中發(fā)揮著重要作用,與其他加工工藝相比,熱處理不是改變整體化學(xué)成分和工件形狀,而是改變鋼鐵的內(nèi)部結(jié)構(gòu),充分挖掘材料性能潛力,降低結(jié)構(gòu)重組,節(jié)省材料,而且能提高機械產(chǎn)品質(zhì)量、大幅度延長機器零件的使用壽命[1-3].采用微觀組織數(shù)值模擬技術(shù)是預(yù)測和研究熱處理機制和機理的有效手段,已成為研究的熱點.

基于金茲堡-郎道(Ginzburg-Landau)理論基礎(chǔ)的相場方法(phase-field method,PFM)是解決微觀組織演變動力學(xué)問題最具優(yōu)勢的研究方法,已經(jīng)在金屬研究領(lǐng)域得到了廣泛的應(yīng)用[4-5].由液相到固相的共晶生長在相場方面取得了許多重要的研究成果[6-8],共析生長是從共晶生長的理論和方程中派生出來的固相到固相的轉(zhuǎn)變,共晶生長的研究為共析生長的研究奠定了基礎(chǔ).

目前,國內(nèi)外利用相場法對Fe-C合金共析轉(zhuǎn)變的模擬研究已經(jīng)取得了一定進展,Mushongera等[9]建立相場模型來模擬Fe-C-Mn珠光體的形態(tài)轉(zhuǎn)變,探討了珠光體發(fā)散生長機理.Nakajima[10]采用多相場法模擬了共析鋼珠光體的協(xié)同生長,計算了不同過冷條件下的層間距和生長速度,進一步討論了奧氏體中碳分布不均勻性和變形應(yīng)變對應(yīng)力的影響.馮力等[11]利用相場法模擬研究了不同界面能和不同擴散系數(shù)條件下珠光體微觀組織生長形貌及生長規(guī)律.張軍等[12]利用相場法研究了Fe-C合金在等溫過程中發(fā)生奧氏體和鐵素體相變過程,進一步研究了溫度、Mn含量等對微觀組織的影響.

以往的研究主要集中在模擬相變的內(nèi)部生長機制,但是在實際的相變過程中會存在雜質(zhì)或界面影響溶質(zhì)的擴散.因此,研究擋板對相變的影響有重要意義.王華明等[13]模擬研究表明,側(cè)向約束條件對單晶Ni基高溫合金凝固方式及其溶質(zhì)再分配都有著十分顯著的影響.杜立飛等[14]利用相場法模擬了純金屬Ni凝固過程中橫向限制的存在對其枝晶微觀形貌演化的影響.Zhang等[15]利用相場法分析了碳化物的位置和尺寸對珠光體協(xié)同生長的影響.可見,研究阻擋對珠光體生長的影響具有重要價值,以期對實際鑄件提供理論指導(dǎo).

基于傳統(tǒng)求解無法解決計算效率低的瓶頸問題,GPU并行計算已經(jīng)成為解決此類問題的重要手段.目前CUDA是非常成熟的GPU并行計算解決方案,已經(jīng)取得大量研究成果[16-17],但CUDA僅支持Nvidia GPU的產(chǎn)品,不同廠商生產(chǎn)的硬件設(shè)備無法使用統(tǒng)一的編程框架實現(xiàn)程序開發(fā),阻礙了高性能計算的廣泛應(yīng)用和快速發(fā)展.開放式免費標(biāo)準(zhǔn)OpenCL應(yīng)運而生,其能在不同異構(gòu)平臺上執(zhí)行,充分利用GPU硬件資源提高程序性能,在眾多領(lǐng)域擁有廣闊的發(fā)展前景[18-20].

本文采用OpenCL編程框架,充分發(fā)揮GPU極強的浮點運算性能,在不損失計算精度情況下求解多場耦合相場模型,分析了GPU相比CPU的計算效率,研究了不同形狀和不同位置的擋板對層片狀珠光體協(xié)同生長的影響.

1 數(shù)學(xué)模型

1.1 相場模型

采用Kim等提出的相場模型,推導(dǎo)過程詳見文獻[21].

其相場控制方程為

(1)

其中：

(2)

式中：

(3)

式中：在奧氏體相中,φ1=φ2=0,φ3=1；在鐵素體相中,φ2=φ3=0,φ1=1；在滲碳體相中,φ1=φ3=0,φ2=1.

1.2 溶質(zhì)場方程

根據(jù)溶質(zhì)守恒定律得溶質(zhì)場方程為

(4)

2 OpenCL求解原理

圖1為基于OpenCL并行求解二維多場耦合相場模型的原理圖.圖1a為利用OpenCL編程時GPU與CPU的混合框架.其中CPU被視為控制整個計算過程的主機,GPU為并行執(zhí)行多個線程的設(shè)備.CPU執(zhí)行用C語言編寫的主程序,創(chuàng)建并管理CPU與GPU之間通信的數(shù)據(jù)結(jié)構(gòu).內(nèi)核在主機上定義,主機程序發(fā)出一個命令將內(nèi)核代碼通過數(shù)據(jù)總線(PCI-E)提交到GPU設(shè)備上執(zhí)行,與主程序不同的是,內(nèi)核代碼使用了OpenCL標(biāo)準(zhǔn)中所定義的高性能功能.GPU調(diào)用內(nèi)核代碼時,一個內(nèi)核函數(shù)會創(chuàng)建一個對應(yīng)維度的索引空間,即網(wǎng)格.它由多個工作組組成,工作組提供了對網(wǎng)格更粗粒度的分解,網(wǎng)格每個維度能整除各個維度中工作組的數(shù)目.多個工作項構(gòu)成一個工作組.由于僅對二維空間模擬,因此網(wǎng)格和工作組均設(shè)置為二維.圖1b為計算域,采用均勻網(wǎng)格有限差分法計算.圖1a、1b間的箭頭①和圖1a、1c間的箭頭②表示計算結(jié)點和GPU線程之間的對應(yīng)關(guān)系.計算空間中每一個網(wǎng)格點由OpenCL中的一個工作項計算,每個子域由一個工作組計算.計算網(wǎng)格中任意一點處的相場值φi,j和溶質(zhì)場值ci,j分別需要獲得周圍8個網(wǎng)格點的值進行計算.每個工作組邊界上的工作項在計算時與相鄰工作項存在數(shù)據(jù)依賴,需要添加gh1ost cell,即相當(dāng)于在周圍各添加1個額外的網(wǎng)格點,用來存儲邊界工作項需要的數(shù)據(jù),如圖1c中的紅色部分所示.圖1d為基于OpenCL編程在GPU上的線程層次.在硬件結(jié)構(gòu)上,OpenCL設(shè)備由多個計算單元組成,每個計算單元之間功能相互獨立,又可以劃分為一定數(shù)量的處理單元,處理單元是設(shè)備端執(zhí)行數(shù)據(jù)計算的最小單元.工作組由GPU的任務(wù)分配單元分配到計算單元上,工作項會在一個計算單元的處理單元上獨立并發(fā)執(zhí)行.圖1e為OpenCL存儲結(jié)構(gòu).由四部分組成,包括：全局內(nèi)存、私有內(nèi)存、局部內(nèi)存、常量內(nèi)存.任意一個工作項可以讀寫由__global聲明的全局內(nèi)存中內(nèi)存對象的任何元素,局部內(nèi)存供工作組中所有工作項共享,私有內(nèi)存只對單個工作項共享.工作組內(nèi)的各個線程通過計算單元上的局部存儲器進行通信,GPU運算單元訪問局部內(nèi)存的速度非常快.

圖1 OpenCL并行求解多場耦合相場模型示意圖

3 模擬結(jié)果及討論

3.1 初始條件及物理參數(shù)

以Fe-3.4946C合金為例,合金物性參數(shù)及模擬參數(shù)見表1[11].Fe-C合金在低于共析溫度1 000 K下會析出滲碳體和鐵素體組織,二者呈現(xiàn)為交替疊壓的層狀結(jié)構(gòu),也稱為層片狀珠光體.

表1 Fe-C合金物性參數(shù)及模擬參數(shù)

在二維空間求解過程中為了保證數(shù)值收斂,對計算模型在二維空間上的空間步長和時間步長的取值滿足以下關(guān)系：

(5)

式中：D為溶質(zhì)擴散系數(shù)；Δx、Δy分別為在二維空間坐標(biāo)系上平行于x軸、y軸的空間步長參數(shù)；Δt為計算模型的時間步長參數(shù).

一對鐵素體和滲碳體的總寬度稱為珠光體層片間距,以λ表示,其隨過冷度變化的規(guī)律如下式所示：

λ=C/ΔT

(6)

式中：C=8.02×103(nm·K);ΔT為過冷度,即珠光體平衡轉(zhuǎn)變溫度低于共析轉(zhuǎn)變溫度的差值.將共析相變初始溫度設(shè)為980 K,即過冷度為20 K,則初始層片間距為0.4 μm.

圖2為珠光體在二維條件下的初始狀態(tài),γ為奧氏體相、α鐵素體相,β滲碳體相,三相均為固相.S1和S2分別對應(yīng)滲碳體和鐵素體在初始條件下的寬度,兩者之和為最小層片間距λ.使用杠桿定律計算得S1和S2的層片寬度之比約為1∶7,交替設(shè)置鐵素體相和滲碳體相在計算區(qū)域底部,與x軸垂直的方向設(shè)為絕熱邊界條件,與x軸平行的方向設(shè)為周期性邊界條件.

圖2 珠光體初始狀態(tài)Fig.2 The initial state of pearlite

3.2 計算效率分析

實驗平臺在中國科學(xué)院近代物理研究所提供的高性能異構(gòu)平臺上實現(xiàn).該集群操作環(huán)境為Linux OS,共有7個物理節(jié)點,每一個物理節(jié)點有16個可使用的CPU內(nèi)核和3個Tesla M2090 GPU,可使用內(nèi)存為7.24 GB,物理內(nèi)存為31.33 GB,總內(nèi)存為47.07 GB,設(shè)備的內(nèi)存大小為6 GB,最大工作頻率為1.15 GHz,處理器核數(shù)為448,GPU功耗是250 W,內(nèi)存帶寬為177 GB/s,雙精度浮點峰值是515 Gflops.

表2為其他條件相同,計算網(wǎng)格大小不同,分別在GPU和CPU上迭代時間步長為t=1×104Δt所需的計算時間.TC表示CPU串行計算的計算時間,TG表示GPU并行的計算時間.可以看出,對比相同模擬網(wǎng)格大小下CPU和GPU計算的計算時間,GPU更具有計算優(yōu)勢.圖3為GPU相對于CPU的加速比.可以看出,在相同情況下,GPU對不同網(wǎng)格大小的相場模型的加速結(jié)果存在較大差別.隨著網(wǎng)格大小的增加,GPU的加速性能也隨之提升.網(wǎng)格在1 600×1 600計算規(guī)模下,在GPU上計算效率提升兩個數(shù)量級,相比串行程序獲得88倍的加速比.可見,足夠的計算量更能充分發(fā)揮GPU的計算運行效率.

表2 CPU和GPU平臺上的執(zhí)行時間

圖3 GPU對CPU的加速比Fig.3 GPU to CPU speedup ratio

3.3 矩形擋板對珠光體生長形貌的影響

圖4和圖5為存在矩形擋板時Fe-C合金珠光體演化過程的二維模擬結(jié)果.網(wǎng)格規(guī)模為800×800,初始層片間距為100Δx,擋板高度為50Δy,擋板長度為200Δx.其中圖4分別為6×104Δt、12×104Δt、30×104Δt、45×104Δt迭代時間步長的珠光體生長形貌.圖5為對應(yīng)時刻的溶質(zhì)分布,標(biāo)尺p和c分別表示相場參量的值和溶質(zhì)百分?jǐn)?shù).可以看出它們輪廓相同,但溶質(zhì)分布更能清晰地反映珠光體的生長特征.珠光體在未接觸到矩形擋板之前規(guī)則對稱生長.為了更清楚地了解共析層片的生長原理,虛線中的局部放大圖對應(yīng)箭頭所指部分,可以看出,與鐵素體相接觸的奧氏體碳原子濃度較高,與滲碳體相接觸的奧氏體碳原子濃度較低.鐵素體與滲碳體相前沿的奧氏體存在濃度差,從而界面處奧氏體中的碳原子從高濃度向低濃度橫向擴散.碳原子在奧氏體中橫向擴散使得鐵素體界面前沿奧氏體碳原子濃度下降,滲碳體界面前沿奧氏體的碳原子濃度上升,破壞了980 K溫度下三相界面處碳原子濃度的穩(wěn)定.為保持這一穩(wěn)定狀態(tài),與鐵素體接觸的奧氏體析出鐵素體,使其碳原子濃度上升至穩(wěn)定濃度；與滲碳體接觸的奧氏體析出滲碳體,使其碳原子濃度下降至穩(wěn)定濃度.珠光體縱向長大,珠光體呈垂直于固-固界面前沿互相耦合的生長模式.當(dāng)珠光體接觸到矩形擋板時,擋板下方的珠光體縱向生長被完全限制.由于擋板所在的區(qū)域沒有溶質(zhì)分布,珠光體前沿沒有生長所需的碳原子,停止生長.部分珠光體穿過擋板間隙繼續(xù)向上生長,矩形擋板間隙鐵素體的生長速率超過橫向擋板上方鐵素體的生長速率.由于穿過擋板后破壞了初始層片間距的穩(wěn)定性,上方區(qū)域的溶質(zhì)沒有達到滲碳體相重新形核的條件,奧氏體向鐵素體和滲碳體轉(zhuǎn)變的平衡被打破,出現(xiàn)發(fā)散生長,層片間距增大,生長速率慢,如圖4c、圖5c所示.隨著繼續(xù)生長,珠光體層片振幅增大.滲碳體前沿碳原子能快速得到補充,而鐵素體前沿的碳原子不能快速降低,橫向形核和縱向長大平衡被打破,使得層片振幅增大.因此,矩形擋板的存在使得繞過擋板后的珠光體形貌發(fā)生明顯變化,通過擋板間隙的珠光體發(fā)散增長,層片間距增大.

圖4 橫向擋板下珠光體的生長形貌

圖5 橫向擋板下珠光體的溶質(zhì)分布

3.4 十字架擋板對珠光體生長形貌的影響

圖6表示采用不同時刻下十字架擋板對珠光體生長的二維溶質(zhì)場模擬結(jié)果.網(wǎng)格大小為700×1 000.圖6a為t=5×104Δt時刻,在未接觸到擋板之前,珠光體的生長有賴于碳原子從鐵素體前沿富碳奧氏體向滲碳體前沿貧碳奧氏體中擴散,滲碳體前沿碳原子含量的增加有利于滲碳體長大,反之,鐵素體前沿碳原子含量降低促使鐵素體長大,珠光體橫向形核和縱向長大保持平衡,珠光體規(guī)則穩(wěn)定生長.當(dāng)接觸到擋板時,擋板下方的珠光體生長會被限制,穿過擋板間隙的珠光體繼續(xù)縱向生長.圖6b為t=18×104Δt時刻,擋板的存在破壞了珠光體的穩(wěn)態(tài)生長形式,使得滲碳體相前沿碳原子的富集程度不足以讓其規(guī)則向上生長,滲碳體向碳原子濃度高的地方生長,靠近擋板的兩個滲碳體傾斜生長,遠離擋板一側(cè)有傾斜生長的趨勢.圖6c為t=37×104Δt時刻,滲碳體在擋板上形成振幅較大的滲碳體,滲碳體周圍的碳原子濃度不足以支持其繼續(xù)生長,致使?jié)B碳體停止生長.由于珠光體橫向形核和縱向長大平衡進一步被破壞,兩側(cè)滲碳體傾斜加劇,并使相鄰的兩個層狀鐵素體凝固成一個層狀.圖6d為t=45×104Δt時刻,兩側(cè)的珠光體傾斜生長,調(diào)整了層片間距,使珠光體層片間距增大,滲碳體前沿碳原子能快速得到補充,珠光體層片振幅增大,隨著進一步生長,碳原子富集程度達到滲碳體形核需求,重新形成不規(guī)則滲碳體.

圖6 十字架擋板存在時珠光體生長的溶質(zhì)場圖Fig.6 Solute field diagram of pearlite growth in the presence of cross baffle

3.5 不同位置的橫向擋板對珠光體生長形貌的影響

圖7為不同位置相同大小條件下的擋板珠光體的生長形貌.圖8為對應(yīng)的溶質(zhì)分布.模擬規(guī)模為800×800,擋板大小均為15×15.圖7a擋板位于滲碳體相正上方,圖7b位于鐵素體相正上方,圖7c位于滲碳體和鐵素體相正上方.從圖中可以看出生長形貌和溶質(zhì)分布擁有相同的輪廓,但溶質(zhì)分布更能清晰反映珠光體的生長特征.

圖7 擋板在不同位置下珠光體的生長形貌

圖8 擋板在不同位置下珠光體的溶質(zhì)分布

對圖7進行分析可知,當(dāng)擋板處于滲碳體相正上方時,正方形擋板邊長大于滲碳體相寬度,珠光體生長至擋板處會受到限制,使得滲碳體無法繼續(xù)生長,與之相鄰的兩個鐵素體相合并為一個鐵素體相,擋板兩側(cè)局部片層間距較初始層片間距大,隨著珠光體生長,調(diào)整層片間距,當(dāng)處于鐵素體相中心正上方時,緊鄰擋板兩側(cè)的滲碳體相合并為不規(guī)則形狀,滲碳體相前沿碳原子不能滿足其生長需求,停止生長,破壞了初始共析層片的穩(wěn)定性,調(diào)節(jié)層片間距,滲碳體層片震蕩加劇.當(dāng)擋板處于滲碳體和鐵素體相正上方,鐵素體片層會繞過擋板,相鄰的滲碳體會出現(xiàn)彎曲,一段時間后,隨著初始層片間距規(guī)則向上增長.可見,擋板位置的不同對珠光體生長有顯著的影響.

4 結(jié)論

采用OpenCL編程框架,在GPU集群上求解多場耦合相場模型,分析了擋板對珠光體生長的影響.實驗結(jié)果表明：1) GPU相比CPU來說,GPU在計算效率方面更具優(yōu)勢,隨著計算規(guī)模的增加,GPU的加速性能也隨之提高.2) 擋板的存在對珠光體生長形貌有顯著影響,擋板下方的珠光體生長會被限制,穿過擋板間隙的滲碳體傾斜生長.3) 當(dāng)擋板核心處于滲碳體和鐵素體相正上方,相鄰的滲碳體會出現(xiàn)彎曲.

因此,擋板的存在可以作為一種材料制備過程當(dāng)中的微觀結(jié)構(gòu)控制方法,來改變珠光體的生長形態(tài)以獲得具有相對優(yōu)異結(jié)構(gòu)的材料,實現(xiàn)材料設(shè)計.