基于OpenCL并行的挡板对珠光体生长的相场法模拟

朱昶胜,李玉杰,马芳兰,冯 力,雷 鹏

(1. 兰州理工大学 计算机与通信学院,甘肃 兰州 730050; 2. 兰州理工大学 省部共建有色金属先进加工与再利用国家重点实验室,甘肃 兰州 730050; 3. 兰州理工大学 网络与信息中心,甘肃 兰州 730050)

目前,以Fe-C合金为基础的钢铁材料被广泛应用于工业领域,在整个国民经济中起着举足轻重的作用.为获得所需性能,热处理在加工过程中发挥着重要作用,与其他加工工艺相比,热处理不是改变整体化学成分和工件形状,而是改变钢铁的内部结构,充分挖掘材料性能潜力,降低结构重组,节省材料,而且能提高机械产品质量、大幅度延长机器零件的使用寿命[1-3].采用微观组织数值模拟技术是预测和研究热处理机制和机理的有效手段,已成为研究的热点.

基于金兹堡-郎道(Ginzburg-Landau)理论基础的相场方法(phase-field method,PFM)是解决微观组织演变动力学问题最具优势的研究方法,已经在金属研究领域得到了广泛的应用[4-5].由液相到固相的共晶生长在相场方面取得了许多重要的研究成果[6-8],共析生长是从共晶生长的理论和方程中派生出来的固相到固相的转变,共晶生长的研究为共析生长的研究奠定了基础.

目前,国内外利用相场法对Fe-C合金共析转变的模拟研究已经取得了一定进展,Mushongera等[9]建立相场模型来模拟Fe-C-Mn珠光体的形态转变,探讨了珠光体发散生长机理.Nakajima[10]采用多相场法模拟了共析钢珠光体的协同生长,计算了不同过冷条件下的层间距和生长速度,进一步讨论了奥氏体中碳分布不均匀性和变形应变对应力的影响.冯力等[11]利用相场法模拟研究了不同界面能和不同扩散系数条件下珠光体微观组织生长形貌及生长规律.张军等[12]利用相场法研究了Fe-C合金在等温过程中发生奥氏体和铁素体相变过程,进一步研究了温度、Mn含量等对微观组织的影响.

以往的研究主要集中在模拟相变的内部生长机制,但是在实际的相变过程中会存在杂质或界面影响溶质的扩散.因此,研究挡板对相变的影响有重要意义.王华明等[13]模拟研究表明,侧向约束条件对单晶Ni基高温合金凝固方式及其溶质再分配都有着十分显著的影响.杜立飞等[14]利用相场法模拟了纯金属Ni凝固过程中横向限制的存在对其枝晶微观形貌演化的影响.Zhang等[15]利用相场法分析了碳化物的位置和尺寸对珠光体协同生长的影响.可见,研究阻挡对珠光体生长的影响具有重要价值,以期对实际铸件提供理论指导.

基于传统求解无法解决计算效率低的瓶颈问题,GPU并行计算已经成为解决此类问题的重要手段.目前CUDA是非常成熟的GPU并行计算解决方案,已经取得大量研究成果[16-17],但CUDA仅支持Nvidia GPU的产品,不同厂商生产的硬件设备无法使用统一的编程框架实现程序开发,阻碍了高性能计算的广泛应用和快速发展.开放式免费标准OpenCL应运而生,其能在不同异构平台上执行,充分利用GPU硬件资源提高程序性能,在众多领域拥有广阔的发展前景[18-20].

本文采用OpenCL编程框架,充分发挥GPU极强的浮点运算性能,在不损失计算精度情况下求解多场耦合相场模型,分析了GPU相比CPU的计算效率,研究了不同形状和不同位置的挡板对层片状珠光体协同生长的影响.

1 数学模型

1.1 相场模型

采用Kim等提出的相场模型,推导过程详见文献[21].

其相场控制方程为

(1)

其中：

(2)

式中：

(3)

式中：在奥氏体相中,φ1=φ2=0,φ3=1；在铁素体相中,φ2=φ3=0,φ1=1；在渗碳体相中,φ1=φ3=0,φ2=1.

