注塑成型Hele-Shaw流动模拟中热对流的异步长求解

王超房，黄　明，石宪章，申长雨，赵振峰



注塑成型Hele-Shaw流动模拟中热对流的异步长求解

王超房，黄明，石宪章，申长雨，赵振峰

（郑州大学橡塑模具国家工程研究中心，河南 郑州 450002）

摘要:在Hele-Shaw流动数值模拟中，速度是压力的后处理结果。如果是点浇口，则浇口附近速度会随单元尺寸缩小而趋于无穷大，导致能量方程作为一个整体求解时，时间步长必须非常小，否则会产生很大误差；而根据热对流物理意义分步求解，则需追踪当前物质在上一时刻位置，当单元速度很高、逆向搜索需穿透多个单元时，搜索可能会失败。鉴于此，基于分步求解法，研究提出一种变长度子时间步长方法处理对流项，确保搜索路径局限在当前单元内，并采用二分法确定子时间步数量，使算法简洁有效。算例表明，该方法在保证计算精度和求解稳定性的同时，可以明显减少计算时间，提高计算效率。

关键词:热对流；传热；Hele-Shaw；注塑成型；数值模拟；聚合物加工

引 言

目前，注塑成型 CAE分析已经成为模具制造前不可或缺的环节，其将模具设计和工艺参数设置建立在科学的分析基础之上，对缩短生产周期、提高产品质量具有重要指导作用[1-5]。然而，相关的CAE计算方法并非完美无缺，许多问题仍有待于深入研究[6]，其中 Hele-Shaw流动中的热对流求解就是其中之一。

在注塑成型模拟中，由于熔体黏度严重依赖温度，因此温度场的计算至关重要。一般而言，温度场的控制方程包含三要素:对流、传导和能量转换。所谓热对流，就是承载一定温度的物质由于自身运动而导致的温度场变化，因此与速度场密不可分。在处理热对流的数值方法中，大致有两大类:一类是直接求解，即热对流和热传导作为一个整体同时求解；另一类是分步求解，即热对流从控制方程中分离出来单独求解[7-8]。前者直接面对控制方程，可以认为是一个纯数学问题，不用考虑其物理意义，基于该方法，文献[9]开展了注塑成型Hele-Shaw流动温度场的模拟，而文献[10-13]则模拟了三维流动过程温度场，为了抑制数值震荡，上述研究均采用隐式格式离散方程，并使用上风类方法[14]处理热对流项，但上风法仅考虑与节点相连的一层单元中的上游单元的贡献，如计算中时间步长选择不当，会导致节点速度与时间步长的乘积超越该节点周围一层单元的空间尺度，造成求解误差被无限放大，因此该时间步必须小于某个特定值，而特定值是由相关上游单元尺寸和节点速度共同决定的。针对上述问题，文献[15]基于分步法，分别求解对流方程和热传导方程，实现了三维注塑充填温度场数值模拟。

对于大部分注塑件，其几何造型为薄壳结构，更适合用Hele-Shaw模型描述[16-18]，而该模型的直接求解对象是压力，速度是通过节点压力计算得到的。对于点浇口，由于浇口通常只能设置在一个有限元节点上，这导致在浇口附近，速度值会随单元尺寸的缩小而趋于无穷大。这给温度场计算带来极大的困扰:除非采用极小的时间步，否则第一类求解方法必然产生很大误差，而采用小时间步则会极大延长计算时间。第二类求解方法则是根据热对流的物理意义，将其从控制方程中独立出来，寻求当前时刻的物质在上一时间步的位置，从而得到热对流的解，然后再计算热传导项。这种方法能有效避免第一类解法产生的问题，但算法相对复杂。如果某个单元的速度很高，则逆向搜索的路径可能会穿透多个单元，使得计算复杂化，并且在筋板一类的地方逆向搜索可能失败。

针对上述问题，本文提出一个基于第二类解法的异步长算法。在逆向搜索物质点运动轨迹时，对于速度高的质点，采用更多而更短的时间步，反之用较少同时较长的时间步，即将时间步分为若干子时间步，确保搜索路径局限在一个单元内，采用二分法确定子时间步数量，使得此类解法更简洁有效。并用算例验证了该算法的有效性和可靠性。

1　Hele-Shaw流动控制方程

