基于ProCAST的熔铸装药顺序凝固工艺边界建模方法研究

黄求安,刘恒著,史慧芳,曹红松,王兆国,牛凯博,杨 阳

(1.中国兵器装备集团自动化研究所有限公司， 四川 绵阳 621000；2.中北大学 机电工程学院， 太原 030051； 3.山东北方民爆器材有限公司, 山东 淄博 255208)

1 引言

熔铸装药工艺极易形成缩孔缩松、气孔、裂纹等装药缺陷[1-2]。通过采用顺序凝固工艺，控制固相界面自下而上顺序生长，可极大地降低这些缺陷在弹体内部形成的概率[3-4]，因而受到广泛关注[5-7]。由于该工艺研发物理试验周期长、成本高、安全风险大，且相关测试设备只能获取装药成型过程时空离散点信息[8]，通过ProCAST仿真软件对药浆凝固过程进行数值模拟，全面获取工艺成型过程各物理场的动态变化过程，为工艺优化提供依据，已成为工艺研发的一项重要手段[9-11]。

利用ProCAST进行顺序凝固工艺仿真的关键步骤之一在于对弹体模具与冷却液面相对运动的工艺边界建模。通过对马松[12]、孙鑫欢[13]、高丰[14]等学者的文献分析发现，目前对该工艺边界建模主要通过对弹体外壁逐层配置耦合运动及热力学信息的跳转曲线，来实现模具外壁随液面升高，其换热对象由空气向冷却液的逐层转换。该方法主要存在建模过程极其繁琐庞杂，建模参数对耦合度高等问题。为此，本文构建了一种基于移动边界的顺序凝固工艺边界建模方法，详细阐述了该方法的主要实现流程，并对其进行了仿真验证。在不显著增加计算量的前提下该方法能够大幅度降低了工艺边界建模工作量，可以为顺序凝固工艺的研发与应用提供技术支撑。

2 顺序凝固工艺传热过程建模

2.1 顺序凝固工艺机理分析

顺序凝固工艺过程如图1所示。

图1 顺序凝固工艺过程示意图Fig.1 Schematic diagram of the sequential solidification process

由图1可知，该工艺将刚注完药的模具缓慢放入冷却液中(或者如图1所示的控制冷却液面缓慢上升)，利用冷却液相对较低的温度以及较高的热量传递效率，使入水部分快速冷却凝固，同时配合冒口保温、环境保温等措施，使模具未入水部分以及冒口区域的药浆不至于完全凝固，而药浆相界面以“U”或“V”字型随着冷却液面上升而逐渐向上生长回缩，从而实现自下而上的顺序凝固[3,12]。整个过程中，由于补缩通道始终保持通畅，模具内侧形成致密紧实、无缺陷的装药，而缩孔、缩松等缺陷通常只发生在冒口区域。

2.2 顺序凝固工艺基本传热过程数学建模

顺序凝固工艺主要涉及2个基本物理过程：其一是热量从药浆通过模具向外界的传递；其二是药浆凝固相变阶段结晶潜热的释放。对此，目前主要通过瞬态热传导方程予以描述[3]，即：

式中：ρ为密度；ch为比热；L为结晶潜热；fs为固相率；T为温度；λ为热传导率；t为时间。

在本文中，对结晶潜热释放的处理，主要通过引入固相率fs对温度T的偏导数项∂fs/∂T予以实现。该偏导数 ∂fs/∂T可通过对类似图2固相率-温度曲线进行数值求导来获取。从图2可以看出，只有当药浆处于凝固相变阶段(77～81 ℃)，∂fs/∂T不为零，从而引入了结晶潜热释放带来的热量。除该方法外，还可利用热焓方法处理结晶潜热释放，该方法需要炸药材料的热焓-温度曲线，具体可参考文献[15]。

图2 某TNT基熔铸炸药固相率-温度曲线Fig.2 Solid fraction curve of a TNT based melt-cast explosive

为实现对式(1)的求解，还需补充初值条件以及边界条件。其中，初值条件主要指是药浆、模具初始时刻的温度场；而边界条件主要是指药浆与模具、模具与环境等不同物质界面之间的热交换，这里主要通过流体与固体壁面的对流换热条件予以描述[3]，有：

(2)

式(2)中：Γ表示边界条件；h为换热系数；Tw与Ts分别表示模具外表面与周围介质的温度。

对于缩孔缩松的预测，直接采用ProCAST提供的临界固相率和补缩距离相结合的方法，具体可参考文献[16]，这里不在赘述。

3 顺序凝固工艺边界建模方法研究

3.1 传统顺序凝固工艺边界建模方法

