基于Giesekus黏弹性本构的注射成型数值算法研究

花少震，曹伟

（1.河南工学院机械工程学院，河南 新乡 453003；2.郑州大学橡塑模具国家工程研究中心 ，郑州 450001）

0 前言

注射成型是塑料成型加工中广泛采用的加工工艺之一，其可高效率经济地实现复杂塑料制件加工成型。正是由于这个原因，注射成型的仿真比其他塑料加工工艺更加成功［1］。然而，注射成型的仿真仍然存在一些挑战。黏弹性计算是注射成型模拟的一个重大挑战之一。黏弹性计算的第一个难点是描述注射成型过程中熔体黏性和弹性效应的本构模型。黏弹性本构模型所预测的应力对于计算分子取向、双折射、收缩、翘曲至关重要［2］。丁军等［3］对国内外橡胶材料黏弹本构关系的研究现状进行了系统的总结和分析，并对橡胶材料黏弹本构关系未来的发展方向提出了建议。陈腾飞等［4］采用流变实验研究了高抗冲聚苯乙烯（PS-HI）黏弹性特征，发现黏弹性Phan-Thien Tanner（PTT）模型的精度明显高于黏性模型。对于注射成型的黏弹性流动模拟，Leonov模型是较早使用的黏弹性本构模型之一［5-6］；Chang等［7］采用标准K-BKZ黏弹性本构模型，模拟了薄壁零件的注射成型；张世勋等［8］、Cao等［9］、Olley［10］采用PTT黏弹性本构模型模拟了熔体充填过程中的流动诱导应力。eXtended Pom-Pom（XPP）［11-12］、Giesekus模型［13-14］、Rolie-Poly模型［15］也被采用来描述注射成型过程中熔体的黏弹性特性。尽管有很多黏弹性模型可以被用来描述注射成型模拟中的熔体流变行为，但到目前仍然没有一个最优黏弹性模型。

注射成型的黏弹性模拟的另一个难点是黏弹性本构模型的求解算法。其难于如何开发准确、稳健、稳定的算法，尤其是针对众所周知的高魏森伯格数（Wi）问题，即当Wi超过一定限度时，数值解会出现振荡、不收敛等问题。在黏弹性数值算法的研究中，有限元法得到了广泛的应用，在早期的研究中常采用混合有限元法［16］。一般来说，使用这种混合有限元法解决问题时，存在两个问题。首先，随着Wi的增加，对流项的影响逐渐增大，Galerkin离散化将导致数值解的不稳定。其次，每个未知物理量的离散空间都需满足Ladyzhenskaya Labuska Brezzi（LBB）条件。针对这些问题，Marchal［17］将streamline upwind /Petrov Galerkin（SUPG）［18］方法应用于黏弹性流体流动，并提出了一个新的流线上风（SU）方案。为了保证数值计算的稳定性，许多研究都集中于在动量方程中加入更多的扩散处理，使动量方程更具椭圆性。Beris［19-20］提出了弹性黏性应力分裂（EVSS）方案，该方案将应力分为黏性应力和弹性应力，并将剪切率张量作为一个新变量。这种方法有很大的局限性，只适合于UCM。更重要的是，增加的项包含剪切率张量的导数，这极大增加了计算量和求解的难度。为了克服这些缺点，Rajagopalan［21］采用最小平方法来近似剪切率的流动导数，使EVSS方案能够适应更多的黏弹性构成方程。后来研究人员在EVSS的基础上发展了一些新的方法，如DEVSS［22］、DEVSS-G［23］等。然而，当魏森堡数增加到一定程度时，仍然会出现计算不稳定。此外，一些无网格法也用于黏弹性流动求解，例如SPH方法［24］、格子玻尔兹曼方法［25］、MPS方法［26］等。

本文的目的是发展一种三维黏弹性求解算法，模拟注射成型过程中的熔体的黏弹性填充过程，以实现注射成型过程中熔体填充过程中的浇口压力、熔体前沿、熔接线、流动诱导应力的预测。

1 数学模型

假设熔体充填过程是不可压缩非等温的，考虑熔体流动过程中的惯性和重力影响，则熔体充填过程中连续性方程、动量方程和能量方程：

其中u、p、T分别为速度、压力、温度。ρ和η、t、g分别代表熔体密度、速度、时间和重力加速度。Cp、k为等压热容和热传导率。τ为偏应力张力，可以分为两部分：

其中，τs=ηs(∇u+∇uT)为黏性所贡献偏应力；τp为弹性所贡献偏应力，由黏弹性本构方程求解所得。对于n个模态黏弹性本构方程，τp满足本文采用Giesekus黏弹性本构模型描述熔体充填过程中的黏弹性。Giesekus黏弹模型不仅能表征聚合物的剪切变稀和拉伸变硬流变特征，其在稳态流场中的第一和第二法向应力差同剪切速率的平方具有高度的非线性关系，即说明Giesekus黏弹本构模型能表征非线性弹性流变行为。τi满足：

