韩世伟1,石多奇1,2,杨晓光1,2,苗国磊1

1 北京航空航天大学能源与动力工程学院, 北京 100191
2 先进航空发动机协同创新中心, 北京 100191

摘要

关键词：高循环疲劳;疲劳分散性;极值概率;多晶微结构;晶体塑性

疲劳导致的破坏在航空发动机典型部件的失效方式中占有很高比例, 其中高循环疲劳引起的失效约占24%[1]. 若能更为准确地掌握材料的疲劳分散性, 不仅可以减少事故的发生, 而且能充分发挥潜在的经济效益. 当前疲劳寿命预测方法多数为经验性, 依赖于大量的实验数据, 不仅耗时费钱, 也使对于疲劳分散性的认识仍停留在定性化的阶段.

疲劳裂纹一般在晶粒内部,晶界,缺陷或自由表面萌生, 依赖于裂纹萌生位置邻域微结构的应力应变状态. 在高循环疲劳(high cycle fatigue, HCF)中, 裂纹萌生阶段寿命约占总寿命80%[2]. 在外部宏观载荷确定时, 构件中最大晶粒尺寸[3]及其位置(表面或内部),局部晶粒集合的取向特征以及缺陷[2]等诸多极端(最弱)因素, 共同决定了高循环疲劳的分散性特征. 因此, 从材料微观层次上合理地控制极限(最弱)微结构, 对于减小高循环疲劳分散性,提高疲劳极限具有重要意义. 目前, 以计算手段定量揭示疲劳分散性规律, 并应用于材料微结构设计, 正在得到重视和发展. McDowell课题组[2,4~13]长期致力于研究多尺度力学, 所跨尺度从原子尺度的分子动力学,细观尺度晶体塑性理论到宏观唯象本构, 应用极值概率理论分析了材料微结构对于疲劳分散性影响的系统方法, 并使用匹配函数关联最差疲劳性能对应的微结构, 提出了更高抗疲劳性能的微结构优化设计方法. 牟园伟和陆山[14]考虑了涡轮盘材料微观结构, 结合Tanaka-Mura疲劳裂纹萌生寿命计算理论[15], 预测了裂纹的萌生和小裂纹扩展寿命. Sharaf等[16]应用极值概率理论分析了微结构对结构钢疲劳分散性的影响规律.

本工作发展了考虑微结构的极值概率方法, 包括几何建模,晶体塑性理论,疲劳指示参数和极值概率理论等. 以航空发动机压气机用TC4合金材料[17]为例, 近似构造合金微结构模型, 计算多晶微结构胞元应力应变响应, 从晶粒细观尺度上的塑性变形得到疲劳指示参数, 获得了不同应变和边界条件下疲劳指示参数的统计规律, 通过极值概率方法给出了疲劳萌生寿命的分散性规律.

1 基于微结构的极值概率计算框架

极值概率计算框架考虑了多晶微结构的随机性特征, 计算流程见图1. 具体包括: (1) 应用Voronoi算法[18]构造出随机微结构胞元, 并划分有限元网格, 赋予材料属性; (2) 对微结构胞元模型施加周期边界条件, 模拟晶粒集合处于构件表面或内部状态; (3) 应用晶体塑性本构模型, 计算多晶胞元的应力应变响应; (4) 计算疲劳指示参数, 表征晶粒受到的裂纹萌生驱动力; (5) 应用极值概率理论对胞元的疲劳指示参数极值进行统计分析.

图1微结构相关的极值概率计算框架

Fig.1Computational framework of microstructure-sensitive extreme value probability

1.1 几何建模

微观上多晶材料由大量晶粒构成, 各晶粒具有唯一的晶体取向. 合金中晶粒的尺寸和取向分布特征受到凝固过程和热机械加工方式的影响. Voronoi算法[18]能够模拟金属凝固过程中晶粒的各向同性生长现象, 快速构造出接近等轴晶微结构的几何模型.

Voronoi算法在空间中随机布置核心点, 通过比较空间点到核心点的距离将其划入最近的核心点, 完成空间的划分. 如果对核心点进行特殊布置, 可以划分出具有不同尺寸分布特征的子空间集合.图2a中为随机分布的核心点, 切分立方空间构成的多晶微结构胞元模型. 与复合材料的代表性体积单元(representive volume element, RVE)不同, 金属多晶材料的微观晶粒结构不存在周期对称特征, 通常取平均响应与宏观试件力学响应相近的晶粒集合作为代表性体积单元, 且这样的多晶微结构胞元能够描述多晶微结构的统计学特征[2,11].

