1 引言
Ti-6Al-4V合金具有良好的比强度、抗拉强度和耐腐蚀性,广泛应用于航空航天、汽车、化工与海洋工程等领域[1]。丝材电弧增材技术(WAAM)因其成本低、增材效率高、适用广泛和材料利用率高等优点,成为制备大尺寸且形状复杂的Ti-6Al-4V合金部件的重要方法[2]。
对增材制造的数值模拟是展示增材制造机理和提升增材制造水平的重要手段,其中对微观结构和力学性能的数值模拟尤为重要[3]。在模拟微观机构的方法中,蒙特卡洛模型是一种常用的方法,具有计算效率高的特点。Zhang等[4]采用双尺度蒙特卡洛方法模拟了激光增材制造工艺下Ti-6Al-4V凝固过程中β晶粒的形成以及冷却后β晶粒向α晶粒转变的过程。在电弧增材过程中,增材截面上β晶粒以等轴晶和柱状晶的方式生长,晶粒尺寸的变化影响了结构的力学性能,材料的屈服强度随晶粒尺寸增大而减小。Wang等[5]研究了电弧移动方向及层间停留时间对钛合金的晶粒形貌和力学性能的影响。Wu等[6]采用光学显微镜、X射线衍射仪、扫描电镜和标准拉伸实验等试验研究了热累积对电弧增材Ti-6Al-4V的微观结构和力学性能的影响。
晶体塑性力学模型根据材料内部各滑移系开动后的累积剪切应变来模拟材料的塑性变形,可以有效计算模型的塑性变形和力学性能。苑红磊等[7]根据蒙特卡洛方法和晶体塑性模型,研究了搅拌摩擦焊工艺下6005A-T6 铝合金的力学行为。Ti-6Al-4V的力学性能由位错密度和晶粒尺寸共同决定[8,9]。Zhang等[10]研究了搅拌摩擦焊工艺下Ti-6Al-4V的微观组织变化和基于位错演化的力学性能,预测屈服强度与试验吻合良好。Babu等[11]建立了基于位错密度和空位浓度演化的本构模型,研究了Ti-6Al-4V材料在不同应变率及不同温度下的塑性变形。
在增材层材料中,柱状晶具有各向异性的力学性能,而基于位错密度和空位浓度演化的本构模型建立在各向同性的假设上,难以准确模拟柱状晶的力学性能。因此,本文通过蒙特卡洛与晶体塑性模型相结合的方法,计算柱状晶区域不同方向上的应力应变曲线,进而修正基于位错密度和空位浓度演化的位错动力学数值模型来表征柱状晶的各向异性力学行为。
2 理论与方法
2.1 移动热源模型
采用ABAQUS软件建立有限元模型并进行热分析,材料为Ti-6Al-4V。有限元模型尺寸如图1所示,共有10层,扫描速度为600 mm/min,电流为140 A,增材形式为往复堆积。
图1 热传导分析有限元模型
Fig.1 Finite Element Model for heat conduction analysis
使用双椭球热源来模拟电弧增材制造过程中的热输入,其前半部分的热流密度qf由方程(1)描述,后半部分的热流密度qr由方程(2)描述[12]
(1)
(2)
式中QUI表示热源热量(U为电压,I为电流,
为热效率,
=0.8),ff和fr为电弧能量在双椭球热源前后半轴上的分布系数,且ff+fr=2。af,ar,b和c为双椭球热源的形状参数。
2.2 蒙特卡洛数值模型
在蒙特卡洛方法中,将计算域离散成格点,每个格点上分配离散变量来表征该处的微观状态,具有相同离散变量的相邻格点共同构成一个微观晶粒。计算区域采用蒙特卡洛模型中的N1×N2=NMC网格进行离散。格点能量可以表示为
(3)
式中J为常数,m为与该格点相邻的晶格数,δSiSj为克罗内克函数,Si为随机选取的当前格点,Sj为当前格点的相邻格点。
随机选取一个晶格,计算格点能量,随后选取任意相邻格点的取向数进行更新,计算改变后能量的变化ΔE,并通过晶粒重取向准则进行判断是否接受取向数的改变。各向异性的晶粒重取向准则表示为
(4)
式中pl和pv分别为电弧移动方向和沉积方向上取向数改变的概率,取决于该方向上的温度梯度,kB为玻尔兹曼常数,T为温度。
蒙特卡洛步数在时间域离散形式下与温度和时间之间的关系可以表示为
(5)
式中l0为初始晶粒尺寸,λ为格点长度,a1和n为模型系数,K1和n1为晶粒生长系数,D可表示为
(6)
式中Vm为原子摩尔体积,A为模型系数,Z为单位面积内的平均原子数,Na为阿伏伽德罗常数,h为普朗克常量,ΔSf为熔熵,γ=0.5 J/m2为边界能。
2.3 晶体塑性力学模型
在增材过程中,增材区域上部以柱状晶为主,在电弧移动方向和沉积方向上具有不同的力学性能,因此根据蒙特卡洛方法计算得到的多晶体模型,建立晶体塑性力学模型,计算柱状晶不同方向上的力学性能。晶体的总变形梯度可以表示为
F=FeFp
(7)
式中Fe为晶格拉伸和刚性旋转的弹性变形梯度,Fp为晶体内部滑移的塑性变形梯度,与滑移速率的关系可表示为
(8)
式中N为滑移系个数,mw和nw分别为第w滑移系的滑移方向和其法线方向,为第w滑移系的剪切应变率,可表示为
(9)
式中为参考应变率,m为率敏感系数,τω为第ω滑移系的剪切应力,sgn(τω)为符号函数,gω为临界分剪切应力,演化方程为
(10)
式中N为滑移系数目,hωυ为潜硬化模量,可表示为
(11,12)
式中f为硬化比例,hww为自硬化模量,h0为初始硬化模量,τ0为屈服应力,τs为饱和临界应力,γ为所有滑移系的累积剪切应变。
2.4 位错动力学模型
钛合金的力学性能归因于晶粒尺寸、位错之间的长程相互作用和位错穿过短程障碍的运动,可以分为三部分。
(13)
式中第一项表示位错移动通过短程障碍的作用力,第二项表示位错的长程相互作用力,σath为剪切强度。kb为玻尔兹曼常数,p=0.3和q=1.8为比例系数,为参考应变率,
为塑性应变率。Δf0Gb3为克服短程障碍需要的活化能,其中Δf0为温度相关的系数,b为博格斯矢量,G为剪切模量,m为泰勒因子,α为比例系数,
i为位错密度,其演化模型分为硬化和动态回复两部分。θ为拉伸方向与热源移动方向的夹角,范围为0～π/2,kθ=1.6×103 MPa·m1/2,Dθ为该拉伸方向下的晶粒尺寸,可以表示为
(14)
在硬化过程中,固定位错密度的增殖与塑性应变率成正比
(15)
式中Λ为平均自由程,g为晶粒尺寸,s为亚晶粒尺寸。亚晶的形成和演化可表示为
(16)
式中Kc为与温度相关的模型参数,列入表1。等温条件下的晶粒生长模型可表示为
(17)
式中g0为初始晶粒尺寸,K为模型参数。
位错密度的动态回复主要由位错滑移、攀移[13]和晶粒球化[14]三部分构成
(18)
式中Ω为模型系数,cr为材料常数,eq为平衡位错密度,
为球化率,Ψ为校正系数。球化分为动态和静态两部分
Xg=Xd+(1-Xd)Xs
(19)
(20)
(21)
式中Xd和Xs分别为动态球化量和静态球化量,M,kp和B为随温度变化的材料参数。
材料的表观扩散率Dapp可表示为
(22)
式中D1为晶格扩散系数,Dp为沿位错线扩散系数,Dp0为材料参数,Qp为活化能,为位错横截面的原子数,nρ近似等于
为单位面积内的原子数。晶格扩散系数会受到α相和β相体积分数的影响
(23)
式中cv为空位浓度,f为β相的体分比,w为模型系数。Da和Db分别为α相和β相的晶格扩散系数。
假设只有长程作用力对空位的形成有影响,加入温度变化的影响,空位浓度产生和湮灭的模型为
(24)
式中Ω0=1.76×10-29 m3为原子体积,=0.1为空位生成的能量分数。
为了计算流动应力的演化过程,将位错密度和空位浓度相关公式转化为迭代形式
(25)
(26)
3 结果与讨论
温度场与不同增材层的节点温度历程如图2所示,灰色部分代表增材过程中的熔池区域。随着增材层数的增大,热量逐渐累积,第二增材层节点的峰值温度为1844.55 ℃,到第十层节点峰值温度升高至2123.63 ℃。
图2 温度分布及不同增材层沉积时节点温度时程曲线
Fig.2 Temperature distribution and node temperature history during deposition of different additive layers
将得到的截面节点温度历程代入蒙特卡洛数值模型,得到增材层截面的微观结构,如图3所示。增材层底部的晶粒以等轴晶为主。随着增材层逐渐升高,热量逐渐累积,等轴晶的平均晶粒尺寸更大。在增材层上部,其散热以向底部单向散热为主,晶粒沿温度梯度增大的方向生长,以柱状晶为主。
图3 增材层截面微观结构
Fig.3 Cross section microstructure of additive layers
根据蒙特卡洛方法得到的微观结构,重构了晶体塑性多晶体模型,如图4所示,其中,等轴晶部分晶粒取向采用随机取向,柱状晶部分晶粒取向选取〈0 0 1〉取向,以及附近的随机取向,进行晶体塑性计算。
图4 晶体塑性有限元模型
Fig.4 Crystal plasticity finite element model
β相晶粒为体心立方(BCC)结构,具有48个滑移系,取为0.001/s。柱状晶区域的应力应变曲线如图5所示。其中模型沿纵向拉伸时,屈服强度为867.3 MPa。沿垂向拉伸时,屈服强度为788.4 MPa,与试验值平均屈服强度803 MPa[15]对比,误差为1.8%,验证了晶体塑性力学模型的有效性。
图5 不同拉伸方向下的应力应变曲线
Fig.5 Stress-strain curves under different tensile directions
根据修正后的位错动力学数值模型计算得到的不同温度下柱状晶纵向和垂向拉伸的应力应变曲线如图6和图7所示。纵向拉伸条件下常温时的屈服强度为865.2 MPa,600 ℃时屈服应力为418.76 MPa,1000 ℃时屈服应力为9.3 MPa。垂向拉伸条件下常温时的屈服强度为787.8 MPa,600 ℃时屈服应力为379.8 MPa,1000 ℃时屈服应力为8.3 MPa。随着温度升高,屈服应力逐渐降低。在600 ℃～1000 ℃时,应力-应变曲线出现明显软化。
图6 不同温度下纵向拉伸的应力应变曲线
Fig.6 Stress-strain curves under longitudinal tension at different temperatures
图7 不同温度下垂向拉伸的应力应变曲线
Fig.7 Stress-strain curves under vertical tension at different temperatures
4 结 论
本文提出了基于蒙特卡洛模型和晶体塑性模型相结合的微观结构-力学行为一体化计算模型,并与试验对比验证了模型的有效性。修正了位错动力学模型,能够更好地表征具有各向异性特征的柱状晶的力学行为。在常温下,纵向拉伸的屈服强度为865.2 MPa,垂向拉伸的屈服强度为787.8 MPa,具有明显差异。随着增材层的升高,晶粒尺寸逐渐增大,晶粒形貌由等轴晶向柱状晶转变。随着温度升高,屈服强度逐渐降低。在600 ℃～1000 ℃时,应力应变曲线出现明显软化。
文章引用：王艺飞,陈静远,张 昭.基于晶粒形貌的增材制造力学行为数值模拟[J].计算力学学
