在数字化、智能化技术浪潮驱动下，传统铸造正加速新旧动能转换，如何通过技术创新推动行业转型升级，已成为提升我国铸造产业国际竞争力的关键[1-3]。
铸造过程的应力场数值模拟可划分为3 个典型阶段： 初期阶段以开发专用应力分析程序为特征，随后进入借助通用有限元软件进行应力分析的阶段， 最后发展至在铸造模拟软件中集成专门的应力计算模块的阶段[4]。 进入21 世纪以来，国外对铸造过程温度场/应力场求解的研究已经进入第Ⅲ阶段，ProCAST、MAGMA 等国外铸造仿真软件凭借成熟的多物理场耦合能力(流动-传热-应力)已形成完整技术生态，主导了全球市场。 2005 年，Lee 等[5]采用ProCAST 的热-力耦合分析功能， 研究了镍铝铜合金船用螺旋桨凝固过程的应力分布与位移情况，并通过测量铸件实际形变量验证了模拟结果的准确性。Yang 等[6]基于MAGMA 软件对Al-Cu 合金气缸盖火面进行热-力耦合分析， 发现火面较大的壁厚过渡差异导致凝固过程中温度和应力梯度较大，从而引发热裂缺陷。 Kim 等[7]利用ProCAST 软件对低压压铸高强度铝合金螺旋桨的多物理场进行分析，开发了耦合热流体应力分析模型来优化低压压铸工艺。 虽然国内已将数值模拟技术广泛应用于铸造工艺优化，但是鲜见能专用于铸造过程的温度场/应力场一体化自主开发仿真软件。 目前我国学界在铸造过程数值模拟研究时，普遍采用“分段式”联合求解策略：基于有限差分法自主开发充型流动与凝固传热模拟程序，将求得的温度场数据导入ANSYS、ABAQUS 等通用有限元平台进行应力计算，或者使用国外铸造模拟软件如ProCAST 或MAGMA 进行多物理场耦合求解[8-10]。这种方式虽能够暂时规避自主应力算法的技术瓶颈，但工艺数据在异构软件间反复迁移，会导致计算效率损失，同时面临参数泄露与系统兼容性风险。 更为关键的是，求解所依赖的本构模型和热-力耦合算法均被国外软件技术封锁，不仅制约我国铸造工艺自主优化能力，更对产业链安全构成潜在威胁。
针对工业软件自立自强的国家战略需求，提出一种基于统一网格的热-力单向耦合算法， 研究基于有限元法的铸造过程温度场/应力场一体化数值模拟技术， 自主开发专用于铸造的仿真工业软件，能够高效、 精确地计算温度变化和应力应变分布，并对铸件凝固过程中的热裂缺陷进行预测。
本文主要介绍基于有限元法的铸造过程温度场、应力场数学模型和铸件边界条件处理的关键技术，自主开发了有限元数值模拟系统。 最后对阶梯试块、应力框架模型和燃烧室壳体模型进行求解计算以验证自主开发软件的算法准确性和工程实用性。
1 基于有限元法的铸造温度场/应力场数学模型
1.1 铸造过程温度场有限元数值模拟的数学模型
1.1.1 三维瞬态温度场的变分过程
非稳态导热偏微分方程是温度场模拟计算的常见数学模型，三维瞬态温度场的场函数温度T(x,y,z,t)满足微分方程[11]：
式中，ρ 为材料密度；c 为材料比热容；t 为时间；λ 为热导率；Q 为热源密度。 为求解方程(1)，需要引入以下3 类边界条件和初始条件：
(1)已知Γr 边界上的温度值
(2)已知Γq 边界上的热流密度
(3)已知Γh 边界与周围介质的热交换
式中，nx,ny,nz 为边界外法线的方向余弦；T0 为第一类边界下指定的温度；q 为第二类边界的热流量；h为界面换热系数， 用于描述边界与周围环境之间的热交换效率；Ta 为环境温度。
单元内温度T 可近似由有限元节点温度Ti 插值得到，此时节点温度是时间的函数，即：
式中，ne 为有限元单元的节点个数；Ni(x,y,z)为空间域内的插值函数。
采用伽辽金加权余量法对温度场的热传导微分方程(1)进行计算，推导出三维瞬态热传导问题的有限元求解方程为：
式中，各矩阵由单元的相应矩阵元素集成：C 为比热容矩阵；K 为热导率矩阵；
为节点温度对时间的导数列阵；P 为温度载荷列阵； T¯为第一类边界中给定温度；nx, ny, nz 为边界外法线的方向余弦；q 为第二类边界的热流量；h 为界面换热系数；Ta 为环境温度。
1.1.2 求解瞬态热传导常微分方程组
式(6)是瞬态热传导问题的有限元求解方程，本文采用直接积分法求解微分方程组[12]：
利用上式， 从t=0 出发可以依次递推求解节点温度矩阵的各个瞬时值Tn，其中：
式中， 
为系统的有效系数矩阵;Tn 和Tn+1 分别是n和n+1 时刻的温度向量;
是有效载荷向量;Pn 和Pn+1 是n 和n+1 时刻的温度载荷向量;参数θ 是差分参数，本文采用伽辽金法(θ=2/3)计算温度场。
1.2 铸造过程应力场有限元数值模拟的数学模型
1.2.1 热弹塑性模型
采用热弹塑性模型， 在弹性阶段应力应变服从广义胡克定律[13]。 数学表达式为：
式中，σ 为应力张量；De 为弹性本构4 阶张量；ε为应变张量；μ 和λ 为拉梅系数。
类似于线弹性本构方程， 基于增量理论的热弹塑性模型可以用数学形式表达为：
式中，dσ 为应力增量；Dep 为弹塑性本构4 阶张量；dε 为应变增量。
引入弹性塑性转变系数m，式(15)中的弹塑性本构4 阶张量Dep 表示为：
式中，弹性阶段m 为1；塑性阶段m 为0；De 为式(14)中的弹性本构4 阶张量；Dp 为塑性本构四阶张量。
增量理论不局限于特定的加载条件， 可以同时满足弹性阶段和塑性阶段， 能够体现变形的历史过程。在应用时，只需沿着加载过程中的变形路径对时间t 进行积分，就可以得到最终的应力和应变状态。用数学表达式表示为：
1.2.2 应力场有限元法的离散及求解
为了方便温度载荷、 网格节点与单元的信息在温度场和应力场之间传递，温度场/应力场耦合的有限元求解系统前后采用同一套网格。 当弹性体在外力作用下发生变形时，根据能量守恒定律，则有：
式中，Πp 为该弹性体的总势能；Wi 为变形后应变能；σ为应力张量；ε 为应变张量；ε0 为初应变张量；Wf为外力做功，等于体积力和表面力做功之和，其中u为位移矢量矩阵；F为体积力矢量；T 为表面力矢量。
由于铸件在各时刻视为准静力平衡， 则根据最小势能原理，当Πp(u)最小时有：
式中，K 为刚度矩阵；P 为载荷项。 对于热弹塑性模型，式(21)需要修改为增量形式：
式中，{K}i 是第i 个增量步的刚度矩阵；{Δu}i 为第i个增量步时发生的位移矢量的增量；{P}i 为第i 个增量步的载荷向量。
2 铸造过程温度场/应力场数值模拟的关键技术
2.1 多材质温度场换热边界处理
图1 是换热边界处理的示意图， 图中蓝色部分是铸件，红色部分是铸型，绿色部分是其他铸造工艺元素，比如冷铁和型芯等，以换热边界最为复杂的铸件为例，铸件的面单元可以分为4 类：①与空气接触的面单元A1；②与铸型接触的面单元A2；③与其他材质接触的面单元A3；④铸件内部的面单元A4。
图1 铸造过程换热边界示意图
Fig.1 Schematic diagram of heat exchange boundaries in the casting process
面单元的区分共有3 步：
(1)根据单元出现频次区分外表面单元(A1/A2/A3)与内部面单元(A4)；
(2)基于面单元外法线方向筛选出与空气接触的面单元A1；
(3)采用立方体搜索算区分A2/A3 单元。
立方体搜索算法如图2 所示， 假设对一个目标外表面单元ABC 进行分析，对流换热边界是针对其单元形心点P 而言，可以根据形心P 的三维坐标位置扩展出一个长宽高为相同大小d 的正立方体，这个正立方体内部包含了哪些其他材质单元， 则判定这个目标外表面单元与哪些其他材质单元相邻。 为了能更快速地判定这个正立方体内部有哪些其他材质单元，以其他材质单元的所有节点坐标Ni(x,y,z)为标准， 根据各节点坐标在3 个方向上的最大最小值， 再构造一个长方体。 以图中其他材质单元JKL为例，该面单元由J、K、L3 个节点构成，根据NJ(x,y,z)、NK(x,y,z)、NL(x,y,z)的坐标，由3 节点相对位置可知x,y,z 坐标的最大值分别为NK(x)，NJ(y),NL(z)，最小值分别为NJ(x),NL(y),NK(z)。 通过形心扩展出的正立方体与其他材质单元节点扩展出来的长方体是否相交， 来判断目标外表面单元与哪些其他材质外表面单元相邻。 然后再进一步找出这些相邻其他材质外表面单元中， 将目标外表面单元最近的其他材质表面单元作为换热边界条件。 这种通过以空间换时间的方式能够快速便捷区分出目标外表面单元与何种材质单元相邻， 为边界条件设置提供必要的前置工作准备。
图2 立方体搜索算法示意图
Fig.2 Schematic diagram of the cube search algorithm
2.2 铸件的力学边界条件处理
采用边界单元法将铸型简化为边界单元， 其基本思想是将铸件与铸型复杂的相互作用替代为弹簧、棘轮和间隙等简单元件组成的表面单元，弹簧单元通过给铸件施加垂直于接触面的边界力来模拟铸型对铸件的支撑力，有效地简化复杂的物理交互，同时保持了数值模拟的精度和效率， 计算量较小且通用性强，能够模拟复杂的铸型结构，对于型芯和冷铁等其他工艺元素也能适用。
引入Chang[14]提出的铸型局部刚度系数概念，定量分析铸型对铸件产生的法向作用力， 并采用陈涛[15]提出的方法，根据铸件的几何模型尺寸估算一个铸型厚度，或者使用一个大致的平均铸型厚度，近似地计算铸型局部刚度系数。 铸型局部刚度系数用数学表达式可以表示为：
式中，E 为铸型的弹性模量；L 为铸件边界面法向上铸型的厚度，即是吃砂量。
通过铸型局部刚度系数， 可以得到铸件边界单元受到铸型的压力，数学表达式为：
式中，δ 为铸件边界单元沿着法向上的位移， 当δ＜0时表示铸件与铸型分离， 当δ＞0 时表示铸件与铸型接触，铸件受到与位移呈正相关的压力。
在1.2.2 节中已经完成铸造过程应力场有限元离散线性方程组的建立， 为了获得符合实际的唯一解，接下来需要对铸件边界单元设置合理约束。结合上一小节介绍的边界单元法，将式(26)代入式(20)中，得到铸件表面单元被施加的离散后的表面力为：
式中，n 为铸件外表面单元的外法线方向，号表示受到表面力的作用与外法线方向相反，{du}i=N{Δu}i为第i 个增量步铸件表面单元的位移增量。
对于βm，同式(26)中一样，需要根据铸件与接触材质的相对位置进行设置， 假设当前增量子步与前一增量子步的相对位置没有发生变化， 那么可以有数学表达式为：
式中，{u}i-1 表示直到i-1 个增量子步时刻为止，之前所有增量子步的累计位移量。
这时式(22)里右端载荷项中含有表面力的一项与位移增量{Δu}i 有关，根据前文中最小势能原理推导过程，需要移到左端刚度矩阵中，最终的刚度矩阵和载荷项分别为：
2.3 基于热弹塑性模型的热裂判据
热裂作为凝固过程的关键缺陷之一， 源于准固态区枝晶骨架形成阶段的补缩失效， 当枝晶网络承受的拉应力超越其强度极限时， 液膜破裂引发裂纹萌生[16-17]。 该区域材料因高温低模量特性，微小应力变化即可诱发显著应变， 故需综合等效应变累积与三向应力状态建立判据。 采用陈涛[15]提出的热裂倾向判据C 热裂，当σi＜0 的时候，表示最大的主应力为负，那么第二、三主应力一定也是负值，此时材料的应力状态为压应力，认为热裂倾向较低。热裂判据C 热裂的数学表达式为：
式中，σ1 为第一主应力； 
为高温阶段准固态区发生的等效应变，最终C 热裂的值越大，表示该处热裂倾向较大，反之C 热裂的值越小，表示该处热裂倾向较小。
2.4 自主开发软件的计算流程
在铸造的凝固过程中， 温度场与应力场之间存在复杂的热-力双向耦合关系，二者相互作用、相互影响。 徐艳等[18]和傅显钧[19]针对铸造过程，对比分析了热-力双向耦合过程和热-力单向耦合过程，发现二者的温度计算结果趋势相同，数值范围差别不大。由于铸造凝固过程通常表现为静态变形， 变形速率较低， 因此应力变形做功及界面热阻变化对热效应的影响相较于温度变化和潜热释放来说较小。
基于有限元数学模型和边界条件处理模型，本文采用热-力单向耦合方法计算铸造过程的温度场和应力场：首先进行铸造过程温度场数值模拟，求得温度场数据之后以此作为热载荷， 再进一步计算铸造过程中的应力场。 具体采用C++ 语言，在Visual Studio 平台上编写有限元软件求解器， 自主开发了有限元铸造过程温度场/应力场数值模拟软件，图3是整个软件计算流程示意图。
图3 基于有限元法的铸造过程温度场/应力场数值模拟模块流程图
Fig.3 Flowchart of the numerical simulation module for temperature/stress fields in the casting process based on the finite element method
3 程序计算结果与分析
3.1 铸造数值模拟系统温度场结果算法验证
3.1.1 阶梯试块模型及参数设置
为了确保有限元温度场程序的精确度， 对阶梯试块的空冷过程计算， 再与ANSYS 模拟结果进行详细对比分析。铸件材质为TA15 钛合金，其成分见表1。 如图4 是对阶梯试块模型进行剖分的六面体一次单元网格，其中节点总数为29481，单元总数为25650，单元尺寸为10 mm。 设置试块初始温度1 350 ℃，环境温度20 ℃，阶梯试块和环境空气的对流换热系数为35 W/(m2·℃)，试块的物性参数设置见表2。
表1 TA15合金化学成分
Tab.1 Chemical composition of the TA15 alloy(mass fraction/%)
AlVZr MoFeSiCONH 6.61 2.33 1.87 1.78 0.065 0.017 0.019 0.11 0.016 0.003 1
表2 TA15钛合金物性参数
Tab.2 Physical properties of the TA15 alloy
Thermal conductivity/(W·m-1·℃-1)Specific heat/(J·kg-1·℃-1)Density/(kg·m-3)46.057004 300
图4 六面体网格剖分图
Fig.4 Hexahedral mesh discretization diagram
3.1.2 与ANSYS 模拟结果对比
如图5a～c 是自主开发软件的温度场计算结果，图5d～f 是ANSYS 计算的温度场结果。 由图可见，两种软件计算结果高度吻合， 在阶梯试块的前端壁厚最薄的位置降温最快，后端降温最慢，整体温度场呈对称分布， 自主开发软件计算结果符合实际的温度场分布规律。
图5 温度场计算结果对比：(a～c) 自主开发软件；(d～f)ANSYS
Fig.5 Comparison of the temperature field results:(a～c)self-developed software;(d～f)ANSYS
选取阶梯试块表面两个网格节点A、B，两个节点在阶梯模型上的位置见图6。 对比自主开发软件得到的瞬态温度场计算结果与使用ANSYS 模拟的结果：A、B 两点在t=0 s 到t=10 000 s 之间的温度变化曲线分别见图7a 和b。 从图中可以看出，二者降温趋势一致，在开始时间降温迅速，之后逐渐平缓，冷却速率随时间逐渐降低。
图6 A、B 点在阶梯试块中的位置
Fig.6 Locations of points A and B in the step-shaped test block
图7 阶梯试块上两点的温度变化曲线：(a)A 点；(b)B 点
Fig.7 Temperature variation curves at two points on the step-shaped test block:(a)point A;(b)point B
表3 记录了自主开发软件和ANSYS 计算得出的各时刻最大温度和最低温度，通过对比结果来看，温度数据最大相对误差为0.789%。 由于ANSYS 使用的求解器与本文不同， 并且在有限元方程的计算求解过程中存在近似处理，作为数值模拟的结果，可以认为二者结果一致， 误差在实际工程应用中为可接受范围内。
表3 最大温度、最小温度各时刻的温度对比
Tab.3 Temperature comparison of the maximum/minimum values at different time instants
Time/sMinimum temperature/℃Self-developed softwareANSYSError/%Self-developed softwareANSYSError/%101 350.001 350.000.0001 260.681 262.620.154 201 349.971 349.810.0121 219.071 221.470.197 501 348.291 347.200.0811 128.061 130.430.209 1001 335.581 333.890.1271 013.401 015.930.249 2001 286.271 285.150.087841.26844.320.363 3001 225.441 224.830.050711.64714.940.463 5001 098.051 098.330.025526.93530.160.608 1 000817.96819.210.158281.88284.120.789 2 000447.79449.290.335112.08112.990.802 3 000246.66247.940.51959.0959.500.686 5 00082.4583.080.76528.9129.020.381 7 50032.2932.480.59821.6421.670.126 10 00022.4122.460.22920.3220.320.033 Maximum temperature/℃
3.2 栅形应力框架模型的温度场和应力场数值模拟
3.2.1 栅形应力框架的几何模型与前处理
选用经典的栅形应力框架铸件模型作为模拟对象， 计算分析铸件在凝固过程中的温度场和应力场分布。 图8 是栅形应力框架铸件的几何模型。 图9a和b 分别是铸件和铸型的四面体网格， 其中节点数量为37650，单元数量为192893，设置单元最大尺寸为5 mm， 并对中间粗杆和两侧细杆部分进行局部自适应加密，设置加密最大尺寸为3 mm。
图8 栅形应力框架铸件的几何模型
Fig.8 Geometric model of the grid-shaped stress framework
图9 栅形应力框架的四面体网格：(a)铸件；(b)铸型
Fig.9 Tetrahedral mesh of the grid-shaped stress framework:(a)casting;(b)mold
3.2.2 栅形应力框架的温度场计算
栅形应力框架的材质为常用结构钢ZG25， 合金成分(质量分数，%)为：0.25 C、0.45 Si、0.9 Mn、0.1 P、0.1 S，铸型材质为自硬呋喃树脂砂。 铸件的热物性参数和力学性能通过JMatPro 软件计算得到，如图10 所示。 铸件初始温度为1 550 ℃，铸型初始温度为30 ℃。
图10 ZG25 物性参数随温度变化的曲线：(a)密度；(b)比热容；(c)泊松比；(d)热导率；(e)热膨胀系数；(f)杨氏模量
Fig.10 Temperature-dependent thermophysical properties of ZG25 steel:(a)density;(b)specific heat;(c)Poisson's ratio;(d)thermal conductivity;(e)coefficient of thermal expansion;(f)Young's modulus
由于栅形应力框架模型的主要研究对象是粗杆和细杆的热应力变化，所以选择粗杆中心处和一侧的细杆中心处作为分析点，分析点位置如图11所示。
图11 栅形应力框架铸件上分析点的选取
Fig.11 Selected analysis points on the grid-shaped stress framework
图12 是粗杆中心和细杆中心的温度变化曲线图。可以看见二者温度变化趋势差异明显：粗杆中心温度曲线呈现平缓下降趋势，整体降温速率较低；细杆中心温度则表现为初期快速下降后趋于稳定。
图12 栅形应力框架铸件上粗杆中心和细杆中心的温度变化曲线
Fig.12 Temperature evolution at centers of thick/thin rods in a grid-shaped stress framework
图13a～f 是栅形应力框架模型在六个不同时刻的温度场分布，从温度场分布来看，模拟结果符合实际傅里叶热传导定律和牛顿冷却定律。
图13 栅形应力框架在不同时刻的温度场分布图：(a)100 s;(b)300 s;(c)600 s;(d)1 000 s;(e)1 600 s;(f)2 500 s
Fig.13 Temperature field distributions of the grid-shaped stress framework at different time instants:(a)100 s;(b)300 s;(c)600 s;(d)1 000 s;(e)1 600 s;(f)2 500 s
3.2.3 栅形应力框架的位移场和应力场计算
如图14 和15 是基于不同的铸型刚度系数得到栅形应力框架在冷却过程中的变形情况， 并对栅形应力框架沿y 方向(即沿杆的轴向)的正应力stress1 进行分析， 可以发现栅形应力框架的变形与应力场分布呈现显著的热-力耦合特征， 其中铸型刚度系数如式(25)所示，用于表示铸型硬度。
图14 较软铸型条件下栅形应力框架中心面的变形情况和y方向上的应力分布
Fig.14 Deformation and y-direction stress distribution at the midplane under the low-stiffness mold
如图14 所示， 此时铸型刚度系数取βm=1×106，对应铸型较软的情况。 凝固初期： 细杆因快速冷却引发显著收缩，而粗杆中心仍处于高温状态，所以细杆相对粗杆收缩明显，变形量远大于粗杆。由于铸型柔性约束，细杆通过向心弯曲释放变形能，内侧因曲率差异形成轴向拉应力集中； 粗杆外侧受框架端部约束产生压应力，内部越接近中心，拉应力越大。 凝固中期：粗杆收缩加剧主导力学行为，细杆因前期塑性变形长度增加，导致其应力状态由拉转压；粗杆外侧内部拉应力持续升高，形成y 方向上“内拉外压”分布特征。凝固末期：框架整体尺寸收缩使细杆完全受压，此时细杆完全转变为压应力状态，弯曲程度逐渐恢复至初始状态，接近最初的笔直状态；粗杆因持续凝固收缩呈现轴向拉应力主导状态， 表面因为框架端部发生了塑性变形而受到压应力作用。
如图15，此时铸型刚度系数取βm=1×109，对应较硬的铸型情况。凝固初期，硬铸型限制细杆自由收缩，导致其轴向拉应力与径向压缩载荷共存。同时应力框架四角降温收缩，产生向心弯曲力矩，该力矩通过节点传递至细杆两端，形成径向压缩载荷，迫使细杆端部向内收缩，进而挤压细杆向外屈曲，弯曲内侧受压应力；粗杆应力分布与软铸型初期相似，呈现外侧压应力与内部拉应力梯度。凝固中期，粗杆主导收缩使框架整体尺寸减小，细杆拉应力反转至压应力。粗杆内部拉应力持续积累， 表层因端部塑性变形维持压应力。 凝固末期，整个应力框架尺寸继续变小，细杆被强制压缩至更高屈曲模态， 弯曲内侧压应力达到峰值。粗杆轴向拉应力进一步增大，表面压应力仍有部分残留。
图15 较硬铸型条件下栅形应力框架中心面的变形情况和y方向上的应力分布
Fig.15 Deformation and y-direction stress distribution at the midplane under the high-stiffness mold
图16a 和b 分别展示了铸型刚度系数1×106 为1×109 和时栅形应力框架的热裂判据值分布， 图17则呈现了两种铸型条件下的等效应力和等效应变云图。在刚度较低的铸型条件下，粗杆与框架连接处因整体收缩引发局部拉应力集中，热裂倾向较大；当铸型刚度较高时， 细杆端部因高拉应力集中成为主要热裂风险区。二者差异源于铸型的约束特性：硬铸型限制变形导致细杆端部持续高应力； 软铸型通过端部塑性变形转移应力集中至粗杆节点。
图16 栅形应力框架在不同铸型条件下的热裂判据值C 热裂：(a)铸型刚度系数βm=1×106；(b)铸型刚度系数βm=1×109
Fig.16 Hot tearing criterion of the grid-shaped stress framework under various mold conditions:(a)βm=1×106;(b)βm=1×109
图17 栅形应力框架在不同铸型条件下的等效应力和等效应变：(a)铸型刚度系数βm=1×106；(b)铸型刚度系数βm=1×109
Fig.17 Equivalent stress and strain under different mold stiffnesses:(a)βm=1×106;(b)βm=1×109
如图18[20]，实际栅形应力框架铸件在细杆端部区域观察到典型热裂纹， 与模拟预测的细杆端部热裂高风险区空间分布高度吻合。 该现象验证了铸造过程热-力耦合数值模拟软件的有效性， 能够复现铸件凝固过程中瞬态应力场演化行为和应力集中现象， 成功实现了刚性与柔性铸型边界条件的差异化建模， 精确表征了细杆端部因约束作用形成的拉应力集中现象。 此外， 模型能够可靠预测复杂铸件的热裂倾向，为铸造工艺提供指导性优化依据。
图18 实际铸造中栅形应力框架出现的热裂缺陷[20]
Fig.18 Actual hot crack defects in the grid-shaped stress framework[20]
3.3 铸钢件燃烧室壳体的温度场和应力场数值模拟
3.3.1 燃烧室壳体的几何模型和前处理
为了验证自主开发软件的工程实用性， 适用于实际铸造场景， 对大型复杂燃机铸钢件的关键部件燃烧室壳体部分进行模拟计算。图19 是燃烧室壳体的实物图和模型，出于保密原则，此处只显示铸件和冷铁部分。
图19 燃烧室壳体：(a)实际铸件；(b)模型
Fig.19 Combustion chamber casting:(a)actual casting;(b)simulation model
3.3.2 燃烧室壳体的温度场、 应力场数值模拟和热裂缺陷预测
图20 是燃烧室壳体在4 个时刻的温度分布，可以看见设置冷铁的区域产生显著激冷效应， 导致冷铁作用区局部冷却速率较大。 薄壁部位表面积与体积的比值大，散热效率高，热量能够快速向周围环境传递，所以薄壁部位降温速率较快。冒口部位由于保温作用和高热容量维持高温，降温较慢，是整个铸件温度最高的部位。
图20 燃烧室壳体不同时刻铸件温度分布：(a)125 s;(b)749 s;(c)1 417 s;(d)3 000 s
Fig.20 Temperature distributions of the combustion chamber at different characteristic time instants:(a)125 s;(b)749 s;(c)1 417 s;(d)3 000 s
图21 是燃烧室壳体热裂倾向较大的部位，主要分为4 类：①铸件支撑结构端部；②冒口和铸件交界处；③铸件通孔处；④铸件薄壁区域。 由于燃烧室壳体在浇注之后会去掉铸件支撑结构和冒口， 并进一步处理通孔结构， 所以本文主要讨论铸件薄壁区域的热裂倾向。
图21 燃烧室壳体出现热裂倾向较大的部位
Fig.21 High hot-tearing tendency zones in the combustion chamber
图22 是燃烧室壳体在t=3 000 s 时第一主应力的大小， 可以看见薄壁区域热裂倾向较大的部位主应力是正值，表明该处承受三向拉应力作用。其成因可归结于薄壁部位与周围结构的冷却速率差异：薄壁部分降温速率较快，率先进行凝固收缩，变形量较大，四周降温速率慢，会对薄壁部分有约束阻碍作用， 导致拉应力累积。 当薄壁区域承受的拉应力超过材料的高温强度极限时， 热裂纹就可能萌生扩展。 图23 是燃烧室壳体的高温阶段等效应变分布，可以看见在冷铁的激冷作用下，薄壁部分出现等效应变集中现象，且应变集中区域与拉应力集中区高度重合。 该现象验证了热裂判据的合理性，可有效表征热裂倾向的力学本质。
图22 燃烧室壳体在t=3 000 s 时第一主应力大小
Fig.22 First principal stress distribution of the combustion chamber at t=3 000 s
图23 燃烧室壳体高温阶段等效应变分布
Fig.23 Equivalent strain distribution during the high-temperature stage
如图24 所示，燃烧室壳体在实际生产中于薄壁弧面区域出现了显著的热裂纹缺陷，这与数值模拟结果预测的高拉应力集中区、等效应变峰值区还有热裂倾向较大的区域吻合。 表明自研仿真软件能够精确捕捉薄壁结构因凝固收缩受阻引发的拉应力累积效应，进而可靠预测热裂缺陷的萌生位置。 根据以上分析，可以判断出薄壁弧面位置出现热裂缺陷的原因主要来自局部过冷导致的拉应力集中现象，所以在实际工艺中可以采用冷铁梯度布局和降低浇铸温度等方式进行优化。
图24 实际生产中燃烧室壳体出现热裂缺陷图片
Fig.24 Actual hot crack defects in the combustion chamber casting
4 结论
(1)基于铸造有限元法的温度场、应力场数学模型，提出一种基于统一网格的热-力单向耦合方法，自主开发了能够求解铸造过程温度场/应力场的数值模拟系统，并对铸件的热裂缺陷倾向进行预测。
(2)对阶梯试块的温度场进行计算，并与ANSYS计算结果进行对比，二者结果基本相同，最大相对误差0.789%，验证了自主开发软件能够对温度场进行准确求解； 对经典栅形应力框架模型的温度场和应力场进行求解计算， 分析应力框架的变形情况和应力应变分布，并预测热裂倾向，与实际栅形应力框架出现热裂部位相同， 验证了本文的铸件力学边界处理能够准确模拟铸件的受力情况； 最后对燃烧室壳体的温度场、应力场进行模拟，预测了热裂倾向，验证了自主开发软件的工程实用性。
[1] 刘小龙. 我国铸造装备的创新、智能、绿色发展之路[J]. 中国铸造装备与技术，2020,55(2):5-9.LIU X L. Innovation, intelligence and green development of foundry equipment in China[J]. China Foundry Machinery &Technology,2020,55(2):5-9.
[2] 曾广峰. 我国智能制造行业发展现状及趋势[J]. 质量与认证，2020(11):46-47.ZENG G F. Current status and trends of china's intelligent manufacturing industry[J]. China Quality Certification, 2020(11):46-47.
[3] 许庆彦.铸造技术路线图:数字化、网络化、智能化铸造[J]. 铸造，2017,66(12):1243-1250.XU Q Y.Foundry technology roadmap:digitalized,networked and intelligent casting[J].Foundry,2017,66(12):1243-1250.
[4] 胡翔宇，魏彦鹏，孙逊，于波. 铸造应力场数值模拟的研究现状与发展趋势[A].2017 中国铸造活动周论文集[C]. 苏州：中国机械工程学会，2017:480-485.HU X Y, WEI Y P, SUN X, YU B. The research status and development tendency of numerical simulation of casting stress field[A]. Foundry Week 2017[C]. Suzhou: Chinese Mechanical Engineering Society,2017:480-485.
[5] LEE S M,LEE W J.Finite-element analysis on thermomechanical behavior of a marine propeller casting in the sand-casting process[J]. Journal of Materials Engineering and Performance, 2005, 14(3):388-394.
[6] YANG C L, XU H, WANG Y, LI F, MAO H K, LIU S Y. Hot tearing analysis and process optimisation of the fire face of Al-Cu alloy cylinder head based on MAGMA numerical simulation[J].Materials Technology,2023,38(1):21652453.
[7] KIM M S, WON S. Development of low-pressure die-cast Al-Zn-Mg-Cu alloy propellers part II:Simulations for process otimization[J].Materials,2024,17(16):4027.
[8] 史东丽，钱坤才. 大型高铬铸铁叶轮凝固过程应力场数值模拟及铸造工艺优化[J].铸造，2020,69(9):972-976.SHI D L, QIAN K C. Casting process optimization of large-scale high chromium cast iron impeller based on stress field simulation of solidification process[J].Foundry,2020,69(9):972-976.
[9] 刘宇鹏，李青春，徐磊，李响，胡小强. 超长耐热钢管离心铸造工艺模拟与优化[J].特种铸造及有色合金，2023,43(7):955-958.LIU Y P, LI Q C, XU L, LI X, HU X Q. Simulation and optimization on centrifugal casting technology of ultra-long heat-resistant steel pipe[J]. Special Casting & Nonferrous Alloys,2023,43(7):955-958.
[10] ZHU Y K,XU G Q,XU C S,LI X L,ZHOU K K,CHEN Y,LI Y T, XU X S, LIU H. Hot cracks in camshaft casting: Initiation and propagation[J].International Journal of Cast Metals Research,2023,37(1):48-70.
[11] 曹流，廖敦明，曹腊梅，谷怀鹏，陈涛，庞盛永. 基于有限元法的熔模铸造过程温度场模拟软件自主开发[J]. 铸造，2014,63(12):1235-1240.CAO L, LIAO D M, CAO L M, GU H P, CHEN T, PANG S Y.Temperature-field simulation software self-development of investment casting based on finite element method[J].Foundry,2014,63(12):1235-1240.
[12] 王勖成.有限单元法[M].北京：清华大学出版社，2003.WANG X C. Finite element method [M]. Beijing: Tsinghua University Press,2003.
[13] CELENTANO D, OLLER S, OÑATE E. A coupled thermomechanical model for the solidification of cast metals[J].International Journal of Solids and Structures,1996,33(5):647-673.
[14] CHANG A, DANTZIG J. Improved sand surface element for residual stress determination[J].Applied Mathematical Modelling,2003,28(6):533-546.
[15] 陈涛. 基于有限元法的铸造热应力数值模拟及其智能化技术的研究[D].武汉： 华中科技大学，2013.CHEN T. Numerical simulation of casting thermal stress based on finite element method and intelligent techniques[D]. Wuhan:Huazhong University of Science and Technology,2013.
[16] HATAMI N, BABAEI R, DADASHZADEH M, DAVAMI P.Modeling of hot tearing formation during solidification[J].Journal of Materials Processing Technology,2007,205(1-3):506-513.
[17] ESKIN D G, KATGERMAN L. A quest for a new hot tearing criterion[J].Metallurgical and Materials Transactions A, 2007, 38(7):1511-1519.
[18] 徐艳，康进武，黄天佑.铸造过程温度场/应力场双向耦合的数值模拟[J].清华大学学报(自然科学版)，2008(5):769-772.XU Y,KANG J W,HUANG T Y.Numerical simulation of two-way coupling between temperature and stress fields in castingprocess[J].Journal of Tsinghua University(Science and Technology),2008(5):769-772.
[19] 傅显钧. 基于ANSYS 的铸造过程热应力场双向耦合的数值模拟[D].武汉： 华中科技大学，2011.FU X J. Numerical simulation of the casting process with coupled thermal and mechanical effects based on ANSYS [D].Wuhan:Huazhong University of Science and Technology,2011.
[20] 廖敦明，刘瑞祥，陈立亮，林汉同. 基于有限差分法的铸件凝固过程热应力场数值模拟的研究[J].铸造，2003(6):420-425.LIAO D M, LIU R X, CHEN L L. Study on numerical simulation of thermal stresses during casting's solidification process based on FDM[J].Foundry,2003(6):420-425.