基于可压缩两相流模型的注射成型填充模拟

2024-03-11王振海贾伟董勤喜于行

工程塑料应用 2024年2期

王振海，贾伟，董勤喜，于行

[1.海南大学土木建筑工程学院，海口 570228； 2.东玺技术(山东)有限公司，济南 250003]

注射成型是批量生产塑料零件最流行和最有用的聚合物加工方法之一，已广泛应用于电子、汽车、医疗技术、航空航天等许多不同领域[1]。填充过程是注塑过程中最重要的阶段，该阶段可能会出现产品的许多缺陷，如填充不足、气穴、熔接线、流痕等。正确选择加工参数对于获得高质量的塑料零件至关重要。注塑填充过程数值模拟主要基于动力学和传热学[2]的数学模型，在计算机上以高效、稳定的数值方法进行求解分析，从而预测熔体在型腔中流动过程和成型各个阶段的物理量变化，如速度、温度、压力、剪切速率[3]和流动诱导应力等，这为优化模具设计和控制成型过程提供了科学合理的理论依据。

注塑填充阶段数学建模的最初发展仅限于一维情况，然后发展到Hele-Shaw (GHS)近似，该近似忽略了沿厚度方向的压力变化，并假设每个位置都有一个完全发展的速度分布[4]。De Miranda等[5]采用广义GHS近似法和全三维牛顿(GNF)法预测注射成型中熔融聚合物流动的特征，验证了GNF方法准确性更高。如今，塑料零件变得越来越复杂，这意味着厚度起着不容忽视的重要作用，GHS近似法的有效性仅限于简单的几何形状条件，为了克服这个问题，需要使用三维建模来模拟注射成型。Hétu等[6]基于标准Galerkin公式提出了一个三维有限元模型，速度和压力由广义斯托克斯方程控制。几年后，Ilinca等[7]用该方法解决气体辅助注射成型和共同注射成型过程，通过模拟C型板的填充过程说明了三维模拟比2.5维的优势。但三维有限元法在处理大形变、材料非线性或者边界条件的变化时，可能遇到数值不稳定性的问题，导致收敛困难以及计算结果的不可靠性。光滑粒子流体动力学(SPH)[8-9]方法是近年来逐渐发展的一种数值模拟方法，该方法不需要规则的网格，可以在自适应粒子分布的情况下进行模拟，具有适用于任意形状和拓扑的流体区域的优点。但该方法粒子是无序排列的，且每个粒子都需要与其他所有粒子进行相互作用的计算，这导致计算复杂度较高，对于稀疏流体区域的模拟效果不佳。有限体积法(FVM)[10]可以将复杂的计算领域划分为小的控制体积进行离散化，并在每个体积中求解物理守恒方程。这种离散化结构使得有限体积法具有高效稳定的特点，可以处理大规模计算。FVM不仅可以较准确地模拟注塑过程中的流动、热传导、固化等物理现象，还能够提供详细的场量分布信息，如温度、压力、速度等，以及有效处理复杂的物理边界条件，如模具壁面、冷却通道等。

随着计算机仿真模拟技术的飞速发展，三维实体模型在注射成型流动模拟中也已日趋成熟，并与中面模型一起处于主导地位，主流的模流分析软件如Autodesk Moldflow[11-14]，Moldex3D[15-16]，ZMold[17-18]和HsCAE[19]均已开发出完善的实体流动模块。然而，三维实体模拟技术也存在自身问题。其控制方程求解的复杂度远高于中面模拟技术和表面模拟技术，计算量也非常庞大，这是制约其推广应用的一个重要因素。另外，目前较为成熟的注塑商业软件均由境外公司开发，而国内的相关研究大多用于科研目的，模型考虑因素与实际工业生产存在一定差异。因此，有必要针对该课题进行进一步研究。