图2几何建模

Fig.2Geometry model
(a) Voronoi tessellation (b) mesh generation (c) {0001} pole figure (RD--rolling direction, TD--transverse direction)

Voronoi算法划分的凸子空间域由多个平面封闭构成. 有限元分网方式有2种[18], 分别为完全适应晶界的自由网格(free mesh)和半适应晶界的映射网格(map mesh), 如图2b所示. 自由分网方式适应晶界平面将模型划分为四面体单元, 缺点是容易在多个晶粒的交界位置或者小晶粒中生成畸变的四面体单元, 而且需要很大网格数量. 相比而言, 映射型分网将多晶集合划分为规则一致的六面体单元, 晶界呈锯齿状, 虽然会造成应力,应变过渡不光滑, 但是在计算晶粒及更大尺度平均物理量时, 映射型分网不仅计算效率更高, 而且能够保持相近的数值精度.

常规铸造工艺可以生成具有随机无序特征的取向分布,图2c所示为钛合金{0001}极图, 满足均匀分布特征. 通过随机算法, 可以构造符合该分布特征的Bunge型Euler角.

此外, 铸件经过锻造,轧辊等机械加工后, 晶粒将沿特定方向变形, 并形成择优取向织构. 对等轴晶几何模型沿特定方向缩放, 赋予晶粒特定的Euler角取向参数, 建立经过机械加工的复杂多晶微结构模型.

1.2 表面和内部边界条件

室温下, 对于低缺陷密度的金属多晶材料, 疲劳裂纹萌生于表面的试件一般具有短寿命特征, 而在长寿命情形中, 疲劳裂纹趋向于在试件内部萌生[7]. 为了模拟处于结构内部(interior)的晶粒集合, 在多晶微结构胞元的3个平行面间施加位移控制的周期边界条件(periodic boundary condition, PBC)[2], 如图3a所示. 然而, 与复合材料胞元模型不同, 多晶微结构胞元的相对表面只满足零阶位移周期性, 而不满足应力周期性特征. 因此, 这种位移周期边界约束方式忽略了结构内部晶粒产生的长程高阶弱相互作用, 近似模拟了全空间(full infinite space)体中多晶微结构胞元的载荷状态. 而处于结构表面(surface)的晶粒集合, 只能忽略部分方向上晶粒间产生的长程高阶弱相互作用, 使得多晶微结构胞元载荷状态近似处于半空间(half infinite space)体的表面, 其边界约束方式如图3b所示.

图3b除了在X和Z方向的相对面间施加周期约束外, 而且在沿着Y方向上的中间面y=0和左侧面y=-1之间也施加周期约束, 从而右侧面y=1为自由表面, 并且y∈(0,1]区间的多晶微结构胞元处于半空间体的表面. 映射型分网方式中, 相对面上的节点间沿正交方向一一对应, 位移约束简便. 对于自由分网方式, 相对面上的节点并不满足一一对应关系, 需要根据节点所在投影单元中的坐标进行插值约束.

图3表面和内部对应的周期边界条件

Fig.3Periodic boundary condition
(a) full infinite space (interior) (b) half infinite space (surface)

1.3 晶体塑性

传统的宏观塑性理论定义了唯象的屈服函数和塑性应变演化方程, 不能客观反映晶粒尺度的变形情况. 晶体塑性理论考虑了金属的微观晶体结构, 将塑性变形解释为晶体沿着密排面的密排方向上发生了剪切变形, 更加接近物理事实[19~21]. 根据乘法分解[22,23], 如图1所示, 变形微元体的总体变形梯度F满足:

(1)

式中,Fe和Fp分别为弹性和塑性变形梯度. 而Fp按照滑移产生的顺序, 可以分解为各滑移系上塑性部分变形梯度乘积, 且始终保持晶格方向不变. 因为滑移引起的变形量很小, 忽略滑移顺序的影响,Fp可近似为[24]:

(2)

式中,I为Kronecker张量,γ(α)为第α个滑移系的滑移量,s0(α)和m0(α)为第α个滑移系在参考构形中的单位法向量和切向量. 联立式(1)和(2), 从F中分解出Fe, 对其进行极分解得到弹性变形和刚体转动张量.

