单晶涡轮叶片定向凝固过程应力场数值模拟

李世峰， 何爱杰， 邱 飞

(1.中国航发四川燃气涡轮研究院,成都 610500; 2.西北工业大学 机电学院,西安 710072)

单晶涡轮叶片是由大量自由曲面和复杂内腔组成的空心结构薄壁件，是现代高推重比航空发动机的核心零件。涡轮叶片长期工作在环境最恶劣的条件下，所以其性能水平，成为一种航空发动机先进程度的重要标志。因此，为提高其性能，当代航空发动机的涡轮叶片设计普遍采用复合气膜冷却式的单晶空心结构，这类叶片的结构形状极复杂、制造技术难度大，气冷叶片壁厚尺寸控制成为研制航空发动机的关键技术。

近年来，国外针对单晶涡轮叶片壁厚尺寸精度控制问题，进行了大量有关叶片定向凝固过程应力场仿真研究，并在控制叶片尺寸精度、优化工艺参数方面起到较大的作用[1-2]。Galantucci等[2]采用仿真的方法模拟了涡轮叶片的凝固过程，并与实验结果进行了对比。Zhu等[3]采用规则与非规则网格混合的方法，着重解决了辐射边界条件的问题。Onate等[4]采用热-力结构耦合的有限元模型对曲轴铸铁件进行了应力分析。美国GE、普·惠、Howmet和PCC以及英国罗·罗、法国斯奈克玛等公司均采用热-力耦合模型仿真技术[5-6]， 分析和优化了单晶涡轮叶片的生产工艺。目前，国内在叶片设计与分析集成技术、仿真技术等方面都处于起步阶段，尤其在单晶气冷叶片的凝固过程数值模拟方面研究较晚、进展较慢，与国外相比差距最大[7-8]；并且，以往的研究都停留在实心涡轮叶片上，而对结构复杂的单晶气冷叶片的定向凝固过程数值模拟研究较少。

本研究针对单晶空心涡轮叶片壁厚尺寸精度偏低和壁厚尺寸漂移大等问题，采用热力双向耦合的凝固应力场数值模拟技术, 得到其冷却过程温度/应力动态变化及应力/变形情况分布规律，以期为避免叶片出现大铸造残余应力区、预防变形工艺、提高叶片壁厚尺寸精度，以及尺寸稳定性，提供量化参考依据。

1 定向凝固过程分析

单晶叶片采用无余量Bridgman定向凝固技术而成，具有优良的抗热冲击性能、较好的蠕变抗力和较长的疲劳寿命，特别是具有较高的承温能力。应用这种技术能使涡轮叶片的使用温度提高10～30℃，涡轮进口温度提高20～60℃，对提高发动机的推力具有重要的实际意义。

在单晶气冷叶片的定向凝固过程中，将单晶高温合金熔化后浇注到固定型壳内(1550±5) ℃，待熔体温度均匀后，经由抽拉单元体带动，型壳随水冷铜盘结晶器一起自上而下移动。当高温金属液浇注到型腔并接触到铜盘时，由于激冷作用(激冷盘水温(10±0.5)℃)而结晶，并迅速生长进入螺旋选晶器，再经螺旋选晶器作用，在其顶端生长出一个晶粒，并最终充满整个叶片。在此过程中，当叶片模组经辐射挡板进入冷却区后，通过叶片模组与冷却区炉壁间的辐射散热作用，在叶片(沿叶身 [001]方向)液固相区建立起特定方向的温度梯度，从而使熔体沿着与热流方向相反的方向凝固，最终获得具有特定取向(晶体取向[001]与主应力轴夹角在18°范围内)的单晶叶片，其凝固过程示意图如图1所示[9]。

图1 单晶叶片定向凝固工艺示意图Fig.1 Unidirectional solidification process of single crystal blade

2 物理模型的建立

2.1 有限元模型

以先进航空发动机涡轮气冷叶片为对象，根据实际铸造工艺，建立了浇注装配模型，如图2所示。采用某软件进行CAD建模后，采用Hyperworks的Hypermesh模块对叶片浇注装备系统进行有限元网格剖分。考虑到涡轮叶片结构的复杂性，有限元网格划分采用稀疏网格技术，在叶片进、排气边及叶根部位采用长度为0.2 mm的网格，叶片叶身部分采用长度为0.5 mm的网格，浇口杯部分采用1.0 mm的网格，水冷结晶器与炉体采用长度为2 mm的网格。并针对定向凝固模型不同结构热物性参数的要求，将整个浇注装配系统分为铸件、壳型与陶芯、激冷盘三部分，并分别进行有限元建模。