1.2 溶质场方程

根据溶质守恒定律得溶质场方程为

(4)

2 OpenCL求解原理

图1为基于OpenCL并行求解二维多场耦合相场模型的原理图.图1a为利用OpenCL编程时GPU与CPU的混合框架.其中CPU被视为控制整个计算过程的主机,GPU为并行执行多个线程的设备.CPU执行用C语言编写的主程序,创建并管理CPU与GPU之间通信的数据结构.内核在主机上定义,主机程序发出一个命令将内核代码通过数据总线(PCI-E)提交到GPU设备上执行,与主程序不同的是,内核代码使用了OpenCL标准中所定义的高性能功能.GPU调用内核代码时,一个内核函数会创建一个对应维度的索引空间,即网格.它由多个工作组组成,工作组提供了对网格更粗粒度的分解,网格每个维度能整除各个维度中工作组的数目.多个工作项构成一个工作组.由于仅对二维空间模拟,因此网格和工作组均设置为二维.图1b为计算域,采用均匀网格有限差分法计算.图1a、1b间的箭头①和图1a、1c间的箭头②表示计算结点和GPU线程之间的对应关系.计算空间中每一个网格点由OpenCL中的一个工作项计算,每个子域由一个工作组计算.计算网格中任意一点处的相场值φi,j和溶质场值ci,j分别需要获得周围8个网格点的值进行计算.每个工作组边界上的工作项在计算时与相邻工作项存在数据依赖,需要添加gh1ost cell,即相当于在周围各添加1个额外的网格点,用来存储边界工作项需要的数据,如图1c中的红色部分所示.图1d为基于OpenCL编程在GPU上的线程层次.在硬件结构上,OpenCL设备由多个计算单元组成,每个计算单元之间功能相互独立,又可以划分为一定数量的处理单元,处理单元是设备端执行数据计算的最小单元.工作组由GPU的任务分配单元分配到计算单元上,工作项会在一个计算单元的处理单元上独立并发执行.图1e为OpenCL存储结构.由四部分组成,包括：全局内存、私有内存、局部内存、常量内存.任意一个工作项可以读写由__global声明的全局内存中内存对象的任何元素,局部内存供工作组中所有工作项共享,私有内存只对单个工作项共享.工作组内的各个线程通过计算单元上的局部存储器进行通信,GPU运算单元访问局部内存的速度非常快.

图1 OpenCL并行求解多场耦合相场模型示意图

3 模拟结果及讨论

3.1 初始条件及物理参数

以Fe-3.4946C合金为例,合金物性参数及模拟参数见表1[11].Fe-C合金在低于共析温度1 000 K下会析出渗碳体和铁素体组织,二者呈现为交替叠压的层状结构,也称为层片状珠光体.

表1 Fe-C合金物性参数及模拟参数

在二维空间求解过程中为了保证数值收敛,对计算模型在二维空间上的空间步长和时间步长的取值满足以下关系：

(5)

式中：D为溶质扩散系数；Δx、Δy分别为在二维空间坐标系上平行于x轴、y轴的空间步长参数；Δt为计算模型的时间步长参数.

一对铁素体和渗碳体的总宽度称为珠光体层片间距,以λ表示,其随过冷度变化的规律如下式所示：

λ=C/ΔT

(6)

式中：C=8.02×103(nm·K);ΔT为过冷度,即珠光体平衡转变温度低于共析转变温度的差值.将共析相变初始温度设为980 K,即过冷度为20 K,则初始层片间距为0.4 μm.

图2为珠光体在二维条件下的初始状态,γ为奥氏体相、α铁素体相,β渗碳体相,三相均为固相.S1和S2分别对应渗碳体和铁素体在初始条件下的宽度,两者之和为最小层片间距λ.使用杠杆定律计算得S1和S2的层片宽度之比约为1∶7,交替设置铁素体相和渗碳体相在计算区域底部,与x轴垂直的方向设为绝热边界条件,与x轴平行的方向设为周期性边界条件.

图2 珠光体初始状态Fig.2 The initial state of pearlite