在现有的三维实体注射成型计算模型中，由于注射成型加工过程中涉及到高温熔融塑料在模具中的填充、冷却和固化等过程，气体往往不可忽略，因此气液两相流模型[20]被应用于注射成型模拟，其是将熔体和空气这两种不同性质的流体，采用同一套控制方程进行计算，用自由界面追踪方法捕捉两相界面的一种数值模型。花少震等[21]采用了三维黏弹性模型来模拟熔体填充过程，该算法借鉴了有限元框架下的黏弹性离散分裂算法。Gao等[22]采用有限元与非连续Galerkin数值算法相结合的方法，建立了一个计算模型来描述非等温聚合物填充过程中错综复杂的气体夹带现象，通过模拟不规则嵌体的空腔填充过程验证了模型的有效性。Deng等[23]通过粒子玻尔兹曼方法建立了一个高效的聚合物注射成型过程仿真两相流模型。但现有的两相流模型大多是不可压缩模型，忽略了密度变化的影响。然而在注塑过程中，随着模具腔体逐渐被填充，熔体和空气互相压缩，其密度都有一定的变化。因此，有必要采用可压缩的两相流模型，并构建相关的计算流程，对注塑过程进行模拟。

笔者基于现有的超算平台，以及有限体积法计算程序，开发用于注塑模拟的三维可压缩两相流数值计算模型，采用压力隐式算子分割算法(PISO)求解压力-速度耦合方程组，结合流体体积(VOF)技术进行界面捕获，使用Cross-WLF黏度模型对聚合物熔体的流变行为进行表征。采用该数值计算模型分别模拟带圆柱嵌件的平板型腔、不均匀厚度的带凹槽平板熔体填充过程，并将模拟结果与试验结果、商业软件模拟结果进行比较。

1 数值计算模型概述

1.1 控制方程

在本文模型中，考虑了密度变化的影响，即空气和聚合物熔体是可压缩的。在填充阶段，流动由质量、动量和能量守恒方程控制，表达式如见(1)、式(2)、式(3)。

式中：ρ，t，u，σ分别代表流体密度、时间、速度矢量和总应力张量；Cp，T，κ，β，p分别代表比热容、温度、热导率、压缩系数和压力。

σ由式(4)给出。

式中：I是单位张量；τ是偏应力张量。

对于可压缩广义牛顿流体的流动，其定义为式(5)的形式。

式中：η为黏度为剪切速率为应变率张量。

熔体相和空气相的流变行为由单流体两相模型定义，对于空气相，假设运动黏度是恒定的，则其为牛顿流体。对于聚合物熔体，采用Cross-WLF黏度模型[24]。该模型被广泛用于研究注塑过程的填充和保压阶段，由式(6)给出。

1.2 自由界面追踪方法

注射成型加工过程中涉及到高温熔融塑料在模具中的填充、冷却和固化过程，在填充阶段的模拟需要求解流体流动方程并在每个时间步长捕获熔体前沿位置。VOF方法[25-28]作为一种常用的多相流模拟方法，可以准确地描述界面变化，特别适用于液体和气体之间的界面问题。为此，使用了VOF方法跟踪熔体流动前沿，其两相流模型示意图如图1所示。

图1 VOF法两相流模型示意图Fig. 1 Schematic diagram of two-phase flow model of VOF method

在此方法中，物质传输方程用于跟踪两相的相对体积分数或相分数α分布，其中标量α用于以速度u传输的液相分数。图1中，α=1表示液态聚合物，α=0表示气态空气，相关方程见式(7)至式(12)。

式中：ur为相对速度矢量；v为运动黏度；cv为相应比热容是加权热导率；下标m和g分别表示熔体相和空气相。

流体的动态黏度μ由密度和运动黏度计算得出，见式(13)。

1.3 数值算法和边界条件