图2 叶片定向凝固浇注系统(a)与有限元模型(b)Fig.2 Casting system(a) and finite element model(b) of directional solidification

2.2 材料热物性参数确定

实验涡轮气冷叶片选用镍基单晶高温合金DD6材料，其熔化温度1340～1380 ℃，液相线温度1370 ℃，固相线温度1320 ℃。该铸件凝固过程中的热物性参数主要包括铸件材质和模壳型材的比热容C，热导率λ和密度ρ=8.78 kg/cm3，这些热物性值均随温度的变化而变化，具体如图3、表1～表3所示。

2.3 边界条件处理

涡轮气冷叶片属小尺寸薄壁铸件，其凝固过程中铸件主要通过铸型向环境散热，因此铸件/铸型间界面换热是整个传热系统的主要方式。铸型在高温区预热至初始浇注温度，与铸件之间存在的热交换微量且可忽略，但当铸件与铸型一起抽拉至低温区后，铸型与低温区炉壁间存在强烈的热辐射。在铸造应力场的计算中，温度场是作为载荷加载到有限元网格节点上，所以，尽可能准确求出温度场是准确预测铸件热应力的保证。在叶片无余量精铸的数值模拟中，铸件与铸型的换热系数很难确定。本研究为精确确定各界面间的换热系数，设计了一个温度场的测量实验，来获得叶片精铸凝固过程温度场数据，然后再利用这些数据，借助大型商用ProCAST软件，对界面换热系数进行反向求解。

图3 DD6合金[001]取向热物性参数 (a)比热容；(b)热传导率Fig.3 Thermal physical parameters of DD6 alloy for [001] orientation (a)specific heat capacity;(b)thermal conductivity

300℃20-400℃20-500℃20-600℃20-700℃20-800℃20-900℃20-1000℃20-1100℃20-1200℃11.9212.5912.9313.1513.5314.1914.3915.0015.7616.75

表2 DD6合金[001]方向不同温度区间的抗拉强度(σb/MPa)

表3 模壳热物理性能参数

2.3.1 温度场实验

温度场实验采用SWP-TSR116型温度记录仪，待热电偶与叶片连接后，通过补偿导线与SWP-TSR116型温度记录仪相连接。共使用10组热电偶，分别测量叶片铸件叶背、叶盆以及浇注系统螺旋选晶器表面处温度-时间变化曲线，如图4所示。

本研究基于有限元软件ProCAST的反求模块[10-12]，以实验实测温度场数据为初始值(见图5)，利用ProCAST的逆运算模块，采用实测温度值，对叶片精铸过程中不同界面间的换热系数进行反向求解，进而获得不同界面的换热系数，具体如图6所示。

图4 测温点分布示意图Fig.4 Distribution of temperature measuring points

图7给出了在热电偶1#处采用铸件/模壳和铸件/水冷结晶器界面换热系数计算的温度与实测值温度的比较。由图7可见，计算温度曲线与实测温度曲线基本吻合，最大误差不超过5 ℃。

2.3.2 边界条件

先将叶片浇注模组装夹在真空炉冷铜盘上，再将真空炉抽真空达到一定真空度(约为0.05 Pa )；待浇注模组在高温区后加热至1550 ℃时，以1410 ℃的温度进行浇注，并保持约10 min后，然后以4.0 mm/min 的速率下降至真空炉低温区，进而完成了单晶叶片的无余量精密成形。

图5 单晶叶片浇注系统实物图Fig.5 Casting system of single crystal blade

3 应力场仿真结果与分析

铸造残余应力是铸件凝固冷却过程中，因发生不均匀膨胀或收缩而引起的应力集中现象，而这种变化往往受到铸件的结构形状与铸造工艺等因素的制约，使其冷却收缩不能自由进行，于是在发生变形的同时产生了内应力。所以，能够准确评估与预测特定结构铸造残余应力的分布形态及其动态特性，将成为控制其铸造残余应力水平的关键[13-14]。为此，本工作基于单晶涡轮转子叶片精铸仿真CAE模型，通过实验测试建立了叶片凝固过程温度-时间边界条件，采用热-力耦合的热弹塑性模型，利用大型商用ProCAST软件，对单晶叶片凝固后阶段的应力场进行仿真计算。

