铸造工艺中传热解析参数的推定研究

为生产高品质的合金铸锭，必须控制铸造工艺中溶质和杂质的偏析，防止发生铸造缺陷。为控制偏析 和防止发生铸造缺陷，正确理解凝固过程极为重要，采用数值模拟的北京中惠普热解析仪传热解析传热解析是正确理解凝固过程的有 效方法。
传热解析一般采用以下公式进行解析。
式中，ρ表示密度，Cp表示比热，T 表示温度，k 表示导热系数，△H 表示潜热，fs表示固相率，qmold/alloy 表示金属熔液-铸模间的热流束，h 表示金属熔液-铸模间的传热系数，Talloy表示与铸模壁接触金属熔液的 温度，Tmold 表示与金属熔液接触铸模的温度。式（1）为铸模和金属熔液中的热传导方程式，式（2）为 铸模和金属熔液边界中的热传导方程式。

为进行适当的北京中惠普热解析仪传热解析，必须输入适合这些方程式用参数的值，但这些输入值的精度不够高，因此， 无法进行适宜的传热解析。与密度和比热相比，导热系数k 和传热系数h 的精度存在着更多的问题。 在合金铸造过程中，固-液界面会产生溶质分配，在浓度场出现不均匀的同时，液相中会产生流动。 浓度场的不均匀性和流动，对凝固过程有很大的影响，因此，除了可以采用原有的式（1）和式（2）进 行解析外，还必须将溶质扩散方程式和流体运动方程式联立起来进行解析。但是，基于计算成本的考虑， 大部分的解析是只采用传热解析对凝固过程进行解析，而没有针对浓度场和流体的解析。在进行这种模 拟解析时，需要包括受浓度场的不均匀性和流动影响在内的表观导热系数k。对于这种表观值，不仅合金 铸造时需要，而且随着铸造条件的改变，这种表观值也会发生变化，因此，需要一种既精度高又简便的 计算方法，来求解 k 的表观值。为使 h 与温度的测定值和传热解析的结果一致，大部分是采用尝试法来 确定h。尝试法一般得不到充分的精度。尤其是在考虑到h 会随时间变化的情况下，必须事先对表示这种 时间变化的近似式进行假设，但遗憾的是该近似式有时不能正确反映h 的时间变化。

因此，需要一种能高精度求解k 和h 的方法。这种方法最好是一种不需要事先对h 的时间变化进行假 设，不依赖尝试法的方法。另外，最好这种方法不是仅限于特殊铸造条件的传热解析，还要具有可扩展 性的传热解析；不仅可简单应用于各种铸造工艺的传热解析，还可应用于在传热解析中，将浓度场和流 体的方程式联立的传热解析。

传热问题中的参数求解尝试，相当于所谓的传热解析的逆问题，以前就进行过各种各样的研究，并提出了各种 h 和 k 的推定方法。例如，有将灵敏度分析和聚合计算并用推定 h 的方法，有根据贝叶斯定理 推定正常状态下的 k 和 h 的方法，还有根据特殊试验测定的热态物性值采用遗传算法推定 k 的方法。这 些推定方法在各自的应用范围内，都是很好的推定方法，但为开发兼有上述扩展性和简便性的推定方法， 还需进一步研究。尤其是能够对 k 和与时间有关的 h 进行同时推定并不容易，希望开发新的推定方法。 作为这些参数的推定方法，本研究针对数据同化方法进行了研究。 数据同化方法就是将实测数据代入模拟解析中的一种科学实测数据的方法，它可用于提高模拟精度和 推定参数，现已广泛应用于大气海洋领域。本研究将数据同化方法应用于铸造工艺传热解析中的参数推 定，并验证了这些参数的有效性。但是，从应用于各种铸造工艺时的实际情况来看，有必要详细弄清以 下三项内容：①在单独推定 k 和 h 时，温度测定点（热电偶的位置）与推定精度的关系；②与时间有关 的h 的推定；③k 和与时间有关的h 的同时推定。但有关这三项内容的详细调查并未见到，因此，本研究 对这三项内容进行了详细调查，并调查了粒子滤波算法的有效性。

计算方法