常用的数值离散方法包括有限差分法、有限元法、有限体积法等。有限差分法在离散过程中会忽略方程的守恒性，并且对于求解区域有局限性，处理不了复杂的几何模型。有限元方法通过应用泛函变分原理，利用插值函数将偏微分方程离散化为代数方程，从而获得方程的近似解。然而，当面对复杂的流动问题时，有限元方法难以确保方程的守恒性，并且需要大量内存。相比之下，有限体积法可以将复杂的计算领域划分为小的控制体积进行离散化，并在每个体积中求解物理守恒方程。因此，有限体积法具有较强的守恒性，能够适应复杂的几何形状，并具有较快的计算速度。但采取有限体积法求解三维注塑流动的问题难点在于对压力-速度耦合方程的求解，这是由于速度和压力相互影响，动量方程中的压力是作为动量源项，缺少关于压力的独立方程。为应对这个问题，常用的耦合算法包括SIMPLE (Semi Implicit Method for Pressure Linked Equations)算法[29]和PISO算法[30]。SIMPLE算法通过简化压力和速度修正来进行迭代计算，但其收敛性较慢。而PISO算法在每个时间步长内都增加了迭代步，并考虑了速度修正。通过循环计算直到收敛，具有数值稳定、精度高、耦合效果好等优点。故笔者采用PISO算法解决压力-速度耦合方程。其算法计算流程如图2所示。

图2 PISO算法框图Fig. 2 Block diagram of PISO algorithm

根据式(2)和式(4)，可以得到黏性条件下塑料熔体在型腔内流动的动量方程，见式(14)。

按照一般化公式的离散方法对式(14)进行离散，其中瞬态项的离散采用欧拉格式，可得时间t+Δt关于速度未知量u的线性方程组，见式(15)。

式中：P和N分别为两个相邻的控制单元；aP和aN分别为该两个控制单元的速度系数；H(u)P系数与速度系数有关；f为该两个控制单元的公共面；S为控制单元P的f面上有向面积；lNP为该两个控制单元中心的距离；Ef为控制单元P的f面上有向面积分向量；VP为控制单元的体积；max{(ρuP·S),0}等于(ρuP·S)和0之间的最大值。

对于计算模型的边界条件，浇口处速度与注塑速度相等，温度采用第一类边界条件，压力采用第二类边界条件；模壁处的速度采用无滑移边界条件，温度为对流换热边界条件，压力为第二类边界条件。本模型支持OpenMPI平台，可用于超算平台上的大规模并行计算。

2 数值模拟算例

2.1 带嵌件型腔填充过程模拟

在本算例中使用的材料是韩国锦湖石油化学公司生产的750级丙烯腈-丁二烯-苯乙烯塑料(ABS)[3]，其热物理性能如比热容、导热系数、密度分别为2 000 J/(K·kg)，0.19 W/(m·K)，960.6 kg/m3，其7个cross-WLF黏度模型参数n，τ*，D1，D2，D3，A1，A2分别为0.246 6，83 300 Pa，4.59×1010Pa·s，373.15 K，0，23.522，51.6 K。该算例几何模型如图3所示。型腔几何以O点为原点建立直角坐标系，x方向为几何体长度方向，y方向为几何体宽度方向，z方向为几何体厚度方向，其中厚度、宽度和长度分别为3.5，50，100 mm，塑件内孔直径为20 mm，其中心位于距入口50 mm处。为了分析不同填充时间(t)下宽度方向和厚度方向上物理量的变化，取直线L1 (x=20 mm，z=0 mm)，L2 (x=20 mm，z=1.75 mm)，L3 (x=65 mm，z=1.75 mm)，L4 (x=20 mm，y=25 mm)，L5 (x=20 mm，y=0 mm)。型腔浇口处的熔体填充速率为8×103mm3/s，浇口熔体温度为235 ℃，壁面温度为80 ℃。

图3 带圆柱嵌件的平板型腔几何模型Fig. 3 Geometric model of a flat plate cavity with cylindrical insert

几何体采用六面体为主导的划分方法，选取3种不同的网格系统M1 (厚度方向共4层，网格数为37 225)，M2 (厚度方向为7层，网格数为106 876)，M3 (厚度方向为12层，网格数为302 582)，图4展示了网格M2的划分情况。图5给出了填充结束时刻在x=20 mm和x=50 mm截面上不同网格的速度分布云图，图6给出了在直线L2 (x=20 mm，z=1.75 mm)和直线L3 (x=65 mm，z=1.75 mm)上不同网格的速度分布曲线，通过定量比较不同网格下的速度，可以发现其结果基本吻合，从而验证了网格的收敛性。考虑到计算精度和效率，本算例模拟采用了网格M2，其网格划分如图4所示。