注塑充模过程属于黏性不可压缩流体的流动，基于Hele-Shaw模型假设[19-21]，熔体流动和传热控制方程可表述为

式中，v、p、η、ρ、cp、t、T、k和γ˙分别为速度、压力、黏度、密度、比热容、时间、温度、热传导系数和剪切速率。对（2）式沿厚度方向(z向)积分，可以得到速度的解析式

式中，h+为型腔半厚度。将式（4）代入式（1），并在厚度方向积分，即可得到充模阶段计算压力的泊松方程

2　热对流项（温度场）求解

2.1 温度场分步求解方法

利用熔体在注塑充模过程为层流且为对流占优的特点，Hele-Shaw流动模拟将熔体沿厚度方向分层，在不考虑熔体前沿的情况下，熔体只在同层内流动。与此相对应地，计算温度时，认为对流只发生在同一层内，同时忽略同层内的热传导，而热传导只沿厚度方向进行。对于热传导，时间步长应该满足步长与速度乘积在速度反方向上不超越该单元边界的条件。否则，计算过程将是不稳定的（对于前差分），或产生较大误差（对于后差分）。对于点浇口，这个步长限制会产生一个比较严重的问题，即当单元尺寸很小时，在点浇口周围会产生很大的流动速度，从而时间步长必须非常小才能得到可靠的解。如果所有的点、所有的层都采用统一的最小时间步，将式（3）作为一个整体进行后差分求解，必然会非常耗时。针对以上问题，本文建议采用分步求解法，即根据能量方程的物理意义分别求解以下方程

即在给定的时间步内:① 不考虑热对流和热传导，假定熔体处于静止的条件下，求解方程（6）得到剪切热引起的温度场变化；② 不考虑热传导，只计算温度场随熔体运动产生的变化，即求解对流方程（7）；③ 然后再假定熔体处于静止状态，只计算热传导，即只对每个节点（或单元）分别求解一个一维热传导方程。上述温度场的计算采用与前面压力场计算相同的有限元网格，该网格沿中面划分，即中面网格。同时，对于每一个单元还要按一定的准则在厚度方向分层，如图1所示。同一单元的不同层以及不同单元的相同层可以有不同的层厚，并假定熔体除了前沿之外，只能在同一层内流动。温度定义在单元的节点和各层的交界处，且速度、比热容、剪切速率和黏度由压力场计算得到。

图1　三角形单元及分层Fig.1 Triangular element and layered

2.2 对流方程异步长求解

对于给定的时间步长Δt，很容易得到式（6）的解

式（9）适用于每一层。

图2　热对流的物理意义Fig.2 Physical meaning of thermal convection

本节重点讨论式（7）的求解。在一个速度为v的局部流场中，在时间步 tΔ内，一个微粒由A点移动到B点，如图2所示，这样在计算B点的当前温度时，如果不考虑其他因素，就需要找到该点的物质在上一时刻所在的位置A，而A点在上一时刻的温度就是B点当前的温度，即

在Hele-Shaw模型中，压力是线性插值，通过式（4）只能计算得到单元中心速度。考虑到对于中面网格，某些转角处的节点速度无法定义，因此本文算法将温度定义在单元中心，再通过将单元温度加权平均到节点，然后对节点温度进行单元插值形成温度场。考虑到流体运动是层流，求解热对流过程主要就是在同一层内沿速度的逆方向追踪当前相关点B在上一时间步的位置A，如图3所示，然而在算法上该过程有两个缺点:一是在速度很大的区域，逆向搜索会穿越很多单元，计算很繁复；二是当搜索穿越单元遇到筋板或熔接线等地方时，搜索方向会遇到歧路。

图3　逆向追踪Fig.3 Tracking in reverse direction of velocity

图4　单元内逆向追踪Fig.4 Reverse tracking inside an element

针对上述问题，本文提出一种变长度子时间步长求解方法，其基本思想为:将搜索路径局限在当前单元内，如果给定的时间步长与速度之积要突破该单元，则将给定的时间步再划分为若干子时间步，确保在子时间步内搜索路径不会突破该单元。不失一般性，以图4为例，在计算B点温度时，根据上面准则将时间步划分为N个等长的子时间步，而在每一个子时间步内的逆向搜索，都会从B点逆向行进到a点。在第n+1个子时间步，B点处的温度为

其中，Li、Lj、Lk为 a点的面积坐标，和由式（12）计算得到

