CN103544341B - 基于有限体积法的气辅注塑成型气体穿透预测方法 - Google Patents

Abstract

本发明公开一种基于有限体积法的气辅注塑成型气体穿透预测方法，包括以下步骤：(1)由型腔实体网格构造三维节点控制体作为有限体积；(2)设置成型工艺参数；(3)设置材料参数；(4)计算塑料熔体流动速度、压力和温度；(5)更新塑料熔体流动前沿，确定时间步长；(6)塑料熔体填充量是否达到指定体积；(7)计算塑料熔体和气体的流动速度、压力和温度；(8)更新塑料熔体流动前沿和气体穿透前沿，确定时间步长；(9)型腔是否充满。本发明基于三维有限体积法，通过改进传统的有限体积法，计算塑料熔体与气体的耦合流动，考虑了粘度、速度、压力、温度之间的强非线性，能准确预测气体辅注塑成型气体穿透前沿及各种物理量场的分布。

Description

基于有限体积法的气辅注塑成型气体穿透预测方法

技术领域

本发明属于计算机数值模拟方法，特别涉及一种基于有限体积法的气辅注塑成型气体穿透预测方法，用于气体辅注塑成型。

背景技术

气体辅注塑成型(gas-assisted injection molding，GAIM)是在传统的注塑成型的基础上，通过气针将气体注入塑料熔体内部，气体压力推动塑料熔体继续流动充满模具型腔，气体在穿透过程中形成塑料熔体内部空洞，由此得到中空塑料制品的先进注塑成型工艺。GAIM工艺具有很多优点，如节省原料、减小合模力、缩短冷却时间、防止制品缩痕、减少制品内应力、减少或消除制品翘曲、提高制品表面性能、提高生产效率、降低生产成本等，增加产品的市场竞争力。GAIM技术从20世纪70年代中期开始发展起来，Ferromatik、Battenfeld、KraussMaffei、Cinpress、Maximator、GAIM技术等公同对GAIM技术进行了不断的开发和完善。目前，GAIM技术在发达国家成功地应用于各种产品如汽车、电脑、手机、家电、家具、包装、农用塑料制品、运动器械等领域。

目前具有代表性的GAIM技术主要有外部气辅注塑成型、振动气辅注塑成型、冷却气辅注塑成型、多腔控制气辅注塑成型、气辅共注射成型技术等。根据预注塑量的不同，气辅注塑成型分为欠料注射和满料注射。气辅注塑成型工艺包括塑料熔体注射、气体穿透、冷却、顶出等步骤，其中气体穿透是至关重要的一步，如何控制和优化气体穿透前沿是GAIM的关键技术。熔体温度、注射速度、预注塑量、延迟时间、气体压力、气针位置、气道形状和长度、产品厚度等诸多工艺和产品因素对气体穿透前沿都有重要影响，可能导致产品缺陷，如气体吹破、迟滞线、放射斑、翘曲、缩痕、气指、短射等，可以通过调整工艺和产品设计来避免缺陷的产生。然而这些问题的解决如果依赖传统的经验法、试错法，即不断地反复试模后再评定产品设计和工艺参数，会使得GAIM工艺的成本高、开发周期长，且产品的废品率高、生产不稳定，因此，会浪费大量的人力物力。而基于数值模拟的方法能预测GAIM的气体穿透过程，为GAIM产品和工艺设计提供了一种经济、实用和科学的手段。

有学者开发了数值模拟软件研究GAIM的气体穿透前沿，但采用的均为有限元法，如南昌大学的周国发、柳和生等开发的三维有限元数值模拟软件，华中科技大学周华民等开发的基于中面和表面模型的有限元数值模拟软件，郑州大学的申长雨等开发的二维有限元数值模拟软件，此外，也有研究者采用商业软件Moldflow模拟GAIM的气体穿透前沿。基于中面、表面模型和二维的数值模拟方法引入简化模型描述气体的穿透前沿，无法模拟厚度方向上，特别是产品较厚处、厚度变化处、结构交叉处复杂的气体穿透前沿。有限体积法兼有有限元和有限差分的优点，能同时满足局部守恒和整体守恒，完全兼容非结构的有限元网格，因此三维有限体积法较有限元法能更好地模拟GAIM的气体穿透前沿。

随着全球市场竞争日趋激烈，要求尽可能科学、准确地设计GAIM产品和设置工艺参数，提高产品质量与生产效率，因此，需要一种全面科学地预测气体穿透前沿的方法。

发明内容

本发明目的在于提供一种基于有限体积法的气辅注塑成型气体穿透预测方法，同时满足局部守恒和整体守恒，完全兼容非结构的有限元网格，用于模拟三维GAIM中气体穿透前沿。本发明采用节点控制体作为限体积，考虑到塑料熔体和气体耦合流动中相邻节点控制体的粘度较大差异，根据相邻节点控制体公共界面上流动剪切应力相等的原则，改进了传统的有限体积法，如节点控制体中心的速度梯度、节点控制体公共界面的速度和速度梯度的计算方法，提高了有限体积法模拟空间粘度差异较大的GAIM过程的数值精度和数值稳定性，避免了数值震荡。