在金属材料发生塑性变形时, 经历了复杂的微观演化过程, 如复杂的位错运动,裂纹萌生等. 但是, 在宏观尺度上并不存在直接对应的可观测物理量能准确且全面地描述上述复杂过程. 通常引入唯象的内部状态变量, 从统计平均意义上表征系统的热力学过程. 通常在滑移系上建立率相关的幂函数形式流动方程[2,9,12,13]:

(3)

式中, 为第α个滑移系上的剪应变率, 为参考应变率, 为分解剪应力, 为背应力, 为各向(双向)同性硬化量,K为拖曳应力,m为表征率相关的指数. 塑性变形中位错网络交错缠绕, 并不存在排列规律, 唯象区分为各向(双向)同性和方向性位错网络. 前者对滑移变形产生双向同性阻力, 即双向同性硬化量, 后者则产生方向性的滑移阻力, 即背应力. 背应力会伴随位错攀移,湮灭等行为发生回复作用, 降低内应力势, 体现为内变量演化方程中具有负反馈作用的回复项. 常用的Armstrong-Frederick模型中, 背应力的演化不但随滑移变形发生线性变化, 而且还伴随有动态回复效应[2,9,12,13]:

(4)

式中, 为第α个滑移系背应力对时间的导数,c为线性硬化模量,d为动态回复系数. 式(4)属于渐近型微分方程. 同样, 各向(双向)同性硬化量也表达为线性硬化和动态回复作用的组合[2,9,12,13]:

(5)

式中, 为第α个滑移系各向同性硬化量关于时间的导数,w为硬化量收敛值,q为模量系数. 与背应力演化方程不同, 各向同性硬化量初始值不为0, 并且随着滑移单调递增, 与滑移的方向无关.

1.4 疲劳指示参数

疲劳裂纹的萌生通常伴随着驻留滑移带的形成和晶粒的剪切变形. 基于临界平面的循环塑性剪应变方法已经广泛应用于关联疲劳裂纹的萌生和早期扩展. Fatemi和Socie[25]提出了以剪切应变主导裂纹萌生的疲劳指示参数(fatigue indicator parameter, FIP). Bennett等[26]将FIP应用在晶粒尺度, 得到了基于晶粒的FIP表达式[2]:

(6)

式中, 为循环载荷一个周期内的最大塑性应变增量Δεp在晶粒内进行体积平均后的最大剪应变; 设最大剪应变所在平面为临界平面, 该平面上作用的正应力为σn,max, 正应力σn,max将促进滑移面间的分离以及裂纹的萌生和初始扩展;σp0.2为0.2%塑性应变对应的屈服应力, 对正应力σn,max进行了归一化;k考虑了疲劳萌生过程中正应力的影响, 取值范围为[0, 1]. FIP值越大, 所对应于的裂纹萌生驱动力也越大. Bennett等[26]指出, FIP与疲劳微结构裂纹分布密度存在很好的相关性, 佐证了FIP表征裂纹萌生驱动力的有效性.

1.5 极值统计分布

疲劳裂纹一般在微结构中的薄弱点萌生, 该点附近邻域的裂纹萌生驱动力也最大, 对应于多晶微结构胞元中的FIP极(大)值. 而任一多晶微结构胞元, 在外部载荷作用下, 应用确定性的晶体塑性本构, 得到唯一的应力应变响应以及确定的FIP极值.

随机生成的多晶微结构胞元, 施加相同的外部载荷, 由晶体塑性本构将会得到不同的力学响应, 且每个胞元具有各自的FIP极值, 对应于各胞元的最大裂纹萌生驱动力. 设相同载荷作用下, 全体随机多晶微结构胞元中晶粒FIP构成的样本总体Y满足函数分布Φ(y). 在该样本总体Φ(y)中, 进行样本容量为n且满足一定概率分布规律的随机抽样, 样本Yn的极值xn表示为:

(7)

式中,yi(i=1, 2, , et al.,n)为多晶胞元中第i个晶粒的FIP;Yn表示yi形成的n维向量;xn为含有n个晶粒的多晶胞元中晶粒FIP的最大值, 且xn符合极值概率分布. 常用的Gumbel极值分布函数P为[27]:

(8)

式中, 下标n对应于样本容量;un表示当前样本容量下随机变量的模;βn表征随机变量的分散程度, 与样本容量相关.

FIP极值作为当前胞元受到的最大裂纹萌生驱动力, 其概率分布特征将对应于高循环疲劳的分散性, 两者之间存在着关联性.

2 结果与分析

2.1 多晶胞元及FIP统计分析