利用ProCAST进行顺序凝固工艺仿真的关键步骤之一在于对弹体模具与冷却液面相对运动的工艺边界建模。而直接模拟弹体缓慢浸入冷却液(或液面上升)的过程，涉及实体的交互以及网格重划分等技术，这在目前ProCAST中是难以实现的。实际上，如果不关注冷却液内部的动量及能量变化，则可将冷却液简化为一项边界换热条件。在边界建模中，只需实现随着液面的上升，模具外壁对流换热对象从空气到冷却液的转换即可。由式(2)可知，这种转换实际就是换热系数h以及周围介质温度的转换。

通过对这种思想的直接应用，便形成了目前应用较多的一种顺序凝固工艺边界建模方法[12-14]。该方法并不建立任何实际的冷却液水浴模型，而是根据弹体对液面的相对运动速度，通过逐层配置换热系数及温度跳转曲线，抽象地实现随液面上升，模具外壁对流换热对象从空气跳转为冷却液的过程。图3为传统顺序凝固工艺边界建模方法示意图。该方法主要流程如下。

首先，在完成网格划分、材料参数设置等常规基本仿真配置后，将弹体外表面划分为如图3所示的数层环形区域。这一步实际是在有限元网格划分基础上，将液面上升过程进一步离散为数个具有一定厚度的环，而每层环的厚度就相当于数值计算意义上液面上升的“步长”。

图3 传统顺序凝固工艺边界建模方法示意图Fig.3 Schematic diagram of the conventional process boundary modeling method

其次，根据冷却水液面上升的相对速度以及环厚度，计算出液面扫过每层环的时间，并据此建立如图4所示的温度及换热系数跳转曲线。这种曲线的构建方法并不唯一，可以将液面到达某层环中间点的时刻作为换热对象跳转时刻。注意到这种跳转曲线的间断性(如图4所示，采用油浴或者水浴，其换热系数可能是103W/(m2·K)或更高,而空气通常为10 W/(m2·K)或者更小)，为防止步长无限缩小导致仿真停止，通常还需要对垂直段添加微小斜率来保证其计算正常进行。

图4 跳转曲线Fig.4 The jumping curve

最后，针对每条换热系数与温度跳转曲线，逐层选择其对应的环形区域所覆盖的弹体外表面网格，从而完成工艺边界建模。

对于这种传统顺序凝固工艺边界建模方法，主要存在以下2个缺陷：

第一，建模工作繁琐且庞杂。该方法中弹体模具外表面换热对象的转换依赖于手动逐层建立的换热及温度跳转曲线。以图3为例，共划分了57层环，亦即，在一次仿真中就需要建立57×2条类似图4的跳转曲线，并且需要通过鼠标框选特定网格等操作，将这些曲线与环形区域逐层关联。其次，由于这些跳转曲线耦合了液面上升的运动信息以及边界换热信息，因此，任一工艺参数的改变都需要对所有层的曲线重新配置。在工艺正交仿真等研究中，其建模工作量是相当具有挑战性的。

第二，工艺边界建模方法与仿真精度强耦合。如果“步长”过大(环形区域厚度过大)，则不可避免地会降低计算精度。一方面，为了提高精度，必须保证一定数量的环形区域层数(难以减轻建模工作量)；另一方面，由于环形区域的划分又是基于有限元网格的，因此，网格的尺寸又对环形区域数量的上限进行了限制。这样，在某些需要利用粗糙网格进行快速估算的场景中，基于该方法的仿真精度可能会受到限制。

3.2 基于移动边界的顺序凝固工艺边界建模方法

为解决3.1节中问题，利用ProCAST的Translate_x(t)功能构建了一种基于移动边界的工艺边界建模方法。该方法主要流程如下：

首先，需要额外建立一个如图5所示的建模方法构建的筒状水浴模型。该模型主要包含一个水托盘以及一个中空的筒状水浴。该模型主要用于模拟图1中所示的冷却液。其内径尺寸等由弹体模具所确定。

图5 基于移动边界的顺序凝固工艺边界 建模方法示意图Fig.5 Schematic diagram of the the process boundary modeling method based on moving boundary

然后，利用位于Geometrical下的Translate_x(t)功能建立筒状水浴的移动边界，从而实现弹体模具与水浴模型的相对运动。建立该移动边界的主要工作是在Translate_x(t)中定义水托盘及水筒的位移曲线(在本文中，为方便后续处理中输出药浆内部某个固定点的温度曲线，选择让模具固定，水托盘与水筒运动)。表1分别为定义水托盘和水筒运动的位移数据示例[12]，图6为筒状水浴模型按照该位移数据运动过程示意图。