本发明目的是通过以下技术方案来实现的：

本发明涉及一种基于有限体积法的气辅注塑成型气体穿透预测方法，包括以下步骤：

A、由型腔实体网格构造三维节点控制体作为有限体积；

B、设置成型工艺参数；

C、设置材料参数；

D、计算塑料熔体流动速度、压力和温度；

E、更新塑料熔体流动前沿，确定时间步长；

F、塑料熔体填充量是否达到指定体积，如果达到，进入步骤G，否则进入步骤D；

G、计算塑料熔体和气体的流动速度、压力和温度；

H、更新塑料熔体流动前沿和气体穿透前沿，确定时间步长；

I、模具型腔是否充满，如果模具型腔已经填充满，则计算结束，否则进入步骤G。

优选地，所述步骤A的由型腔实体网格构造三维节点控制体作为有限体积步骤进一步包括：导入型腔实体网格，连接网格单元的体心、面心、边中点构造三维节点控制体作为有限体积。计算中的变量，如速度、压力、温度等均位于节点控制体的中心，即节点上。

优选地，步骤B中，所述成型工艺参数包括塑料熔体的注射压力、注射时间或注射速度，气体的注射压力、注射时间或注射速度，以及模具温度、熔体温度。

优选地，步骤C中，所述材料参数包括气体粘度，塑料熔体的密度、比热容、热传导系数，以及塑料熔体的7参数WLF-Cross粘度模型的材料常数。

进一步优选地，所述材料常数包括 A₁、A₂、D₁、D₂、D₃。

优选地，所述步骤D的计算塑料熔体流动速度、压力和温度步骤中，忽略塑料熔体的惯性力、体积力，塑料熔体为不可压缩流体，其速度、压力、温度等变量位于节点控制体中心(即节点处)，塑料熔体流动的速度、压力和温度的计算步骤进一步包括：

(1)根据7参数WLF-Cross粘度模型计算节点控制体中心的粘度，塑料熔体粘度为温度、压力、剪切应变速率的函数：

$\mathrm{\η}(T,\stackrel{\·}{\mathrm{\γ}})=\frac{{\mathrm{\η}}_0(T,p)}{1+{\left(\frac{{\mathrm{\η}}_0\stackrel{\·}{\mathrm{\ϵ}}}{\stackrel{\·}{\mathrm{\τ}}}\right)}^{1-\tilde{n}}},{\mathrm{\η}}_0=D_1\exp \left(\frac{-A_1(T-(D_2+D_{3}p))}{A_2+T-D_2}\right)$

式中，η为粘度，η₀为零剪切粘度，T为温度，p为压力，为非牛顿指数，为剪切应变速率，为材料剪切常数，A₁、A₂、D₁、D₂、D₃为材料参数，其中T、p、为节点控制体中心的变量值，由速度、压力和温度计算得到；

(2)塑料熔体速度、压力、温度的控制方程分别采用如下动量守恒方程、质量守恒方程和能量守恒方程：

$-p_{,i}+2\mathrm{\η}{\stackrel{\·}{\mathrm{\ϵ}}}_{\mathrm{ij},j}=0$

u_i，j＝0

ρC_p(T_，t+u_iT_，i)＝Φ+λT_，ii

式中，为剪切应变速率，下标“’”表示偏导，i，j＝1，2，3为空间局部笛卡儿坐标系分量，u_i为速度分量，p为压力，η为粘度，T为温度，t为时间，ρ、C_P分别为塑料熔体的密度、比热容，λ为塑料熔体的热传导系数，为剪切耗散功，是塑料熔体充模过程的热源项；

(3)采用有限体积法对动量方程进行离散时，考虑到塑料熔体流场中相邻节点控制体的粘度较大差异，根据相邻节点控制体公共界面上流动剪切应力相等的原则，改进传统的有限体积法，得到如下计算公式：

a、某节点控制体中心P₀的速度梯度计算式为：

${\▿u}_{P_0}=D^{-1}\·c$

式中， $D=\underset{k=1}{\overset{n_f}{\mathrm{\Σ}}}{(r_{P_k}-r_{P_0})}^T(r_{P_k}-r_{P_0}),c=\underset{k=1}{\overset{n_f}{\mathrm{\Σ}}}\frac{{(r_{P_k}-r_{P_0})}^T(u_{P_k}-u_{P_0})}{1-w_k+q_k\×w_k},$ u为节点速度矢量，下标“k”表示相邻节点控制体，P_k为相邻节点控制体中心，n_f为相邻节点控制体数，r为控制体中心的位置矢量，上标“T”表示转置向量，为加权系数，r_k为与相邻节点控制体的公共界面中心的位置矢量，为两个节点控制体中心的粘度比；