图4 带圆柱嵌件的平板型腔的网格划分Fig. 4 Grid division of a flat plate cavity with cylindrical insert

图5 不同网格下在x=20 mm和x=50 mm截面上的速度分布云图Fig. 5 Velocity distribution cloud diagram at x=20 mm and x=50 mm cross section with different meshes

图6 不同网格情况下直线L2、直线L3的速度分布曲线Fig. 6 Velocity distribution curves of straight lines L2 and L3 in different grid cases

图7给出了填充结束时刻3种网格情况下在x=20 mm和x=50 mm截面上的温度分布云图。图8和图9给出了填充结束时刻网格M2下在直线L1 (x=20 mm，z=0 mm)、直线L2 (x=20 mm，z=1.75 mm)和直线L4 (x=20 mm，y=25 mm)、直线L5 (x=20 mm，y=0 mm)上的温度分布曲线。由图8至图9可知，在高温聚合物熔体接触到低温的模壁和嵌件后，其温度会迅速下降，这是由于热传导的作用，导致模壁和嵌件周围的温度相对较低。

图7 不同网格下在x=20 mm和x=50 mm截面上的温度分布云图Fig. 7 Temperature distribution cloud diagram at x=20 mm and x=50 mm cross section with different meshes

图8 网格M2下直线L1和直线L2的温度分布曲线Fig. 8 Temperature distribution curves of straight lines L1 and L2 under grid M2

图9 网格M2下直线L4和直线L5的温度分布曲线Fig. 9 Temperature distribution curves of straight lines L4 and L5 under grid M2

图10为用本算法模拟的直线L4位置在不同时刻的剪切速率分布曲线。图10中，第一条曲线从t=0.4 s开始到最后一条曲线t=1.5 s结束，每条曲线的时间间隔Δt为0.05 s，箭头方向则代表着时间递增。从图10可以看出，剪切速率在模壁处相对较小，随填充时间增加，剪切速率逐渐接近零，这是因为塑料熔体在模壁处开始固化，时间越长，固化程度越高，故在凝固层处的剪切速率越小；在非模壁处，随着填充时间的增加，型腔剪切速率逐渐增大；此外剪切速率分布在厚度方向上近似呈“M”形，这与文献[31]中分布规律相同。该结果也展示了三维数值模型的优势。

图10 L4位置模拟的厚度方向z上不同时刻剪切速率大小分布Fig. 10 Distribution of shear rate magnitude at different moments in thickness direction z simulated at position L4

熔体流体前沿追踪是注射成型填充模拟的一个关键问题，图11给出了熔体绕圆柱嵌件填充过程的模拟结果与实验结果[32]对比。可以发现，在t=0.6，0.8，1.4 s时刻下熔体流动前沿以及填充量的实验和模拟结果基本一致，满足流动前沿的曲率随着填充时间的增加而减小的规律。即在t=0.8 s时刻下熔体接触圆柱嵌件壁，熔体温度会下降，导致黏度增加，流速变慢，从而曲率减小。验证了笔者计算模型的有效性和准确性，说明模型能够有效追踪型腔内熔体流动前沿的位置和形状。

图11 熔体前沿位置(上：实验结果；下：数值模拟结果)Fig. 11 Melt front position (top: experimental results；bottom: numerical results)

在注射成型中由于质量输送、温度场和压力场的变化等因素的耦合作用，塑料熔体填充过程是一个典型的瞬态多场耦合问题，这类问题在模拟中易出现数值振荡、求解溢出等问题。图12显示了在不同时刻z=1.75 mm和y=25 mm的中截面上压力分布云图。从图12可以清晰地观察到，在填充过程中，带圆柱嵌件型腔中的压力在任何时刻都具有光滑、无振荡和关于中轴线对称的特点，说明本文模型可以得出较为稳定的计算结果。浇口处的压力在填充时间为2.28 s时达到最大值，且压力会沿着熔体流动方向逐渐降低，型腔末端处的压力最小。此外，在薄壁型腔中，压力沿着型腔厚度方向分布比较均匀，其压力分布规律与文献[33]结果基本吻合。图13a为P1，P2，P3，P4和P5检测点位置分布，图13b为5个检测点在填充过程中具体压力数值随时间的变化。从图13b看出，各个检测点在不同时刻下压力随填充时间增加而增大，P3和P4处压力曲线基本保持一致，更加验证在填充过程中压力是光滑、中轴线对称的特点。