1.1 传热解析和验证方法
为验证铸造工艺中数据同化的有效性，所有参数都应是试验所得，在进行铸造过程中的热分析后，根 据冷却曲线对 k 和 h 进行推定。本研究采用被称作“成对比较试验”的简便验证方法替代已准备好的试 验所得的已知参数。该方法采用事先设定的 k 和 h 值进行传热解析，然后把根据解析计算的冷却曲线， 视为实测数据对k 和h 进行推定。 通过把推定结果和事先假设的k 和h 值进行比较，可对精度进行验证。本研究传热解析用的铸模和金 属熔液的模式示于图1。图1（a）表示金属熔液和铸模的轮廓图，图 1（b）表示实际计算所采用的一维 计算面积。假设 y 轴方向的传热是可以忽略不计，就只对 x 轴方向进行一维解析。铸模内和金属熔液内 的热传导可采用差分法对式（1）进行离散，然后用焓法进行解析。金属熔液的中心部，即图1（b）的右 端可采用密勒定理的使用条件进行更新，铸模和金属熔液的边界温度可在式（2）的基础上进行更新。由 于铸模会与一定温度的大气接触，因此，图1（b）左端的边界温度，也可以采用热传导方程式进行求解。 另外，假设铸模为铸铁，金属熔液为Al-4mass%Cu 合金，可将表1 所示的参数作为输入值。在这些参数 中，金属熔液的导热系数k 和铸模-金属熔液间的传热系数h 是参数推定的对象。以下对成对比较试验中 实测数据的求解方法进行说明。 假设k=87Wm-1K-1，h=600Wm-2K-1时的传热解析结果示于图2。图中的实线、点线和虚线表示金属熔 液中的冷却曲线，分别为距离铸模34mm、16mm 和1mm 位置上的温度变化。本研究以1s 的间隔，从这 些冷却曲线中抽出温度值。实际铸造工艺中的温度测定会产生误差，由于考虑到这种测定误差，因此对 于计算出的冷却曲线给出了平均导热系数为0k，标准偏差按照σT的正态分布的噪声，并将该数据作为测 定值。对表 2 所示的三种值进行了研究，成对比较试验就是采用这种实测数据，对事先设定的真值进行 调查，即调查k=87Wm-1K-1，h=600Wm-2K-1是否可以推定。

1.2 基于数据同化方法的参数推定
本研究采用粒子滤波算法进行参数推定。粒子滤波算法是采用蒙特卡洛法，对贝叶斯定理的状态和参 数进行推定的方法。在这里，k 为已知数，以没有发生时间变化的 h 推定为例，对粒子滤波算法的顺序 进行说明。 首先，在开始推定时，由于不知道 h 的真值，因此，随机抽取了数个 h 值。图 3（a）表示抽取的试 样。纵轴为 h，橫轴为时间 t 的图，虚线表示真值。该例在 t=0 时，在 0-3000Wm-2K-1范围内，随机抽取 128 种值。在粒子滤波算法过程中，使用这些h 值进行传热解析，直至获得实测数据为止。也就是说，在 这种情况下，要在1s 内进行128 种传热模拟。对这些传热解析集进行了分析，每个解析都称为粒子，在 本例中处理了128 个粒子。所得的冷却曲线用点线示于图3（b）。本研究为便于观看，只列出128 条冷却 曲线中的20 条。这些曲线是在靠近铸模的金属熔液中测定的温度实测数据，计算结果也表示相同位置的 冷却曲线。在 t=1s 时，实测数据只表明热电偶的数量，与每个热电偶对应的解析结果有 128 种。此时， 采用粒子滤波算法进行了滤波操作，计算出各自模拟结果的“似然”值，并以该“似然”值为基础对 h值进行更新（再抽样）。所谓“似然”值是表示模拟结果与实测数据具有某种程度吻合的值。
在计算各粒子的“似然”值后，按照该“似然”值对h 进行再抽样。再抽样的目的在于毁弃“似然” 值较低的粒子，复制“似然”值高的粒子。再抽样后的h 值，可以做成图3（a）的t=1s 的线图。本研究 中粒子i=55、h=2000 的计算被放弃，粒子i=123、h=700 的计算被复制时，还包含了粒子i=123 的温度场 的计算结果，结果可复制为粒子i=55。但是，为防止粒子i=55 和i=123 成为完全相同的计算，对粒子i=55 的h 给出了噪声。本研究视这种噪声平均值为0，标准偏差按照σT的正态分布随机抽取。另外，对于σT 值的变化也进行了调查，但在本解析范围内，推定精度与σT值并没有很强的相关性。 本文对再抽样后的h 和温度场进行了 128 种的传热解析，获得了 t=2s 内的冷却曲线（如图3（b））。 因此，在对“似然”值进行计算后再抽样。这种计算一直进行到取得实测数据为止。结果如图 3（b）所 示，所有 h 值都会随这种计算而接近真值。但是，采用粒子滤波算法时，在粒子数少的情况下，有时由 于粒子会从推定途中，只向 h 窄的范围内集中，因此，有必要增加粒子数进行推定，但增加粒子数，意 味着会增大计算成本。因此，有必要设定适当的粒子数来解决这一问题。