b、与相邻节点控制体的公共界面中心的速度梯度计算式为：

${\▿u}^{*}=\stackrel{\‾}{\▿u}+[\frac{(u_{P_k}-u_{P_0})}{\left|d_k\right|(1-w_k+q_k\×w_k)}-\frac{\stackrel{\‾}{\▿u}\·d_k}{\left|d_k\right|}]\frac{\left|d_k\right|s_k}{d_k\·s_k}$

式中，u^*为公共界面中心速度矢量， $\stackrel{\‾}{\▿u}=w_k{\▿u}_{P_0}+\frac{(1-w_k){\▿u}_{P_k}}{q_k},d_k=r_k-r_{P_0},$ s_k为相邻节点控制体公共界面的面积外法矢量(由P₀指向P_k)；

c、与相邻节点控制体的公共界面中心的速度计算式为：

$u^{*}=\stackrel{\‾}{u}+\frac{\stackrel{\‾}{V}}{a}[-\frac{(p_{P_k}-p_{P_0})}{\left|d_k\right|}+\frac{\stackrel{\‾}{{\▿p}_k}\·d_k}{\left|d_k\right|}]\frac{\left|d_k\right|s_k}{d_k\·s_k}$

式中， $\stackrel{\‾}{u}=\frac{q_k{w}_k{u}_{P_0}+w_k{u}_{P_k}}{1-w_k+q_k\×w_k},\frac{\stackrel{\‾}{V}}{a}=\frac{q_k{w}_k\frac{V_{P_0}}{a_{P_0}}+w_k\frac{V_{P_k}}{a_{P_k}}}{1-w_k+q_k\×w_k},$ 其中V为节点控制体体积，a为节点速度分量对应方程组矩阵的主对角线上的值；

上述速度梯度和节点控制体公共界面速度、速度梯度的离散方法很好地解决了不同节点控制体粘度的较大差异导致的数值计算困难和数值振荡，保证了节点控制体公共面两侧剪切应力的一致性，提高了有限体积法求解塑料熔体流动速度、压力问题的精度和稳定性；

d、采用有限体积法对能量守恒方程进行离散时，相邻节点控制体公共界面中心温度采用“延迟”算法提高计算精度，保证能量守恒方程计算的精度和稳定性，界面中心的温度计算式为：

T^*＝T^SU+γ(T^SO-T^SU)

式中，T^SU为根据一阶精度的迎风格式得到的温度，T^SO为采用中心差分的得到的温度，具有二阶精度，即γ为0.8-95之间的常数；

e、采用节点型有限体积法耦合求解塑料熔体的动量守恒方程、质量守恒方程和能量守恒方程得到速度、压力和温度。速度、压力和温度之间通过粘度互相影响，采用非线性耦合迭代的方法求解速度、压力和温度直到粘度、速度、压力和温度迭代变化在指定的很小范围内。

优选地，所述步骤E具体为：根据计算得到的速度场，采用流动网络分析法，选择当前流动前沿中各节点控制体所需的最短填充时间作为时间步长，更新流动前沿中其它节点控制体的填充率。

优选地，在所述步骤G的计算塑料熔体和气体的流动速度、压力和温度步骤中，气体穿透推动塑料熔体继续填充模具型腔，塑料熔体和气体之间存在动量交换和质量交换，塑料熔体和气体的流动速度、压力和温度的计算步骤进一步包括：

a、气辅注塑成型过程中，气体粘度小、密度小，不影响塑料熔体的温度分布，且各部分压力变化小，以近似均匀的压力推动塑料熔体继续流动，因此气体穿透只考虑其动量守恒和质量守恒，忽略惯性力、体积力及可压缩性，该简化在不损失计算精度的前提下，大大减少计算量，提高计算速度和稳定性，气体穿透过程满足动量守恒、质量守恒，对应的控制方程分别为：

$0=-p_{,i}+2\mathrm{\η}{\stackrel{\·}{\mathrm{\ϵ}}}_{\mathrm{ij},j}$

u_i，j＝0

式中，为剪切应变速率，下标“’”表示偏导，i，j＝1，2，3为空间局部笛卡儿坐标系分量，u_i为速度分量，p为压力，η为气体粘度；

b、采用节点型有限体积法求解耦合的塑料熔体和气体的动量守恒方程、质量守恒方程和能量守恒方程，忽略气体温度的变化，将气体作为具有一定粘度的假想“塑料熔体”，气体与塑料熔体粘度差异较大，所有节点控制体中心速度梯度、相邻节点控制体公共界面中心速度和速度梯度的计算方式与上述塑料熔体的速度、速度梯度计算方式相同；