图6 换热系数随时间变化的曲线 (a)铸件/模壳间;(b)铸件/冷铜间Fig.6 Curves of heat transfer coefficient with time (a)casting / mold shell;(b)casting/cold copper

图7 1#热电偶温度实测值与仿真值的对比Fig.7 Temperatures of measured value and simulation value for 1# thermocouple

3.1 叶身应力场分布

选取冷却结构极复杂，气动外形、壁厚尺寸精度要求极高的气冷涡轮转子叶片，其叶身由大量自由曲面和复杂内腔组成，现以叶片的叶身Ⅴ和Ⅷ截面作为分析对象，其结构及截面尺寸如图8所示。

图8 叶片截面尺寸Fig.8 Blade section size

在凝固时间5950 s/7000 step时，其叶身Ⅴ和Ⅷ截面的应力场分布情况如图9所示。对最大应力点的残余热应力进行定量分析，其对应的应力与温度之间的关系如图10所示。

可以看出，叶身应力最大点出现在进气边圆弧中点，这是因为进气边圆弧曲率变化较大，散热较慢，导致其冷却收缩减慢，且在凝固后阶段又受到邻近圆弧面的约束，使其冷却收缩不能自由进行，于是产生了较大的内应力。另外，Ⅴ和Ⅷ截面的叶盆、叶背中部和靠尾缘区均出现大应力，由于受到强度、冷却设计的约束，叶片的叶盆、叶背中部截面厚度较大，横向冷却肋密集，导致热应力集中；而在排气边区域沿径向存在大量扰流柱，且受到冷却设计的限制，该扰流柱尺寸较小，且与叶盆叶背交界处倒角较小，从而极易引起该区域铸造残余应力过大。

图9 叶身截面XY面热应力分布 (a)Ⅷ截面；(b)Ⅴ截面Fig.9 Thermal stress distribution of XY surface with leaf body section (a)section Ⅷ;(b)section Ⅴ

图10 Ⅴ和Ⅷ截面的最大热应力与温度之间的关系Fig.10 Relationship between maximum thermal stress and temperature on section Ⅴ and Ⅷ

图11 榫头进气窗口区域XY截面热应力分布Fig.11 Thermal stress distribution of XY section for the tenon inlet window area

3.2 榫头应力场分布

选取榫头进气窗口倒圆角为R=0.5 mm时的截面作为分析对象，其结构及截面尺寸如图11所示。在凝固时间5950 s/7000 step时，叶片榫头进气窗口区域的热应力沿径向分布情况如图12所示。

榫头进气窗口不同截面最大应力点Pi1，Pi2，Pi3和Pi4处的残余热应力与温度间的关系曲线如图13所示。

图12 榫头进气窗口区域XY截面热应力分布云图(7000 step) (a)H1截面；(b)H2截面Fig.12 Thermal stress distribution of XY section for tenon inlet window (a)section H1;(b)section H2

图13 H1和H2截面不同位置处的残余热 应力与温度关系曲线Fig.13 Relationship between residual thermal stress and temperature at different locations of H1 and H2

从图13可以看出，榫头进气窗口倒圆角R=0.5 mm时，进气窗口锐角区域形成的残余应力最大，其最大残余热应力出现在PH1,1点，而残余应力沿叶片径向变化较小，且随着冷却的进行，残余应力数值不断增加，分布区域也不断扩展。主要是因为：锐角区域与钝角区域结构形状差别较大，从而导致非同步冷却收缩，但铸件又是一个整体，再加陶瓷型芯约束与铸件彼此间互相制约，使收缩受到一定阻碍，于是在铸件内产生铸造热应力；同时，发现先冷却的钝角区域形成压应力，而后冷却的锐角区域形成较大的拉应力。另外，由于榫头进气窗口区域结构壁厚尺寸大，残余热应力很难通过弹性变形释放，于是残余应变能累积起来；当残余应变能累积到一定程度时，残余应力将在铸件内部分布极不均匀，致使叶片表面形成一定的应力梯度。

由此看来，整个叶片最大铸造热应力出现在榫头底面四大进气窗口区域，该区域最大应力比叶身截面最大应力高28.4%，主要因为叶身壁厚较薄，凝固前期形成的热应力均以弹性变形方式释放，而榫头进气窗口区域结构壁厚尺寸大，残余热应力很难通过弹性变形释放，反而累积起来，最终导致叶片榫头进气孔区域的铸造热应力最大。

4 结论