3.2 计算效率分析

实验平台在中国科学院近代物理研究所提供的高性能异构平台上实现.该集群操作环境为Linux OS,共有7个物理节点,每一个物理节点有16个可使用的CPU内核和3个Tesla M2090 GPU,可使用内存为7.24 GB,物理内存为31.33 GB,总内存为47.07 GB,设备的内存大小为6 GB,最大工作频率为1.15 GHz,处理器核数为448,GPU功耗是250 W,内存带宽为177 GB/s,双精度浮点峰值是515 Gflops.

表2为其他条件相同,计算网格大小不同,分别在GPU和CPU上迭代时间步长为t=1×104Δt所需的计算时间.TC表示CPU串行计算的计算时间,TG表示GPU并行的计算时间.可以看出,对比相同模拟网格大小下CPU和GPU计算的计算时间,GPU更具有计算优势.图3为GPU相对于CPU的加速比.可以看出,在相同情况下,GPU对不同网格大小的相场模型的加速结果存在较大差别.随着网格大小的增加,GPU的加速性能也随之提升.网格在1 600×1 600计算规模下,在GPU上计算效率提升两个数量级,相比串行程序获得88倍的加速比.可见,足够的计算量更能充分发挥GPU的计算运行效率.

表2 CPU和GPU平台上的执行时间

图3 GPU对CPU的加速比Fig.3 GPU to CPU speedup ratio

3.3 矩形挡板对珠光体生长形貌的影响

图4和图5为存在矩形挡板时Fe-C合金珠光体演化过程的二维模拟结果.网格规模为800×800,初始层片间距为100Δx,挡板高度为50Δy,挡板长度为200Δx.其中图4分别为6×104Δt、12×104Δt、30×104Δt、45×104Δt迭代时间步长的珠光体生长形貌.图5为对应时刻的溶质分布,标尺p和c分别表示相场参量的值和溶质百分数.可以看出它们轮廓相同,但溶质分布更能清晰地反映珠光体的生长特征.珠光体在未接触到矩形挡板之前规则对称生长.为了更清楚地了解共析层片的生长原理,虚线中的局部放大图对应箭头所指部分,可以看出,与铁素体相接触的奥氏体碳原子浓度较高,与渗碳体相接触的奥氏体碳原子浓度较低.铁素体与渗碳体相前沿的奥氏体存在浓度差,从而界面处奥氏体中的碳原子从高浓度向低浓度横向扩散.碳原子在奥氏体中横向扩散使得铁素体界面前沿奥氏体碳原子浓度下降,渗碳体界面前沿奥氏体的碳原子浓度上升,破坏了980 K温度下三相界面处碳原子浓度的稳定.为保持这一稳定状态,与铁素体接触的奥氏体析出铁素体,使其碳原子浓度上升至稳定浓度；与渗碳体接触的奥氏体析出渗碳体,使其碳原子浓度下降至稳定浓度.珠光体纵向长大,珠光体呈垂直于固-固界面前沿互相耦合的生长模式.当珠光体接触到矩形挡板时,挡板下方的珠光体纵向生长被完全限制.由于挡板所在的区域没有溶质分布,珠光体前沿没有生长所需的碳原子,停止生长.部分珠光体穿过挡板间隙继续向上生长,矩形挡板间隙铁素体的生长速率超过横向挡板上方铁素体的生长速率.由于穿过挡板后破坏了初始层片间距的稳定性,上方区域的溶质没有达到渗碳体相重新形核的条件,奥氏体向铁素体和渗碳体转变的平衡被打破,出现发散生长,层片间距增大,生长速率慢,如图4c、图5c所示.随着继续生长,珠光体层片振幅增大.渗碳体前沿碳原子能快速得到补充,而铁素体前沿的碳原子不能快速降低,横向形核和纵向长大平衡被打破,使得层片振幅增大.因此,矩形挡板的存在使得绕过挡板后的珠光体形貌发生明显变化,通过挡板间隙的珠光体发散增长,层片间距增大.

图4 横向挡板下珠光体的生长形貌

图5 横向挡板下珠光体的溶质分布