c、采用节点型有限体积法求解耦合的塑料熔体和气体的动量守恒方程、质量守恒方程和能量守恒方程得到速度、压力和温度。速度、压力和温度之间通过粘度、界面互相影响，采用非线性耦合迭代的方法求解速度、压力和温度直到粘度、速度、压力和温度迭代变化在指定的很小范围内。

优选地，所述步骤H具体为：根据计算得到的速度场，采用流动网络分析法，选择当前塑料熔体流动前沿和气体穿透前沿中各节点控制体所需的最短填充时间作为时间步长，更新流动前沿中节点控制体的塑料熔体填充率和气体填充率。

与现有技术相比，本发明有以下优点：

1、与有限元法相比，有限体积法更适用于气辅注塑成型过程中的流体数值模拟，与传统有限体积法相比，本发明通过实体网格构造三维节点控制体作为有限体积，速度、压力、温度等变量位于节点处，考虑到塑料熔体和气体耦合流动中相邻节点控制体的粘度较大差异，根据相邻节点控制体公共界面上流动剪切应力相等的原则，改进传统的有限体积法，如节点控制体中心的速度梯度、节点控制体公共界面的速度和速度梯度的计算方法，由此提高有限体积法模拟空间粘度差异较大的GAIM过程的数值精度和数值稳定性，避免数值震荡。能量守恒方程的计算中，相邻节点控制体公共界面中心温度采用“延迟”算法提高计算精度，保证了能量守恒方程计算的精度和稳定性。以上措施能更准确地预测塑料熔体流动和气体穿透前沿及各种变量的分布；

2、本发明引入合理的假设与简化，将气体作为假想“塑料熔体”，塑料熔体和气体只是具有不同粘度的流体，从而可以将塑料熔体和气体作为统一的流体进行处理，且忽略流体的惯性力、体积力、可压缩性，由此既可以准确地预测GAIM中气体穿透前沿，还大大减少数值模拟的计算量，可以为合理地设计GAIM产品和设置工艺参数提供科学依据，避免产品缺陷的产生，从而提高产品质量和生产效率。

附图说明

通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：

图1为某产品的实体网格模型图；

图2为本发明的流程框图；

图3为本发明构造的节点三维子控制体；

图4为本发明的有限体积示意图；

图5为塑料熔体填充型腔80％时的流动前沿；

图6为气体穿透推动熔体充满型腔时的气体穿透前沿。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干调整和改进。这些都属于本发明的保护范围。

请参阅图1，为某制品的几何模型实体网格图，图中左下角端面中心为塑料熔体和气体的入口。预测气体穿透前沿的步骤如图2所示，具体步骤如下。

(1)由型腔实体网格构造三维节点控制体作为有限体积。请参阅图3，以四面体为例，连接四面体体心、面心、边中点形成的表面将四面体划分为4个节点的4个子控制体，节点在所有四面体中的子控制体积合并为三维节点控制体。

(2)设置成型工艺参数步骤。成型工艺参数包括注射压力或注射时间或注射速度(3者任选1)、模具温度、熔体温度等；设置熔体温度为240℃，模具温度为40℃。

当设置的工艺参数为注射压力时，将注射压力作为压力方程计算的边界条件；

当设置的工艺参数为注射速率时，通过不断调整注射压力，使注入的流量等于或者接近于设置的注射速率来满足注射速率要求；

当设置的工艺参数为注射时间时，按照制品的体积将其转化为注射速率，再按照注射速率的方式求解压力方程；

在本实例中，设置熔体注射时间为1s充满型腔，由此根据型腔体积计算出注射速率为0.0002118m³／s，注射0.8s后停止熔体注射，气体注射压力为6MPa。

(3)设置材料参数步骤。在本实例中，采用PP6533塑料材料和氮气气体，材料参数如表1-2所示，此外，氮气粘度为1.79×10^-5Pa·s。

表1、塑料熔体材料参数

表2、塑料熔体粘度参数

塑料熔体的WLF-Cross粘度模型中，T为温度，P为压力，为剪切应变速率；

(4)计算塑料熔体流动速度、压力和温度步骤。根据步骤(3)中的粘度模型计算不同位置的粘度，塑料熔体粘度为温度、压力和应变剪切速率的函数：

$\mathrm{\η}(T,\stackrel{\·}{\mathrm{\γ}})=\frac{{\mathrm{\η}}_0(T,p)}{1+{\left(\frac{{\mathrm{\η}}_0\stackrel{\·}{\mathrm{\ϵ}}}{\stackrel{\·}{\mathrm{\τ}}}\right)}^{1-\tilde{n}}},{\mathrm{\η}}_0=D_1\exp \left(\frac{-A_1(T-(D_2+D_{3}p))}{A_2+T-D_2}\right)$