2 结果

2.1 推定精度随热电偶位置的变化
首先，把导热系数k=87Wm-1K-1，传热系数h=600Wm-2K-1作为固定值，只对h 进行推定。假设温度 测定点只有一处，粒子数为128 个，可根据σT=1K 的实测数据对h 进行推定。将热电偶设置在金属熔液 中心（距离铸模壁34mm 的位置）时的结果示于图4（a）。虚线表示h 的真值，实线表示“似然”最高粒 子h 的值（极大“似然”推定值），点线表示h 的最大值和最小值。初期极大“似然”推定值h 会产生大 的振动，并在t=10s 左右时，聚合成真值。另外，最大和最小值在初期的振幅宽，但在t=15s 左右时会在 真值附近聚合。图 4（b）示出将热电偶放在距离铸模壁 1mm 的铸模内时的结果。与图 4（a）相比，极 大“似然”推定值h、最大值和最小值会在非常短的时间内，向真值附近聚合，表明在参数推定时，粒子 的变化会因热电偶设置位置的不同，而产生明显差异。 因此，对热电偶在铸模和金属熔液的不同位置设置时的推定精度进行了验证。在这里，把在t=15s 至 t=30s 的范围内每1s 所得的极大“似然”推定值的平均推定时间作为推定结果，尤其是在相同条件下进行 了 5 次推定，并把这些推定结果的平均值用曲线表示，用估计值表示极大“似然”推定值的最大值和最 小值。另外，虚线表示 h 的真值。横轴为根据与铸模壁的距离表示温度测定点的值。当 x 为正时，表示 在金属熔液内部测定的结果；为负时表示在铸模内部测定的结果。根据这些测定结果可知，把热电偶放 置在铸模内时，估计值与真值的误差大约在 10%以内，表明推定精度高。尤其是热电偶越靠近铸模壁， 推定精度越高。另外，即使把热电偶放置在金属熔液内时，如果热电偶的放置位置距离铸模壁大约5mm， 也能获得较高的推定精度。 以上推定结果表明，采用数据同化可以分别对h 和k 进行推定。尤其是在5 次推定的估计值大的情况 下，其平均值也会偏移真值。也就是说，即使对真值不清楚，根据数次参数的推定，也能对估计值的大 小进行评价，由此可以对温度测定点与推定精度的关系进行预测。因此，在实际铸造工艺中，确定热电 偶的放置位置前进行这种解析，就可事先验证测定点与推定精度的关系，这在实际应用上是有益的。

2.2 与时间有关的传热系数的推定
对h 没有随时间变化的情况进行了假设，但在大多数情况下，h 是与时间有关的量。通常，在求h 的 时间变化时，要事先设定其近似式，确定该系数推定值的估计方法。此时，要充分注意该近似式能否能 体现 h 的时间变化。在推定 h 时，采用铸模壁附近的冷却曲线，就可在非常短的时间内进行推定。根据 这一结果，对于与时间有关的h，也希望极大“似然”推定值的h，能在各时点快速向真值聚合。也就是 说，如果按照数据同化进行推定，希望在即使不导入任何假设的情况下，也能对h 的时间变化进行推定。 因此，本研究调查了极大“似然”推定值 h 的时间变化能否记述真值的时间变化。另外，还希望能根据 各时点的“似然”值，求出h 随时间变化的期待值，但本研究采用的是更加简便的极大“似然”推定值。 铸造时的传热系数常常会受凝固收缩等的影响而随时间变化减小。例如，以α和 m 作为常数，可以 采用 h（t）=αt -m的式子。因此，在本研究中假设真值为 h（t）=3000t-0.5Wm-2K-1。使用该 h（t）的式子和 k=87Wm-1K-1再次进行传热解析，把由此获得的冷却曲线作为实测数据。在这里，只列出σT=1k 的结 果。设定热电偶的测定点为1 个，粒子数为128 个，σh是原值的10%，如表2 所示。由此可知，推定结 果与真值非常一致。也就是说，采用本方法时，只要根据一条冷却曲线，就可高精度推定出随时间变化 的h。尤其是它不需事先对h 的时间变化进行假设，这是非常重要的。 关于30s 的冷却，在相同条件下进行了5 次推定，在把热电偶放置在铸模内的情况下，距离铸模和金 属熔液边界线3mm 左右时的推定误差为10%左右，可知推定精度较高。另外，即使在把热电偶放置在金 属熔液内，如果热电偶距离铸模壁非常近为1mm 左右时，推定误差也在10%左右。