3.4 十字架挡板对珠光体生长形貌的影响

图6表示采用不同时刻下十字架挡板对珠光体生长的二维溶质场模拟结果.网格大小为700×1 000.图6a为t=5×104Δt时刻,在未接触到挡板之前,珠光体的生长有赖于碳原子从铁素体前沿富碳奥氏体向渗碳体前沿贫碳奥氏体中扩散,渗碳体前沿碳原子含量的增加有利于渗碳体长大,反之,铁素体前沿碳原子含量降低促使铁素体长大,珠光体横向形核和纵向长大保持平衡,珠光体规则稳定生长.当接触到挡板时,挡板下方的珠光体生长会被限制,穿过挡板间隙的珠光体继续纵向生长.图6b为t=18×104Δt时刻,挡板的存在破坏了珠光体的稳态生长形式,使得渗碳体相前沿碳原子的富集程度不足以让其规则向上生长,渗碳体向碳原子浓度高的地方生长,靠近挡板的两个渗碳体倾斜生长,远离挡板一侧有倾斜生长的趋势.图6c为t=37×104Δt时刻,渗碳体在挡板上形成振幅较大的渗碳体,渗碳体周围的碳原子浓度不足以支持其继续生长,致使渗碳体停止生长.由于珠光体横向形核和纵向长大平衡进一步被破坏,两侧渗碳体倾斜加剧,并使相邻的两个层状铁素体凝固成一个层状.图6d为t=45×104Δt时刻,两侧的珠光体倾斜生长,调整了层片间距,使珠光体层片间距增大,渗碳体前沿碳原子能快速得到补充,珠光体层片振幅增大,随着进一步生长,碳原子富集程度达到渗碳体形核需求,重新形成不规则渗碳体.

图6 十字架挡板存在时珠光体生长的溶质场图Fig.6 Solute field diagram of pearlite growth in the presence of cross baffle

3.5 不同位置的横向挡板对珠光体生长形貌的影响

图7为不同位置相同大小条件下的挡板珠光体的生长形貌.图8为对应的溶质分布.模拟规模为800×800,挡板大小均为15×15.图7a挡板位于渗碳体相正上方,图7b位于铁素体相正上方,图7c位于渗碳体和铁素体相正上方.从图中可以看出生长形貌和溶质分布拥有相同的轮廓,但溶质分布更能清晰反映珠光体的生长特征.

图7 挡板在不同位置下珠光体的生长形貌

图8 挡板在不同位置下珠光体的溶质分布

对图7进行分析可知,当挡板处于渗碳体相正上方时,正方形挡板边长大于渗碳体相宽度,珠光体生长至挡板处会受到限制,使得渗碳体无法继续生长,与之相邻的两个铁素体相合并为一个铁素体相,挡板两侧局部片层间距较初始层片间距大,随着珠光体生长,调整层片间距,当处于铁素体相中心正上方时,紧邻挡板两侧的渗碳体相合并为不规则形状,渗碳体相前沿碳原子不能满足其生长需求,停止生长,破坏了初始共析层片的稳定性,调节层片间距,渗碳体层片震荡加剧.当挡板处于渗碳体和铁素体相正上方,铁素体片层会绕过挡板,相邻的渗碳体会出现弯曲,一段时间后,随着初始层片间距规则向上增长.可见,挡板位置的不同对珠光体生长有显著的影响.

4 结论

采用OpenCL编程框架,在GPU集群上求解多场耦合相场模型,分析了挡板对珠光体生长的影响.实验结果表明：1) GPU相比CPU来说,GPU在计算效率方面更具优势,随着计算规模的增加,GPU的加速性能也随之提高.2) 挡板的存在对珠光体生长形貌有显著影响,挡板下方的珠光体生长会被限制,穿过挡板间隙的渗碳体倾斜生长.3) 当挡板核心处于渗碳体和铁素体相正上方,相邻的渗碳体会出现弯曲.

因此,挡板的存在可以作为一种材料制备过程当中的微观结构控制方法,来改变珠光体的生长形态以获得具有相对优异结构的材料,实现材料设计.