(1)叶片最大铸造热应力出现在榫头底面四大进气窗口区域，该区域最大应力比叶身截面最大应力高28.4%。

(2)当榫头进气窗口倒圆角R=0.5 mm时，进气窗口锐角区域形成的残余应力较大。同时，发现先冷却的钝角区域形成压应力，而后冷却的锐角区域形成较大的拉应力。

(3)采用温度场与应力场耦合的仿真方法，通过实际温度场数据采集，可准确反映单晶涡轮叶片凝固过程应力动态变化，将为提高叶片壁厚尺寸精度及尺寸稳定性，提供了量化参考依据。

[1] MODUKURU S C, RAMAKRISHNAN N，SRIRAMAMURTHY A M. Determination of the die profile for the investment casting of aerofoil-shaped turbine blades using the finite-element method [J]. Journal of Materials Processing Technology, 1996,58(2/3): 223-226.

[2] GALANTUCCI L M, TRICARICO L. A computer-aided approach for the simulation of the directional solidification process for gas turbine blades [J]. Journal of Materials Processing Technology, 1998,77: 160- 165.

[3] ZHU J D， OHNAKA I．Three dimensional computer simulation on mold filling of casting by direct finite difference method[J]．J Jpn Foundry Eng 1996，68：668-676．

[4] 廖敦明.基于有限差分法的铸件凝固过程热应力场数值模拟的研究[D].武汉：华中科技大学,2002.

(LIAO D M. Study on numerical simulation of thermal stresses during casting’s solidification process based on FDM[D]. Wuhan：Huazhong University of Science & Technology,2002 .)

[5] MODUKURU S C , RAMAKRISHNAN N，SRIRAMAMURTHY A M. Determination of the die profile for the investment casting of aerofoil-shaped turbine blades using the finite-element method [J]. Journal of Materials Processing Technology, 1996, 58(2/3): 223-226.

[6] CHIO D S, IN Y E, Prediction of shrinkage and warp age in consideration of residual stress in integrated simulation of injection molding[J]. Composite Structures,1999, 47(1/2/3/4):655-665.

[7] 王振军.基于ANSYS的铸件凝固过程瞬态温度场的有限元模拟[J].现代制造工程,2003(12):40-42.

(WANG Z J. Finite element simulation of transient temperature field in casting solidification process based on ANSYS[J]. Modern Manufacturing Engineering, 2003,12 :40-42.)

[8] 柳百成. 铸件充型凝固过程数值模拟国内外研究进展[J].铸造,1999(8):40-44.

(LIU B C. Research progress on numerical simulation of mold filling and solidification process at home and abroad[J]. Foundry, 1999(8):40-44.)

[9] 徐艳,康进武,黄天佑.铸造过程温度场-应力场双向耦合的数值模拟[J].清华大学学报(自然科学版), 2008,48(5):769-772.

(XU Y, KANG J W, HUANG T Y. Numerical simulation of coupling of temperature field and stress field in casting process[J]. Journal of Tsinghua University(Science and Technolog), 2008,48(5):769-772.)

[10] 王跃平.高温合金精铸叶片热应力力与变形的数值模拟[D].哈尔滨：哈尔滨工业大学,2007.

(WANG Y P. Numerical simulation of stress and deformation of the superalloy blade[D].Harbin： Harbin Institute of Technology,2007.)

[11] 董一巍. 净成形空心涡轮叶片精铸模具型腔优化设计方法研究[D]. 西安：西北工业大学, 2012: 557-81.

(DONG Y W. Study on the design method of cavity optimization for the hollow turbine blade[D]. Xi′an:Northwestern Polytechnical University, 2012: 557-81.)

[12] 李永毅. 精铸涡轮叶片凝固过程温度场实验及仿真验证[D]. 西安：西北工业大学, 2009： 3.

(LI Y Y. Simulation and experiment on solidification of the turbine blade[D]. Xi′an: Northwestern Polytechnical University,2009：3.)

[13] 杨屹, 蒋玉明, 刘力菱. 铸件凝固过程中热应力场及热裂的数值模拟研究分析[J]. 铸造技术, 2000(2): 36-38.

(YANG Y, JIANG Y M, LIU L L. Numerical simulation of shermal stress fields and crack in casting solicitation process[J]. Foundry Technology, 2000(2): 36-38.)

[14] KANG J W, L IU B C, XIONG S M. Numerical simulation of shaped casting based on rheological modal [J]. Journal Mater, 1999, 15(3): 267-269.