式（5）中，i代表n个模态的编号，λi为松弛时间；ηi为黏度，且满足等于零剪切黏度；αi为分子流动常数。为上随体导数：

为了追踪熔体在型腔内的流动前沿，采用VOF方程描述熔体在型腔内的充满程度：

其中，ϕ为充填因子，0≤ϕ≤1。当ϕ=1代表熔体充满控制体，ϕ=0代表控制体内尚未有熔体流入，0<ϕ<1时代表控制体内充有空气和熔体，即为熔体充填流动所在前沿。

为了准确求解熔体在型腔内的黏弹性充填，假设在型腔模壁上为无滑移无渗透边界，法向偏应力张量梯度为零，入口速度为恒定值，入口偏应力张量为零。假设充填过程中，型腔排气良好，即出口处速度和偏应力张力的法向梯度均为零为零，压力恒为大气压。对于温度，入口温度即为熔体的加工温度，模壁和流动前沿温度采用第三类边界条件：

其中，Tb是模壁或空气温度，ha是熔体与模壁或熔体与空气之间的传热系数。

2 数值算法

2.1 FVM-DEVSS方法

高魏森堡数（high-Weissenberg number）问题已经困扰了计算流变学30多年了。它极大地限制了黏弹性材料加工中模拟的应用。为了克服这个问题，这些年来已经开发了几种数值方法，如黏性形成、EVSS、DEVSS、AVSS、EEME、对数构象张量法等。本文将Guénette和Fortin［22］在有限元方法下的DEVSS方法发展到有限体积框架。DEVSS方法在动量方程的离散形式中引入了一个稳定的椭圆算子以保证稳定性。而在有限体积法离散框架下，本文首先在动量方程扩散项中引入椭圆算子，然后采用有限体积法离散，其中扩散项中椭圆算子做隐式处理，源项中椭圆算子作为显式处理。即：

其中，等号右边ϕ(ut)=-∇p+∇⋅τp+ρg-∇⋅(ηs∇ut)作为源项处理，t和t+Δt代表充填时间，Δt为时间步长。采用有限体积法离散式（9），在t+Δt时间点有：

其中，nf和Af分辨代表控制单元f面上的单位外法向量和面积。Vp为控制体p的体积。

2.2 Giesekus黏弹本构方程离散

省略下标i将式（5）整理为瞬态项、对流项和源项，可得：

式（11）中，S(τ)为源项，S(τ)=(∇u)T⋅τ+τ⋅

采用有限体积法离散式（11），在时间t+Δt有：

式（12）中，τf为控制单元f面上的偏应力张量。在每一时间步长，上式中等号右边作为已知量，其计算采用上一时间步长的求解结果。

3 结果与讨论

3.1 算法验证

为了验证本文算法的有效性，采用所发展黏弹性算法计算半4∶1平面收缩流基准问题。首先将Giesekus黏弹性本构方程无量纲化：

式（13）中We即为魏森堡数，其数值越高，说明黏弹性算法稳定性越好；β=然后依照章节2.2中算法离散上式，计算等温条件下充分发展的半4∶1平面收缩流。图1（a）为α=0.2、β=0.2、We=8.25时半4∶1平面收缩流的流线分布。We=8.25时，根据流线分布发现本文算法能计算出半4∶1平面收缩流的涡流。图1（b）为a、b两条线上第一法向应力差计算结果和文献［27］中结果对比，其中，l1/H=0.0925，（l1+l1）/H=0.425，a、b两条线上第一法向应力差计算结果和文献［27］中误差小于5%，离涡越近误差越大。因此，可以认为本文算法能有效计算黏弹性流动。

图1 4∶1半平面收缩流数值结果Fig.1 Numerical results of half 4∶1 planar contraction flow

3.2 流动前沿及浇口压力模拟

注射成型过程中熔体流体前沿追踪是填充模拟的一个关键问题。本节采用实验和数值手段相结合的方式验证本文所发展黏弹性算法能有效追踪型腔内熔体流动前沿的位置和形状。图2为模具流道和型腔的几何形状及尺寸。实验设备和材料分别为Demag注塑机（System 80/420-430）、聚苯乙烯（PS，PG-33，奇美公司）。用旋转流变仪测得材料参数为：ηpηs=8，η0=258 Pa⋅s，λ=0.18 s，α=0.2。熔体和模具温度分别为240、35 ℃，充填时间为1 s。然后用上文发展算法模拟相同条件下熔体的充填过程。

图2 流道和型腔几何结构及尺寸Fig.2 Geometry and dimension of the runner and cavity