图12 填充过程不同时刻压力分布云图Fig. 12 Pressure distribution cloud diagram at different moments of filling process

图13 检测点在填充过程中压力数值的变化Fig. 13 Variation of pressure values at detection point during filling process

2.2 不均匀厚度型腔填充过程模拟

在本算例中，使用的材料是李长荣化学工业股份有限公司生产的牌号为6733的聚丙烯(PP)，这是一种广泛使用的注塑材料，其7个cross-WLF黏度模型参数n，τ*，D1，D2，D3，A1，A2分别为0.219，56 351 Pa，1.3×1014Pa·s，263.15 K，0，28.44，51.6 K。该算例几何模型如图14所示。型腔几何形状的中间厚度为4 mm，两边厚度为8 mm，长度为120 mm，入口边界也显示在图中。入口熔体温度为240 ℃，壁面温度为50 ℃，填充时间为0.5 s。网格划分如图15所示，采用六面体网格，其网格总量为50 353个。

图14 不均匀厚度凹槽平板型腔的几何模型Fig. 14 Geometric model of a flat plate cavity with uneven thickness groove

图15 不均匀厚度凹槽平板型腔的网格划分Fig. 15 Grid division of a flat plate cavity with uneven thickness groove

图16为HsCAE3D，Moldflow MPI_3D[34]和本文计算模型模拟在t=0.125 s，t=0.25 s和t=0.375 s上的流动前沿。根据图16的比较结果可以看出，3种模型得到的流动前沿位置和形状基本保持一致。当塑料熔体从浇口注入型腔时，会首先在薄壁凹槽内流动，随后流进厚壁区域。这是由于厚壁区域内的阻力较小，更多的塑料熔体可以像“快车道”一样向前流动，流动速度也会迅速增加，从而使得厚壁内的流动前沿位置超过薄壁区域的流动前沿位置[35]，而厚壁区域中心处的流体速度也会远远高于平均值，从而使得塑料熔体的流动更加顺畅，进而达到最佳的流动效果。由于壁面黏性高于中心层，导致中心层流动也偏快，其会产生形状像喷泉的流锋，这种现象被称为“喷泉流”。这种流动现象无法被2.5维模型捕捉到，但可以通过三维模型预测。

图16 Moldflow MPI_3D，HsCAE3D和本文计算模型模拟的流动前沿对比Fig. 16 Comparison of flow fronts simulated by Moldflow MPI_3D，HsCAE3D and calculation model in this paper

将笔者开发的模型模拟得到的浇口压力与Moldflow MPI_ Fusion，Moldflow MPI_3D以及HsCAE3D的模拟结果对比，得到图17中的曲线。从图17可以看出，表面流MPI_ Fusion与其他3种三维模型的浇口压力值随着填充时间增加，差距不断变大，开发的计算模型模拟得到的浇口压力则与Moldflow MPI_3D模拟结果基本接近。图18展示了不同填充流量下所需注射压力的变化，列出了开发的计算模型与Moldflow MPI_3D模型以及HsCAE3D的计算结果。图18结果表明，开发的模型模拟结果与已有的三维模型计算结果基本吻合。随着填充流量的增大，注射压力数值与Moldflow MPI_3D结果更接近。

图17 开发的计算模型与Moldflow MPI_ Fusion，Moldflow MPI_3D及HsCAE3D模拟的浇口压力对比Fig. 17 Comparison between the simulation gate pressures of developed calculation model，Moldflow MPI_ Fusion，Moldflow MPI_3D and HsCAE3D

图18 开发的计算模型与Moldflow MPI_3D及HsCAE3D预测的不同流量所需的注射压力Fig. 18 Required injection pressure for different flow rates predicted by developed calculation model，Moldflow MPI_3D and HsCAE3D