表1 筒状水浴模型位移数据示例

结合表1与图6可以看出，该数据以模具底部为零点，水托盘和水筒初始位置定义在-10 mm处。测试发现，即使从弹体模具与液面接触时刻开始仿真，也不能在建模时让弹体与水面在初始时刻就保持接触，否则可能会在模具、水筒与水托盘接触区域形成共节点，并且随水筒的运动而被向上拖拽，从而产生网格畸变。为此，表1中设置0.1 s后筒状水浴就迅速运动到模具底部，近似模拟弹体与水面在初始时刻保持接触。当筒状水浴与模具底部接触后，水托盘则保持静止，而水筒则继续匀速向上移动。经过4 500 s后,水筒上沿升高至400 mm的位置(弹体模具高度),从而完全包裹住弹体模具。由表1数据可知,这里液面上升相对速率为320 mm/h。

图6 筒状水浴运动过程示意图Fig.6 Schematic diagram of the movement of tubular water bath model

接下来是对筒状水浴模型设置温度及换热边界条件。在温度方面，只需根据工艺参数对水托盘与水筒的外表面(如图7所示)配置某一恒定温度或者温度-时间曲线即可。在换热系数设置方面(见图8)，弹体模具与水筒、水托盘初始时刻没有接触，因此需要手动添加对流换热条件，并设置成NCOINC类型。当筒状水浴模型与弹体外壁发生接触后，弹体模具外表面计算单元的换热对象可自动转换为冷却液。对于水筒和水托盘，由于两者外表面均为某一恒温，因此，可以保持默认的COINC类型，并且将换热系数设置为0。

图7 设定温度条件的水筒及水托盘表面Fig.7 Setting the surface temperature of the tube and tray

图8 换热系数设定界面Fig.8 Setting the heat transfer coefficient

以上即为基于移动边界的顺序凝固工艺边界建模方法的主要流程。从以上建模流程可以看出：

首先，该方法中弹体外壁换热对象的转换通过移动边界自动实现。同时，液面相对运动与边界换热条件的设置相互独立，建模工作量小，在多方案仿真研究中，只需根据仿真试验设计建立对应的位移曲线以及边界换热条件即可。

其次，该方法中液面上升的“步长”是由ProCAST内部自适应时间步长以及水浴模型位移曲线所决定，没有额外的人为设定“步长”问题。因此，基于该方法的仿真计算精度主要由有限元网格划分精细程度以及时间步长所决定，而该方法本身只提供一种可移动的边界换热模式。

此外，由于该方法中只利用水浴模型外表面单元，而内部单元不参与热交换，因此筒状水浴模型的壁厚对仿真精度的影响小，在建立筒状水浴几何模型时可以将壁厚尽可能减薄，从而使得筒状水浴模型在计算量上不会产生过多的额外负担。

4 仿真验证与结果讨论

该部分主要利用文献[12]的相关数据对本文中所构建的基于移动边界的顺序凝固工艺边界建模方法进行仿真验证。在文献[12]中，其工艺边界建模方法为2.1节所述的传统方法。

在该仿真验证中，所选择炸药为HBX-3TNT基含铝炸药，模具为Al，冒口为PTFE。空气与模具外壁的换热系数为0.01 kW/(m2·K)，冷却液选择水作为介质，冷却液与模具外壁的换热系数为3.0 kW/(m2·K)。其他材料参数等可参考文献[12]，这里不在列出。仿真试验的主要工艺参数如表2所示[12]。其中，保温结束后，再经过4 500 s匀速下降至25 ℃。

表2 仿真验证工艺参数

图9为仿真验证所使用的几何模型。模型总高度为500 mm，模具部分高400 mm，外径为100 mm，壁厚为2 mm。在工艺边界条件方面，建立壁厚为4 mm的筒状水浴模型(水托盘厚度也为4 mm)。水筒模型高度为450 mm，以保证能完全包裹住模具。水托盘和水筒的位移曲线按表1设置。

图9 几何模型示意图Fig.9 Geometries of the model

为同文献[12]的仿真结果进行对比，从文献[12]所设置的温度监测点中选择了4个更具有代表性的位置。其中，#3位于中心线上距地面300 mm处，#5位于中心线上距地面100 mm处，#8和#10在#3和#5基础上向右偏移25 mm。

在网格划分方面，通过网格无关性测试，将网格平均尺寸设置为3.5 mm，利用Viusal-mesh模块完成有限元网格的划分，总计生成108万个计算单元。最后，仿真计算截止条件设置为冷却到30 ℃即停止。