式中，η为粘度，η₀为零剪切粘度，T为温度，p为压力，为非牛顿指数，为剪切应变速率，为材料剪切常数，A₁、A₂、D₁、D₂、D₃为材料参数，其中T、p、为节点控制体中心的变量值。

塑料熔体速度、压力、温度的控制方程分别采用如下动量守恒方程、质量守恒方程和能量守恒方程：

$-p_{,i}+2\mathrm{\η}{\stackrel{\·}{\mathrm{\ϵ}}}_{\mathrm{ij},j}=0$

u_i，j＝0

ρC_p(T_，t+u_iT_，i)＝Φ+λT_，ii

式中，为剪切应变速率，下标“’”表示偏导，i，j＝1，2，3为空间局部笛卡儿坐标系分量，u_i为速度分量，p为压力，η为粘度，T为温度，t为时间，ρ、C_P分别为塑料熔体的密度、比热容，λ为塑料熔体的热传导系数，为剪切耗散功，是塑料熔体充模过程的热源项；

根据相邻节点控制体公共界面上流动剪切应力相等的原则，改进传统的有限体积法，得到某节点控制体中心P₀的速度梯度计算式为：

${\▿u}_{P_0}=D^{-1}\·c$

式中， $D=\underset{k=1}{\overset{n_f}{\mathrm{\Σ}}}{(r_{P_k}-r_{P_0})}^T(r_{P_k}-r_{P_0}),c=\underset{k=1}{\overset{n_f}{\mathrm{\Σ}}}\frac{{(r_{P_k}-r_{P_0})}^T(u_{P_k}-u_{P_0})}{1-w_k+q_k\×w_k},$ u为节点速度矢量。如图4所示，下标“k”表示相邻节点控制体，P_k为相邻节点控制体中心，n_f为相邻节点控制体数，r为控制体中心的位置矢量，上标“T”表示转置向量，为加权系数，r_k为与相邻节点控制体的公共界面中心的位置矢量，为两个节点控制体中心的粘度比。

与相邻节点控制体的公共界面中心的速度梯度计算式为：

${\▿u}^{*}=\stackrel{\‾}{\▿u}+[\frac{(u_{P_k}-u_{P_0})}{\left|d_k\right|(1-w_k+q_k\×w_k)}-\frac{\stackrel{\‾}{\▿u}\·d_k}{\left|d_k\right|}]\frac{\left|d_k\right|s_k}{d_k\·s_k}$

式中，u^*为公共界面中心速度矢量， $\stackrel{\‾}{\▿u}=w_k{\▿u}_{P_0}+\frac{(1-w_k){\▿u}_{P_k}}{q_k},d_k=r_k-r_{P_0},$ s_k为相邻节点控制体公共界面的面积外法矢量(由P₀指向P_k)。

与相邻节点控制体的公共界面中心的速度计算式为：

$u^{*}=\stackrel{\‾}{u}+\frac{\stackrel{\‾}{V}}{a}[-\frac{(p_{P_k}-p_{P_0})}{\left|d_k\right|}+\frac{\stackrel{\‾}{{\▿p}_k}\·d_k}{\left|d_k\right|}]\frac{\left|d_k\right|s_k}{d_k\·s_k}$

式中， $\stackrel{\‾}{u}=\frac{q_k{w}_k{u}_{P_0}+w_k{u}_{P_k}}{1-w_k+q_k\×w_k},\frac{\stackrel{\‾}{V}}{a}=\frac{q_k{w}_k\frac{V_{P_0}}{a_{P_0}}+w_k\frac{V_{P_k}}{a_{P_k}}}{1-w_k+q_k\×w_k},$ 其中V为节点控制体体积，a为节点速度分量对应方程组矩阵的主对角线上的值。

与相邻控制体的公共界面中心的压力为：

$p^{*}=0.5\×(p_{P_0}+p_{P_k})$

根据有限体积法建立每个节点控制体的离散的动量守恒方程和质量守恒方程：

$\underset{k=1}{\overset{n_f}{\mathrm{\Σ}}}{\mathrm{\η}}_{P_0}[{\▿u}^{*}+{\left({\▿u}^{*}\right)}^T]\·s_k-\underset{k=1}{\overset{n_f}{\mathrm{\Σ}}}{p}^{*}\·s_k=0$

$\underset{k=1}{\overset{n_f}{\mathrm{\Σ}}}{u}^{*}\·s_k=0$

式中，为的转置张量，下标“k”表示相邻节点控制体，n_f为相邻节点控制体数，由此可得速度方程组矩阵的对角线上的值为：

$a_{P_0}=\underset{k=1}{\overset{n_f}{\mathrm{\Σ}}}\frac{{\mathrm{\η}}_{P_0}}{(1-w_k+q_k\×w_k)}\frac{s_k\·s_k}{d_k\·s_k}$