但上述做法对注塑成型而言效率很低，因为在成型过程中不同区域的速度差异极大，有的点可能需要划分上百个子时间步，而有的点可能根本不需要。都采用统一的子时间步，会导致计算时间不必要的浪费。为解决这个问题，本文建议对不同的点采用不同的子时间步。众所周知，在注塑流动模拟中，总的时间步长Δt并不是由温度计算自己规定的，而是取决于流动前沿节点的充填时间。如果对流计算中某个点能够允许的最大时间步小于Δt，则必须设定自己的子时间步

其中，d为待求温度点到该点速度反方向与单元边界交汇处的距离，N为一整数。仍以图4为例，假定编号1的点（B点）需划分的子时间步长为Δt/4，编号2、3、5、6需要的子时间步长为Δt/2，编号4的点不需划分，即Δt。初始时刻t0的温度在各节点已知，在经第一次计算后，得到了、和。当根据式（12）计算节点i温度时，需知道2、 3、4、5、6各点在(t0+Δt/4)时刻温度，该温度通过对时间的插值求得

上述计算必须要遵守次序，即必须先处理时间滞后的点。每计算完一个点的一个子时间步，就搜索下一个时间滞后的点，虽然算法不复杂，但会降低计算效率。为解决此问题，本文采用二分法确定子时间步，即子时间步数满足 2N0，相应的子时间步长为Δt/2N0，然后根据N0的大小对节点分组。实践表明，N0通常不会超过10。这样计算次序就变成先处理时间滞后的组，而组数是十分有限的。

通过上述算法得到某个时间步的热对流解后，还要求解式（8），即对各个节点在厚度方向求解一次时域为(t,t0+Δt)的一维热传导方程。在厚度方向，将各层作为一个一维单元并进行有限元离散，对时间后差分，则得到标准的三对角有限元方程

其中，K为刚度阵，M为质量阵。处理完边界项后，应用追赶法或三角分解很易得到方程的解，具体过程参见文献[22]，这里不多赘述。

3　算 例

基于上述理论和算法，本文改进了自主研发“注射成型模拟仿真系统”（Z-Mold）的中面流动传热模块，并以平板和某型号航天面窗为例开展模拟分析。需指出，前两个算例只是为了突出本文算法的可行性，并不意味着一定符合工程实际设定（如浇口位置和注射时间）。此外，在极短时间内对模具内高速流动的熔体进行准确的温度在线测量，目前还有一定困难，所以只定性地讨论其合理性。

算例1 长、宽和厚度分别为200 mm、100 mm 和2 mm的平板，网格尺寸4 mm，点浇口注射，浇口设在长边中点，如图5所示。材料为PP（Adflex KS-221P），注射时间1 s，熔体温度230℃，模具温度50℃，由流率到压力（V/P）控制设置为充填体积达到98%。

图5　算例1的网格F i g . 5 M e s h o f E x a m p l e 1

图6给出了本文方法（Z-Mold）模拟的体平均温度。对于这样的点浇口注塑，在开始一段时间熔体前沿总是呈半圆形状，热的熔体与模壁接触时会在型腔表面形成冷凝层，熔体像“三明治”一样被夹在冷凝层之间，熔体速度在厚度方向呈抛物线分布，即中心层速度最高，而且受模具传热的影响十分有限。同时由于剪切热的存在，这样会将一个略高于浇口温度熔体带到前沿，图6（a）清楚显示了这一过程。当熔体流动遇到边界时，相应的部分流速会降低或停止，此时热对流的作用降低或消失，热传导起主要作用，熔体温度快速下降。热对流在仍然流动的部分继续以前的作用，熔体温度在前沿保持最高，在前沿后方温度下降也比静止部分缓慢，图6（b）清楚地展现了这一过程，说明本文算法给出的结果是合理的。

图6　算例1体平均温度Fig.6 Bulk temperature of Example 1

算例2 一长500 mm、宽200 mm、厚2 mm平板，点浇口注射，浇口设在短边中点，网格剖分如图7所示。材料为ABS，熔体温度200℃，注射时间1 s，模具温度40℃。本算例中，人为地构造一个高剪切场，用以考察本文算法对计算效率的改进。

图7　算例2的网格Fig.7 Mesh of Example 2