图10为4个监测点的温度仿真结果。从图10中可以看到，在熔融态(77～81 ℃)阶段，均出现了显著程度不一的温度平台期，这正是结晶潜热释放导致的。由于#3和#5处于轴心位置，自身热量释放较为缓慢，因此在结晶潜热作用下会产生较为明显的平台期，而#8和#10靠近外壁，自身热量释放本身较快，因此，结晶潜热对#8和#10的作用相对较弱。并且越靠近底部的，由于换热效率更大(越早地利用水浴冷却)，因此#10的温度平台期几乎无法辨识。在符合度方面，除了3#与10#在平台期结束区域稍有偏差外，基于本文方法的仿真结果与文献[12]的计算结果，符合度非常高。

图10 本文方法与传统方法温度仿真结果曲线Fig.10 Comparison of the simulated temperature by present method and the conventional method

取3个不同时刻的温度与固相率分布云图，如图11。在本算例中，固相界面总体以V字型，随液面升高而自下而上的回缩上升。在中前期，中心区域始终没有完全凝固，从而保证了补缩通道的通畅，符合文献[3,12]中关于顺序凝固工艺中固相界面上升特性的描述。

图11 3个时刻的温度与固相率分布云图Fig.11 Distribution of temperature and solid fraction at three different time points

从图11的仿真结果可以发现，液面上升速度略快于固相界面的最低位置，且由于没有对未入水部分采取保温措施，中心补缩通道较小，这在实际生产中可能会增加孤立区形成的风险，因此，有必要进一步对工艺参数与炸药物性材料的匹配开展进一步深入研究。

图12反映了不同工艺下缺陷预测结果。其中，完全自然冷却条件下，由于其“回”字形界面生长特性使得补缩通道提早关闭，易形成孤立区，在中心区域分布着大小不等的缩孔。在采取冒口保温条件下，中心区域的缩孔有所减少，但由于其固相界面本质上还是以类似“回”字形的方式上升，因此中心区域依然存在缩孔。在顺序凝固工艺中，通过控制界面自下向上顺序生长，很好地将缩孔缩松缺陷控制在冒口区域，符合2.1节顺序凝固工艺机理分析的相关描述。

图12 不同工艺下缩孔预测结果云图Fig.12 Predicted shrinkage cavity with different processes

为测试该方法中筒状水浴模型壁厚对仿真精度的影响，建立了5组壁厚分别为4 mm、8 mm、20 mm、15 mm以及30 mm的筒状水浴模型并进行了仿真计算，仿真结果如图13所示，从结果中提取了3#监测点的温度曲线。从图13中可以看到，除了凝固结束后的一段区域(约4 000 s至5 500 s)存在一定的偏差外，在其他大部分区域中，不同壁厚下得到的温度曲线几乎完全一致。从图13可以看到，在55 ℃左右，20 mm 的情况相对4 mm的情况，其温度曲线存在最大约360 s的延迟，仅占总冷却时间(约11 000 s)的3%左右。对于这种微小偏差，其可能的原因是：在该测试中，对不同壁厚的水浴模型采用了相同的平均网格尺寸，导致弹体模具及药浆的网格并不完全一致，从而在计算精度上产生了一定误差。对此，可通过对筒状水浴部分与弹体模具部分根据各自的几何尺寸独立划分网格予以消除。

图13 不同壁厚水浴模型3#点的温度仿真结果曲线Fig.13 Simulated temperature at 3# with different thickness of the water bath

总的来说，该测试表明水浴模型壁厚尺寸对仿真精度的影响是非常小的。在实际仿真中，可在满足ProCAST对网格尺寸与模型尺寸的相关限制条件下，尽可能建立相对较薄的水浴模型，从而保证该方法在计算量上相对传统工艺边界建模方法不产生过多的额外负担。以本文中的情况为例，在3.5 mm总体网格的条件下，总体计算单元为108万个。其中额外添加的筒状水浴部分约32万个计算单元，约占30%。

5 结论

1) 构建了一种基于移动边界的顺序凝固工艺边界建模方法。建模过程简明直接，自动化程度高，大幅度降低了工艺边界建模工作量。

2) 该方法能够实现对熔铸装药顺序凝固工艺的精确仿真计算，其温度场仿真结果与文献[12]的结果高度符合，缩孔缩松缺陷预测与工艺原理分析相一致。

3) 该方法中水浴模型的壁厚对计算精度的影响小，不同的壁厚水浴模型仿真结果存在的时间偏差最大约为3%。在实际仿真建模中，可通过减小模型壁厚来减少水浴模型的额外计算量。在4 mm壁厚条件下，筒状水浴部分计算量单元约占总体计算单元30%。