根据有限体积法建立每个节点控制体的离散的能量守恒方程：

${\∫}_V{\mathrm{\ρC}}_p\frac{T^{t+\mathrm{\Δt}}-T^t}{\mathrm{\Δt}}\mathrm{dV}+\underset{k=1}{\overset{n_f}{\mathrm{\Σ}}}{\mathrm{\ρC}}_p{m}_k{T}^{*}-\underset{k=1}{\overset{n_f}{\mathrm{\Σ}}}{\mathrm{\λ}\▿T}^{*}\·s_k={\∫}_V\mathrm{\ΦdV}$

式中，V为节点控制体的体积，T^t+Δt为t+Δt时刻的温度，T^t为t时刻的温度，Δt为时间步长。上述离散的能量守恒方程中左边第二项为对流项，第三项为扩散项。

为提高计算精度和稳定性，对流项采用“延迟”算法，即相邻节点控制体公共界面中心温度为：

T^*＝T^SU+γ(T^SO-T^SU)

式中，T^SU为根据一阶精度的迎风格式得到的温度，T^SO为采用中心差分的得到的温度，具有二阶精度，即γ为0.8-95之间的常数。m为相邻节点控制体公共界面上的流量：m_k＝u^*·s_k，根据m_k确定T^SO的值：

$T^{\mathrm{SO}}=\left\{\begin{array}{cc}{T}_{P_0}& m_k\&GreaterEqual;0\\ T_{P_k}& m_k<0\end{array}\right.$

扩散项中相邻节点控制体公共界面中心的温度梯度计算式为：

${\&dtri;T}^{*}=\stackrel{\&OverBar;}{\&dtri;T}+[\frac{(T_{P_k}-T_{P_0})}{\left|d_k\right|}-\frac{\stackrel{\&OverBar;}{\&dtri;T}\&CenterDot;d_k}{\left|d_k\right|}]\frac{\left|d_k\right|s_k}{d_k\&CenterDot;s_k}$

式中， $\stackrel{\&OverBar;}{\mathrm{\&Delta;T}}=0.5\&times;({\&dtri;T}_{P_0}+{\&dtri;T}_{P_k}).$

采用节点型有限体积法耦合求解塑料熔体的动量守恒方程、质量守恒方程和能量守恒方程得到速度、压力和温度。速度、压力和温度之间通过粘度互相影响，采用非线性耦合迭代的方法求解速度、压力和温度直到粘度、速度、压力和温度迭代变化在指定的很小范围内；

(5)更新塑料熔体流动前沿，确定时间步长步骤。根据计算得到的速度场，采用流动网络分析法，选择当前流动前沿中各节点控制体所需的最短填充时间作为时间步长，更新流动前沿中其它节点控制体的填充率；

(6)塑料熔体填充量是否达到指定体积，如果达到，进入步骤(7)，否则进入步骤(4)，塑料熔体填充型腔80％时的流动前沿如图5所示；

(7)计算塑料熔体和气体的流动速度、压力和温度步骤。气体穿透推动塑料熔体继续填充模具型腔，塑料熔体和气体之间存在动量交换和质量交换。气辅注塑成型过程中，气体穿透只考虑其动量守恒和质量守恒，忽略惯性力、体积力及可压缩性。该简化在不损失计算精度的前提下，大大减少计算量，提高计算速度和稳定性。气体穿透过程满足动量守恒、质量守恒，对应的控制方程分别为：

$-p_{,i}+2\mathrm{\&eta;}{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{\mathrm{ij},j}=0$

u_i，j＝0

式中，为剪切应变速率，下标“’”表示偏导，i，j＝1，2，3为空间局部笛卡儿坐标系分量，u_i为速度分量，p为压力，η为粘度。

采用节点型有限体积法求解耦合的塑料熔体和气体的动量守恒方程、质量守恒方程和能量守恒方程，忽略气体温度的变化，将气体作为具有一定粘度的假想“塑料熔体”，气体与塑料熔体粘度差异较大，所有节点控制体中心速度梯度、相邻节点控制体公共界面速度和速度梯度的计算方式与上述塑料熔体的速度、速度梯度计算方式相同；

采用节点型有限体积法求解耦合的塑料熔体和气体的动量守恒方程、质量守恒方程和能量守恒方程得到速度、压力和温度。速度、压力和温度之间通过粘度、界面互相影响，采用非线性耦合迭代的方法求解速度、压力和温度直到粘度、速度、压力和温度迭代变化在指定的很小范围内；

(8)更新塑料熔体流动前沿和气体穿透前沿，确定时间步长步骤。根据计算得到的速度场，采用流动网络分析法，选择当前塑料熔体流动前沿和气体穿透前沿中各节点控制体所需的最短填充时间作为时间步长，更新流动前沿中节点控制体的塑料熔体填充率和气体填充率；

(9)型腔是否充满，如果模具型腔已经填充满，则计算结束，否则进入步骤(7)，气体穿透推动熔体充满型腔时的气体穿透前沿如图6所示。

与有限元法相比，有限体积法更适用于气辅注塑成型过程中的流体数值模拟，与传统有限体积法相比，本发明通过实体网格构造三维节点控制体作为有限体积，速度、压力、温度等变量位于节点处，考虑到塑料熔体和气体耦合流动中相邻节点控制体的粘度较大差异，根据相邻节点控制体公共界面上流动剪切应力相等的原则，改进传统的有限体积法，如节点控制体中心的速度梯度、节点控制体公共界面的速度和速度梯度的计算方法，由此提高有限体积法模拟空间粘度差异较大的GAIM过程的数值精度和数值稳定性，避免数值震荡。能量守恒方程的计算中，相邻节点控制体公共界面中心温度采用“延迟”算法提高计算精度，保证了能量守恒方程计算的精度和稳定性。以上措施能更准确地预测塑料熔体流动和气体穿透前沿及各种变量的分布。本发明引入合理的假设与简化，将气体作为假想“塑料熔体”，塑料熔体和气体只是具有不同粘度的流体，从而可以将塑料熔体和气体作为统一的流体进行处理，且忽略流体的惯性力、体积力、可压缩性，由此既可以准确地预测GAIM中气体穿透前沿，还大大减少数值模拟的计算量，可以为合理地设计GAIM产品和设置工艺参数提供科学依据，避免产品缺陷的产生，从而提高产品质量和生产效率。

以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变形或修改，这并不影响本发明的实质内容。

Claims (8)

1.一种基于有限体积法的气辅注塑成型气体穿透预测方法，其特征在于，包括如下步骤：

A、由型腔实体网格构造三维节点控制体作为有限体积；

B、设置成型工艺参数；

C、设置材料参数；

D、计算塑料熔体流动速度、压力和温度；

E、更新塑料熔体流动前沿，确定时间步长；

F、塑料熔体填充量是否达到指定体积，如果达到，进入步骤G，否则进入步骤D；

G、计算塑料熔体和气体的流动速度、压力和温度；

H、更新塑料熔体流动前沿和气体穿透前沿，确定时间步长；

I、模具型腔是否充满，如果模具型腔已经填充满，则计算结束，否则进入步骤G；

所述步骤D中，忽略塑料熔体的惯性力、体积力，塑料熔体为不可压缩流体，其速度、压力、温度变量位于节点控制体中心，所述塑料熔体流动的速度、压力和温度的计算具体包括如下步骤：

(1)根据7参数WLF-Cross粘度模型计算节点控制体中心的粘度，塑料熔体粘度为温度、压力、剪切应变速率的函数，如下：

$\begin{array}{cc}\mathrm{\&eta;}(T,\stackrel{\&CenterDot;}{\mathrm{\&gamma;}})=\frac{{\mathrm{\&eta;}}_0(T,p)}{1+{\left(\frac{{\mathrm{\&eta;}}_0\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}{\stackrel{\&CenterDot;}{\mathrm{\&tau;}}}\right)}^{1-\tilde{n}}},& {\mathrm{\&eta;}}_0=D_1\exp \left(\frac{-A_1(T-(D_2+D_{3}p))}{A_2+T-D_2}\right)\end{array}$

式中，η为粘度，η₀为零剪切粘度，T为温度，p为压力，为非牛顿指数，为剪切应变速率，为材料剪切常数，A₁、A₂、D₁、D₂、D₃为材料参数，其中T、p、为节点控制体中心的变量值；

(2)塑料熔体速度、压力、温度的控制方程分别采用如下动量守恒方程、质量守恒方程和能量守恒方程：

$-p_{,i}+2\mathrm{\&eta;}{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{ij,j}=0$

u_i,i＝0

ρC_p(T_，t+u_iT_，i)＝Φ+λT_，ii

式中，为剪切应变速率，下标“,”表示偏导，i,j＝1,2,3为空间局部笛卡儿坐标系分量，u_i为速度分量，p为压力，η为粘度，T为温度，t为时间，ρ、C_P分别为塑料熔体的密度、比热容，λ为塑料熔体的热传导系数，为剪切耗散功，是塑料熔体充模过程的热源项；

(3)采用有限体积法对动量守恒方程进行离散，

a、某节点控制体中心P₀的速度梯度计算式为：

$\&dtri;u_{P_0}=D^{-1}\&CenterDot;c$

式中，u为节点速度矢量，下标“k”表示相邻节点控制体，P_k为相邻节点控制体中心，n_f为相邻节点控制体数，r为控制体中心的位置矢量，上标“T”表示转置向量，为加权系数，r_k为与相邻节点控制体的公共界面中心的位置矢量，为两个节点控制体中心的粘度比；

b、与相邻节点控制体的公共界面中心的速度梯度计算式为：

$\&dtri;u^{*}=\stackrel{\&OverBar;}{\&dtri;u}+\&lsqb;\frac{(u_{P_k}-u_{P_0})}{\left|d_k\right|(1-w_k+q_k\&times;w_k)}-\frac{\stackrel{\&OverBar;}{\&dtri;u}\&CenterDot;d_k}{\left|d_k\right|}\&rsqb;\frac{\left|d_k\right|s_k}{d_k\&CenterDot;s_k}$

式中，u^*为公共界面中心速度矢量，s_k为相邻节点控制体公共界面的面积外法矢量，是由P₀指向P_k；

c、与相邻节点控制体的公共界面中心的速度计算式为：

$u^{*}=\stackrel{\&OverBar;}{u}+\stackrel{\&OverBar;}{\frac{V}{a}}\&lsqb;-\frac{(p_{P_k}-p_{P_0})}{\left|d_k\right|}+\frac{\stackrel{\&OverBar;}{\&dtri;p_k}\&CenterDot;d_k}{\left|d_k\right|}\&rsqb;\frac{\left|d_k\right|s_k}{d_k\&CenterDot;s_k}$

式中，其中V为节点控制体体积，a为节点速度分量对应方程组矩阵的主对角线上的值；

d、采用有限体积法对能量守恒方程进行离散，界面中心的温度计算式为：

T^*＝T^SU+γ(T^SO-T^SU)

式中，T^SU为根据一阶精度的迎风格式得到的温度，T^SO为采用中心差分的得到的温度，具有二阶精度，即γ为0.8-95之间的常数；

e、采用节点型有限体积法耦合求解所述塑料熔体的动量守恒方程、质量守恒方程和能量守恒方程得到速度、压力和温度。

2.如权利要求1所述的基于有限体积法的气辅注塑成型气体穿透预测方法，其特征在于，所述步骤A具体为：导入型腔实体网格，连接网格单元的体心、面心、边中点构造三维节点控制体作为有限体积。

3.如权利要求1所述的基于有限体积法的气辅注塑成型气体穿透预测方法，其特征在于，步骤B中，所述成型工艺参数包括塑料熔体的注射压力、注射时间或注射速度，气体的注射压力、注射时间或注射速度，以及模具温度、熔体温度。

4.如权利要求1所述的基于有限体积法的气辅注塑成型气体穿透预测方法，其特征在于，步骤C中，所述材料参数包括气体粘度，塑料熔体的密度、比热容、热传导系数，以及塑料熔体的7参数WLF-Cross粘度模型的材料常数。

5.如权利要求4所述的基于有限体积法的气辅注塑成型气体穿透预测方法，其特征在于，所述材料常数包括A₁、A₂、D₁、D₂、D₃。

6.如权利要求1所述的基于有限体积法的气辅注塑成型气体穿透预测方法，其特征在于，所述步骤E具体为：根据计算得到的速度场，采用流动网络分析法，选择当前流动前沿中各节点控制体所需的最短填充时间作为时间步长，更新流动前沿中其它节点控制体的填充率。

7.如权利要求1所述的基于有限体积法的气辅注塑成型气体穿透预测方法，其特征在于，所述步骤G中，气体穿透推动塑料熔体继续填充模具型腔，将气体视为假想塑料熔体，塑料熔体和气体之间存在动量交换和质量交换，所述塑料熔体和气体的流动速度、压力和温度的计算具体包括如下步骤：

a、气辅注塑成型过程中，气体穿透只考虑动量守恒和质量守恒，忽略惯性力、体积力及可压缩性；对应的动量守恒、质量守恒控制方程分别为：

$0=-p_{,i}+2\mathrm{\&eta;}{\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{ij,j}$

u_i,i＝0

式中，为剪切应变速率，下标“,”表示偏导，i,j＝1,2,3为空间局部笛卡儿坐标系分量，u_i为速度分量，p为压力，η为气体粘度；

b、忽略气体温度的变化，气体的节点控制体中心速度梯度、相邻节点控制体公共界面中心速度和速度梯度的计算方式同所述塑料熔体的节点控制体中心速度梯度、相邻节点控制体公共界面中心速度和速度梯度的计算方式；

c、采用节点型有限体积法求解耦合的塑料熔体和气体的动量守恒方程、质量守恒方程和能量守恒方程得到速度、压力和温度。

8.如权利要求1所述的基于有限体积法的气辅注塑成型气体穿透预测方法，其特征在于，所述步骤H具体为：根据计算得到的速度场，采用流动网络分析法，选择当前塑料熔体流动前沿和气体穿透前沿中各节点控制体所需的最短填充时间作为时间步长，更新流动前沿中节点控制体的塑料熔体填充率和气体填充率。