TC4合金广泛应用于航空发动机压气机部件, 具有α/β两相微结构特征. 初生α相体积分数约为60%, 且晶粒间相互连接, 平均晶粒尺寸约为15 μm[28], 晶粒取向服从图2c所示的均匀分布. 由于初生α相体积分数很高, 元素偏聚严重, 导致片层组织强度降低. 通过中间退火处理, 使得α相稳定元素(如Al, O等)从初生α相向片层组织α相中扩散, 增加片层组织中稳定元素含量, 提高片层组织中α相强度[21]. 初生α相晶粒相连以及聚集, 位错运动自由程增大, 容易产生塑性滑移, 尤其当其处于自由表面时, 使得初生α相晶粒易成为高循环疲劳的萌生源[2,7].

为了研究与微结构随机性相关的统计规律, 应用Neper软件[18]模拟TC4合金微结构特征, 构造了100个(多晶胞元个数N=100)随机多晶微结构胞元, 划分为图2b所示映射型网格. 每个微结构胞元中包含125个晶粒(样本容量n=125), 能够近似描述合金的宏观力学行为和相关统计特征[2]. 全体多晶胞元中晶粒体积的总体分布如图4中所示, 近似满足Gamma分布. 这100个多晶胞元中全体晶粒(125×100=1.25×104)构成有限容量的样本总体, 体积均值9839.64 μm3, 标准差4816.3 μm3. 相应地, 各多晶微结构胞元的样本容量为125, 样本均值与总体一致. 各多晶胞元构成样本的体积标准差却与总体并不一致, 频数分布如图5所示. 样本体积标准差分布在区间(3500 μm3, 6000 μm3)上, 表明当前应用Voronoi算法划分空间时, 并未考虑子空间(晶粒)体积的标准差. 因此, 在模拟晶粒体积分布影响时, 当前生成的微结构胞元只满足一阶均值近似, 而不满足二阶标准差及更高阶统计特征. 同样地, 赋予所有多晶微结构胞元如图2c所示的均匀分布取向.

图4晶粒体积概率密度分布

Fig.4Probability density distribution of grain volume (k--shape parameter,q--scale parameter)

图5微结构胞元中晶粒体积的标准差频数分布

Fig.5Frequency distribution of standard deviation in grain volumes within each microstructure representative volume element (RVE)

初生α相为hcp结构, 具有基面<>{0001}, 柱面<>{}, 2种锥面<>{}和<>{}共4类包括24个滑移系, 其中基面和柱面滑移系最易开动; 片层组织中hcp的α相与bcc的β相间满足Burgers取向关系, 其中滑移面或者滑移方向与两相界面平行的滑移系具有更低的滑移阻力, 其余滑移系较难开动[10]. 根据前文本构方程, 由实验数据优化出合适的材料参数. 室温下应变率为1×10-3s-1时, TC4合金的拉伸曲线如图6所示. 为了比较不同的应变幅,内部或表面约束方式对多晶微结构胞元中裂纹萌生驱动力的影响, 按照表1中参数进行批量计算.表1中第3列给出在单调拉伸时各应变幅对应的宏观塑性应变增量. 每种情形均计算相同的100个多晶微结构胞元, 应变比为-1, 计算前3个循环.

图12为拟合Gumbel分布得到的相对分散性参数关于应变幅的变化曲线. 对比表2中数据可得, 应变幅值的增大使得Gumbel分布函数的相对分散性参数βn/un下降. 而且, 应变幅小于0.50%时, 分散性参数转变速度最为急剧, 当应变幅超过0.50%的弹性极限时, 分散性参数下降速度趋于平缓. 参见表1或图6可知, 当应变幅超过0.50%时, 多晶微结构胞元发生宏观塑性变形, 向低循环疲劳范畴转变, 呈现低分散性特征. 同时, 因为极值FIP升高, 寿命也将普遍降低.

3 结论

(1) 实现了微结构对多晶合金高循环疲劳分散性影响的可计算性.

(2) 对高循环疲劳, 当应变幅于弹性范围内, 多晶微结构胞元的FIP均值和标准差跨越2~3个量级, 产生了很宽的分散带.

(3) 多晶微结构胞元的FIP极值分散性随应变幅升高而降低, 且在接近弹性极限点附近转变剧烈, 与实验观察到的低应力下疲劳寿命分散性大, 高应力下疲劳寿命分散性小的规律一致.

(4) 位于表面的微结构比内部微结构更加容易萌生裂纹, 且分散性有所下降, 与实验中短寿命情形对应于裂纹表面萌生机制的规律一致.