2.3 导热系数和传热系数的同时推定
根据前面的推定结果可知，在对k 和h 单独推定时，用一个热电偶就能进行高精度推定。但是，在推 定 k 时，热电偶应放置在金属熔液内；在推定 h 时，热电偶应放置在铸模内或靠近铸模的金属熔液内。 在同时推定k 和h（t）时，要在铸模和金属熔液内设置热电偶。 本研究对热电偶设置的两种方法进行了试验。第一种方法是在铸模和金属熔液内，分别设置一只热电 偶进行推定，共计设置了两只热电偶。第二种方法是在铸模内设置两只热电偶，在金属熔液内设置一只 热电偶进行推定，共计设置了三只热电偶。热电偶设置如图 5 所示。在共计设置两只热电偶的推定中， 铸模中的热电偶位置和金属熔液中的热电偶位置都与铸模壁的距离相同为x。即使在设置三只热电偶的推 定中，热电偶在铸模中和金属熔液中的设置位置与铸模壁的距离也相同。但是，设置在铸模内的第三只 热电偶与另一铸模内热电偶的设置距离为 10mm。k 和 h（t）的真值分别为 k=87Wm-1K-1，h（t） =3000t-0.5Wm-2K-1。由于在粒子滤波算法中，必须对k 和h（t）两者进行抽样，因此必须考虑的粒子比以 往都多。因此，设定粒子数为1282 。 推定结果表明，适当调整热电偶的位置，就可同时并高精度推定k 和h（t）。关于h（t）的推定，如 果热电偶的设置位置靠近铸模壁，就能推测出在各时点上 h（t）的真值。在热电偶与边界距离小的情况 下（x≤4），不论是两只热电偶还是三只热电偶的推定，推定误差都在10%以内。另一方面，当距离大于 5mm 时，热电偶在两只和三只的情况下出现的差别较小。由此可知，x 的值无论在何种情况下，使用三 只热电偶的推定精度更高。 另外，k 的推定结果与h（t）的推定误差趋势一样，在热电偶与铸模壁距离小的情况下，能得到高精 度的推定结果，它与热电偶的设置数量无关。但是，在热电偶为两只的情况下，当热电偶与铸模壁的距 离大于3mm 时，推定精度会急剧下降。而且，k 的推定结果也表明，使用三只热电偶可以抑制推定精度 的下降。

结语

本研究把作为数据同化方法的粒子滤波算法应用于铸造工艺中导热系数和传热系数的推定，尤其是① 热电偶位置对推定精度影响的验证；②随时间变化的传热系数的推定；③导热系数和与时间有关的传热 系数的同时推定，并就粒子滤波算法的有效性进行了验证。结果表明，在单独推定导热系数的情况下， 把热电偶设置在金属熔液中,能进行高精度推定；在单独推定传热系数的情况下，把热电偶设置在铸模中， 能进行高精度推定。而且，利用传热系数的“似然”，无需事先对函数形式进行假设，就可对随时间变化 的传热系数进行推定。尤其是在铸模内，只设置一只热电偶的情况下，能进行高精度推定。另外，还可 对导热系数和传热系数进行同时推定。由此表明，粒子滤波算法是铸造工艺传热解析时参数推定的有效 方法。 粒子滤波算法的精度和计算机成本会随粒子数的增加而增加。本研究在对上述①和②的推定时，使用 的粒子数非常少为 128 个，但两者的推定精度都很高。在进行上述③的同时推定时，使用 1282 的粒子数 就能进行高精度推定。无论何种情况，此次传热模拟本身的计算成本低，因此，能在极短的时间内进行 推定。 近年来，以GPU 为代表的并列计算技术取得发展。粒子滤波算法是接近并列计算的方法，因此，采 用并列计算，还能进行更加复杂传热模拟时参数的推定。