图3为t=0.22、0.28、0.42 s时充填的实验和数值结果对比图。可以发现，3个时刻下熔体在型腔内流动前沿以及充填量的实验和数值结果均非常吻合。图3（a）中当t=0.22 s时，流动前沿的曲率较大，随着充填时间的增加，流动前沿曲率逐渐减小。对比可以发现，数值模拟结果［图3（b）］中流动前沿同样满足随时间增大曲率逐渐减小的规律。即本文发展的黏弹性充填算法可以准确预测流动前沿的位置和形状。

图3 不同时间点型腔内熔体充填形态Fig.3 Experimental and umerical results at different filling time

为了进一步分析本文算法的有效性，相同条件下采用Cross-WLF黏度模型下三维（3D）充填算法［28］以及采用模流分析软件Moldflow的中性面（2.5D）算法模拟熔体的充填流动，并且提取浇口压力以进行对比，如图4所示。从图中可以发现，3种情况下浇口压力误差小于2%。在仅考虑熔体黏性条件下，3D算法计算的浇口压力始终大于2.5维，而且随着充填时间的增加，二者的差值逐渐增大，在充填末端达到0.45 MPa。而当考虑熔体的黏弹性质后，在充填前期（t<0.44 s）浇口压力大于仅考虑黏性条件下的数值结果，最大达到0.60 MPa，此阶段熔体的弹性和黏性影响浇口的压力；当t>0.44 s，由于熔体的松弛效应，此时浇口压力小于仅考虑黏性条件下的数值结果，而且差值随着充填增大，此阶段熔体的黏性占主导地位，如图中所示在充填末期三维黏性和黏弹性算法下所计算浇口压力差别较小。所以本文所发展的黏弹性算法不仅能准确模拟出熔体流动前沿位置的形状，而且能表达出熔体充填过程中黏性和弹性对浇口充填压力的影响规律。

图4 3种情况下浇口压力对比Fig.4 Comparison of gating pressure in three cases

3.3 熔接线与流动诱导应力模拟

为了模拟熔接线位置以及流动诱导应力，采用文献［29］中相同的型腔结构和工艺条件。型腔结构及尺寸如图5所示，“▲”为浇口所在位置。由于矩形型腔带有孔洞，所以熔体从浇口射入后在充填过程中会围绕矩形障碍物充填，即在型腔末端会有两股熔体相遇，两股熔体相遇后即会产生熔接线。图6为采用本文黏弹性算法所模拟的熔接线形成过程。可以发现随着两股熔体绕过矩形障碍相遇后，在两股熔体相融处即会产生熔接线。根据型腔结构的对称性，t=0.9、1、1.05 s这3个时刻充填图清晰地揭示了两股熔体相遇时所生成的直线型熔接线。图7为两股熔体相遇时的熔接线生产试验和数值模拟结果，可以发现本文发展的黏弹性算法模拟的熔接线和实验结果［29］非常吻合。即说明可以用本文算法预测注射成型过程中两股熔体相遇时熔接线的生成。

图5 带孔洞的矩形型腔结构与尺寸Fig.5 Geometry and dimension of rectangular cavity with hole

图6 充填过程的数值模拟Fig.6 Numerical simulation of filling

图7 熔接线的形成过程：（a）数值结果、（b）实验结果［29］Fig.7 Forming of weld line ：（a） numerical results；（b）experimental results

流动诱导应力主要影响塑料制件的光学性能。为了分析本文所发展算法对于流动诱导应力的适用性，对比了双折射实验结果［29］和数值模拟结果，如图8所示，其中图8（a）为双折射实验结果［29］，（b）为数值模拟结果。在熔接线处，由于熔体融合程度低，除了影响力学性能，也影响透明件的光学性能，所以透明件双折射实验中熔接线附近的光弹条纹密集，图8中实验和数值结果均说明了这一规律。此外，在模壁处，由于法向偏应力张量梯度为零，所以型腔宽度和长度方向的主应力差值梯度较大，即在模壁处光弹条纹较为密集，从图8中实验和数值结果同样论证了这个结论。而且对比图8中实验和数值结果，可以发现二者光弹条纹的分布规律非常吻合，说明了本文算法能有效模拟注射成型过程中熔体充填的流动诱导应力，这对于提高塑料制件的光学性能、预测收缩、翘曲等具有重要意义。

图8 流动诱导应力数值模拟与实验结果对比：（a）双折射实验结果［29］、（b）数值模拟结果Fig.8 Comparison of numerical simulation and experimental results of flow-induced stresses ：（a） Birefringence experimental results［29］；（b） Numerical simulation results

4 结论

（1）本文所发展算法可以计算魏森堡数为8.25时充分发展的半4∶1平面收缩流；

（2）本文所发展算法可以准确预测注射成型过程中熔体的流动前沿的形状和位置、浇口压力大小、熔接线位置和形状，并且可以用之计算充填过程中熔体的流动诱导应力。