如前所述，在注塑数值模拟中，每个时间步长是由充填程序决定的，如果此步长大于对流计算的步长要求[式(13)]，则需要在对流计算中选取合适的N减小步长。在本算例中，取决于每个步长的长短，式(13)中的N为2～32不等。

（1）如果将N在整个计算中取为定值1～8，计算会产生溢出。取为16时，计算可以完成。这表明如果在少量的时间段不满足式(13)，可以进行计算。但是，这取决于具体问题和计算机系统，是一个很难定量给出的准则，为了计算的稳定可靠，对流计算的时间步长应该满足式(13)。

（2）在每一个时间步中比较各个单元，选取满足式(13)的最大的 N，即以最小的Δt/N作为本步计算的子时间步长。在充填完成时，体平均温度和压力分布如图8所示。截止充填结束，总用时369 s，对流计算用时228 s。

（3）采用本文的异步长算法，充填完成时体平均温度和压力分布如图9所示。截止充填结束，总用时206 s，对流计算用时65 s。

图8　使用统一子时间步长充填完成时计算结果Fig.8 Computing results of Example 2 at filling end time by using uniform sub time step

图9　使用异步长充填完成时计算结果Fig.9 Computing results of Example 2 at filling end time by using variable time step

由于存在高剪切，熔体前沿温度要远远高于熔体入口温度，计算结果充分反映了这一点，与文献[5]的论述是一致的。采用异步长计算结果与采用统一的最小时间步长结果相仿，但计算时间大大缩减，说明本方法在保证计算精度和求解稳定性的同时，可以明显提高计算效率。

算例 3 本文方法与其他计算模块一起，成功分析、指导了某型号航天面窗的模具设计和工艺设定，其形状类似于半径145 mm的1/4球壳，网格尺寸10 mm，点浇口注射，如图10所示。材料为光学级PC，注射时间4 s，熔体温度268℃，模具温度90℃，由流率到压力（V/P）控制设置为充填体积达到98%。

图10　面窗中面网格Fig.10 Mid-plane mesh of transparent window

图11　充填结束时刻体平均温度Fig.11 Bulk temperature at filling end time

图11为模拟得到的体平均温度分布。图中流动前沿温度最高，而浇口附近略低于熔体温度。实际上，正如文献[5]所论述的，注塑过程中流动前沿呈泉涌状，其平均温度在大多情况下是最高的，略高于入口温度；而浇口附近，由于受模壁冷却时间最长，同时不断有新的熔体进入，该处平均温度并不是最高。算例表明本文方法对于计算热对流是有效、可靠的。

4　结 论

针对当前CAE在Hele-Shaw流动热对流求解过程中存在的问题，本文开展了以下研究工作。

（1）采用分步法求解能量方程剪切项、对流项和热传导项，提出了采用变长度子时间步长方法处理对流项，确保搜索路径局限在当前单元内，解决了逆向搜索路径穿透多个单元时搜索可能失败的问题，并采用二分法确定子时间步数量，使得算法简洁有效。

（2）基于本文提出的算法，改进了自主研发“注射成型模拟仿真系统”（Z-Mold）的中面流动传热模块。算例表明，Z-Mold模拟的体平均温度在流动方向上有明显的热对流痕迹，剪切速率高、流速大的位置并未出现热量聚集。在采用统一时间步计算对流时，如时间步大，不能保证计算的稳定（前差分）或精度（后差分），而时间步太小，又会耗时，本文方法有效地解决了该问题。这对于构造复杂、计算耗时的注塑模仿真，有着重要意义。

符 号 说 明

cp——比热容，J·kg-1·K-1

d ——待求温度点到该点速度反方向与单元边界交汇处的距离，m

h——型腔半厚度，m

k——热传导系数，w·m-1·K-1

L ——单元面积坐标

N ——子时间步数

p ——压力，Pa

T ——温度，℃

t ——时间，s

Δt ——时间步长，s

v ——速度，m·s-1

x, y, z ——笛卡儿坐标，m

γ˙ ——剪切速率，s-1

δ ——权值

η ——黏度，Pa·s

ρ ——密度，kg·m-3

上角标

n——子时间步

下角标

i, j, k——单元节点编号

References

[1] MATIN I, HADZISTEVIC M, HODOLIC J, et al. A CAD/CAE-integrated injection mold design system for plastic products[J]. Int. J. Adv. Manuf. Technol., 2012, 63: 595-607.

[2] 黄明, 石宪章, 刘春太, 等. 基于统一网格的塑件成型与模具结构一体化分析[J]. 化工学报, 2012, 63(8): 2617-2622.

HUANG M, SHI X Z, LIU C T, et al. Integrated analysis of part molding and mold structural mechanics based on identical mesh[J]. CIESC Journal, 2012, 63(8): 2617-2622.

[3] YANG B, OUYANG J, LIU C T, et al. Simulation of non-isothermal injection molding for a non-Newtonian fluid by level set method[J]. Chin. J. Chem. Eng., 2010, 18 (4): 600-608.

[4] ZHOU H M, SHI S X, MA B. A virtual injection molding system based on numerical simulation[J]. Int. J. Adv. Manuf. Technol., 2009, 40: 297-306.

[5] 王利霞. 基于数值模拟的注塑成型工艺优化及制品质量控制研究[D].郑州: 郑州大学, 2004.

WANG L X. Process optimization and part quality control for plastic injection molding based on numerical simulation[D]. Zhengzhou: Zhengzhou University, 2004.

[6] WANG C F, HUANG M, SHEN C Y, et al. Warpage prediction of the injection-molded strip-like plastic parts[J]. Chin. J. Chem. Eng., 2015, doi: 10.1016/j.cjche.2015.12.012.

[7] DEMKOWICZ L, ODEN J T, RACHOWICZ W. A new finite element method for solving compressible navier-stokes equations based on an operator splitting method and h-p adaptivity[J]. Computer Methods in Applied Mechanics and Engineering, 1990, 84(3): 275-326.

[8] MINEV P, VABISHCHEVICH P N. An operator- splitting scheme for the stream function-vorticity formulation of the unsteady Navier–Stokes equations[J]. Journal of Computational and Applied Mathematics, 2016, 293: 147-163.

[9] 王利霞, 申长雨, 陈静波, 等. 注塑成型充模过程的温度场计算[J].计算力学学报, 2002, 19(2): 173-178.

WANG L X, SHEN C Y, CHEN J B, et al. Calculation of the temperature field in filling stage of injection molding[J]. Chinese Journal of Computational Mechanics, 2002, 19(2): 173-178.

[10] ZHOU H M, YAN B, ZHANG Y. 3D filling simulation of injection molding based on the PG method[J]. Journal of Materials Processing Technology, 2009, 204: 475-480.

[11] 李阳, 严波, 赵朋, 等. GLS/GGLS/SUPG 在三维注射成形充填模拟中的应用[J]. 化工学报, 2010, 61(2): 510-515.

LI Y, YAN B, ZHAO P, et al. GLS/GGLS/SUPG method in 3D numeical simulation of filling stage of plastic injection molding [J]. CIESC Journal, 2010, 61(2): 510-515.

[12] 曹伟, 王蕊, 申长雨. 注塑成型三维流动的数值模拟 [J]. 计算力学学报, 2005, 22(6): 792-795.

CAO W, WANG R, SHEN C Y. Numerical simulation for 3D flow in injection molding[J]. Chinese Journal of Computational Mechanics, 2005, 22(6): 792-795.

[13] 耿铁, 李德群, 周华民, 等. 注塑充模过程中温度场的全三维数值模拟[J]. 化工学报, 2010, 56(9): 1612-1617.

GENG T, LI D Q, ZHOU H M, et al. Simulation of temperature field in filling stage of injection molding using three-dimensional model[J]. CIESC Journal, 2010, 56(9): 1612-1617.

[14] BROOKS A N, HUGHES T J R. Streamline upwind/ Petrov-Galerkin formations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[J]. Computer Methods in Applied Mechanics and Engineering, 1981, 32(1/2/3): 199-256.

[15] 刘永志, 赵振峰, 申长雨. 基于分步法的三维注塑充填温度场数值模拟[J].化工学报, 2011, 62(5): 1455-1459.

LIU Y Z, ZHAO Z F, SHEN C Y. 3D numerical simulation based on operator splitting method for temperature of plastic injection molding in filling stage[J]. CIESC Journal, 2011, 62(5): 1455-1459.

[16] 申长雨.注射成型模拟及模具优化设计理论与方法[M].北京: 科学出版社, 2009.

SHEN C Y. Injection Molding Simulation and Mold Optimization Design: Theory and Methods[M]. Beijing: Science Press, 2009.

[17] 刘永志, 赵振峰, 马兰, 等. 注塑充填过程变厚度截面速度压力场数值模拟[J].化工学报, 2009, 60(7): 1818-1822.

LIU Y Z, ZHAO Z F, MA L, et al. Numerical simulation for velocity and pressure in filling stage of injection molded parts with uneven thickness[J]. CIESC Journal, 2009, 60(7): 1818-1822.

[18] 柳和生, 匡唐清, 周国发, 等. 基于粘弹本构模型的注射成型充填过程模拟研究[J]. 机械科学与技术, 2008, 27(5): 579-582.

LIU H S, KUANG T Q, ZHOU G F, et al. Numerical simulation of mold filling in injection molding of viscoelastic polymers[J]. Mechanical Science and Technology for Aerospace Engineering, 2008, 27(5): 579-582.

[19] 林建忠, 阮晓东, 陈邦国, 等. 流体力学[M]. 北京: 清华大学出版社, 2005: 85-86.

LIN J Z, RUAN X D, CHEN B G, et al. Fluid Mchanics[M]. Beijing: Tsinghua University Press, 2005: 85-86.

[20] HIEBER C A, SHEN S F. A finite element/finite difference simulation of the injection molding filling process[J]. J. Non-Newtonian Fluid Mech., 1980, 7: 1-32.

[21] 刘春太, 申长雨, 陈静波, 等. 注射模充模流动和传热过程的理论与算法[J]. 高分子材料科学与工程, 2002, 18(6): 101-105.

LIU C T, SHEN C Y, CHEN J B, et al. Theory and algorithm of polymer melt flow and heat-transfer in injection-molding filling[J]. Polymer Materials Science and Engineering, 2002, 18(6): 101-105.

[22] 申长雨, 郭恒亚, 赵振峰, 等. 模腔表面平均温度边界对注塑仿真的影响[J]. 郑州大学学报(工学版), 2005, 26(1): 9-12.

SHEN C Y, GUO H Y, ZHAO Z F, et al. Influence of the average temperature boundary condition on mold cavity surface to injection molding simulation[J]. Journal of Zhengzhou University (Engineering Science), 2005, 26(1): 9-12.

2015-12-14收到初稿，2016-03-19收到修改稿。

联系人:黄明。第一作者:王超房(1982—)，男，博士研究生。

Received date: 2015-12-14.

中图分类号:TQ 320；TK 124

文献标志码:A

文章编号:0438—1157（2016）07—3047—08

DOI:10.11949/j.issn.0438-1157.20151900

基金项目:国家自然科学基金重点项目（11432003); 河南省高等学校重点科研项目(15A430009, 15A430045)。

Corresponding author:HUANG Ming, huangming@zzu.edu.cn supported by the Key Program of National Natural Science Foundation of China (11432003) and the Key Research Project for Henan Universities (15A430009, 15A430045).

Thermal convection calculation with variable time step in Hele-Shaw flow simulation of injection molding

WANG Chaofang, HUANG Ming, SHI Xianzhang, SHEN Changyu, ZHAO Zhenfeng
(National Engineering Research Center for Advanced Polymer Processing Technology, Zhengzhou University, Zhengzhou 450002, Henan, China)

Abstract:In Hele-Shaw flow simulation, the directly solved variable is pressure and the velocity is only the post-treating result of pressure. Around the injection gate, the velocity may be very high along with reducing elemental size. This means that when the energy conservation equation is solved as a whole, the time step must be very short, otherwise, the error in the heat convection is unavoidable. The above problem can be overcome by using the operator-splitting method, in which the material at the current computing points needs to be tracked back to its position in the last time-step. However, this may lead to a new difficulty. If the elemental velocity is very high, the tracking needs to pass through a few elements and the reverse tracking may fail. To solve the problem, a new algorithm named sub time step with variable size was suggested to deal with the thermal convection in this paper, in which the sub time step using dichotomy was introduced that bounded the tracking path in a certain element. The new algorithm made the computation more simple and effective. The numerical examples showed that the new method had same accuracy as one using uniform small time step and high solving stability, but calculating time was dramatically reduced.

Key words:thermal convection; heat transfer; Hele-Shaw; injection molding; numerical simulation; polymer processing