塑料制品具有质量轻、易于成型、生产效率高、成本低等诸多优点，被广泛应用于消费品、医疗、汽车、航空航天等各个领域．为进一步提升塑料制品的外观和性能，发展出多种塑料注射成形与功能结构相结合的复合成形工艺．如模内层压工艺[1]可将皮革等材料与塑料注射成形结合，从而获得免粘接的汽车内饰件；连续纤维复合材料注射压缩-热压成形工艺[2-4]可一步成形具有连续纤维增强结构的复杂塑料制品．上述复合成形工艺存在塑料熔体与功能结构的复杂双向耦合问题．功能结构在塑料熔体压力的作用下发生大变形行为，同时功能结构的变形又会改变塑料熔体的充填区域，这一现象使得功能结构与塑料熔体之间存在一个柔性移动边界．柔性边界的移动使得塑料熔体的计算网格随之变形，且每一时刻的网格边界均须动态计算获得，因此充填模拟须从初始网格逐步变形至终态网格，显著区别于注射压缩成形模拟的网格变形方式[5]，数值模拟难度也更大．网格变形导致其质量下降，从而引起塑料熔体充填模拟的稳定性问题．当网格发生严重畸变时，会产生网格单元重叠、反转等错误，导致数值模拟失败．针对塑料熔体充填区域显著改变而引起的充填模拟稳定性问题，本研究基于四面体网格，采用任意拉格朗日欧拉法求解大变形问题和捕捉移动边界，使用四面体单元中引入中心速度稳定的速度压力耦合有限元计算方法和网格动态优化的反转单元罚函数方法，实现此类成形工艺中塑料熔体充填的稳定模拟．1 移动边界下的充填模拟数学模型本研究以连续纤维复合材料注射压缩-热压成形为例，通过分析此类工艺的成形特点，建立起塑料熔体充填模拟的数学模型．如图1所示，该工艺过程分为如下三个阶段：a．注射阶段，塑料熔体通过流道系统注入型腔；b．压缩阶段，模具合模对塑料熔体和复合板材进行压缩成形；c．保压冷却阶段，塑料熔体和复合板材固化后得到最终制品．10.13245/j.hust.230620.F001图1注射压缩-热压成形图2(a)和(b)分别为初始状态和结束状态下模具、塑料熔体与复合板材的位置．图2中可见模具闭合使得塑料熔体向型腔四周扩展，由于塑料熔体的压力较大，因此导致复合板材同步发生大变形．塑料熔体与动模、复合板材的接触界面随时间变化，导致其计算域发生动态更新，因此本研究使用任意拉格朗日欧拉法进行建模．基于任意拉格朗日欧拉法的计算域如图2中红色虚线区域所示，在计算域上应用类似于欧拉法的控制方程，因此可以处理塑料熔体的大变形问题．同时移动边界上的网格又可像拉格朗日法一样变形，使得塑料熔体的网格能与动模和复合板材的网格协调变形，这种网格协调变形的能力使得任意拉格朗日欧拉法能够对柔性移动边界的运动进行精确捕捉．10.13245/j.hust.230620.F002图2任意拉格朗日欧拉法本研究使用基于任意拉格朗日欧拉法的连续性方程式、动量方程式和能量方程式来定义塑料熔体在型腔中的充填过程[6-7]．不同于欧拉法中的控制方程，由于任意拉格朗日欧拉法中存在网格移动速度，因此下式的对流项中出现了对流速度矢量c，以体现网格速度对流体对流造成的影响．即∂ρ/∂t+c⋅∇ρ=-ρ∇⋅v；(1)ρ(∂v/∂t+(c⋅∇)v)=-∇p+∇⋅τ+ρg；(2)     ρcp(∂T/∂t+(c⋅∇)T)=∇⋅(k∇T)+ηγ2+βT(dp/dt), (3)式中：∇⋅为散度算子；∇为梯度算子；v为材料速度矢量；c为对流速度矢量；g为重力加速度矢量；τ为偏应力张量；p为流体压力；T为温度；t为时间；ρ为材料密度；cp为材料比热容；k为材料导热系数；η为材料动力学黏度；γ为剪切速率．与欧拉法相比，任意拉格朗日欧拉法的数值实现须要额外对网格变形和对流项进行处理，并在变形后的网格上求解速度压力等物理场量．2 速度压力场耦合计算方法求解式(1)～(3)最大的难点是速度场和压力场的计算．基于热塑性塑料熔体充填过程中的蠕流假设[7]和不可压缩假设[8]，式(1)和式(2)可进一步简化为：∇⋅v=0；(4)∇⋅τ-∇p=0．(5)塑料熔体的速度场和压力场通过求解式(4)和式(5)获得，其数值计算方法分为迭代法和耦合法．迭代法通过数学处理将压力与速度解耦，然后使用数值迭代的方法来逼近压力和速度的真实解．迭代法的优点是节省内存，缺点为速度压力多次迭代的求解过程较长．迭代过程中一般须使用亚松驰方法来提升算法的收敛性，但不合适的亚松驰因子取值会进一步恶化收敛速度，甚至导致数值发散．耦合法则是将连续性方程与动量方程联立为一个整体方程组，直接求解出速度场和压力场，由于无须速度-压力迭代过程，因此其数值稳定性和精度均显著优于迭代法．当使用耦合有限元法求解速度压力方程时，若直接使用线性四面体单元，并对速度自由度和压力自由度使用相同的插值函数，则可得到如下系统矩阵ABTB0vp=00，式中A和B为速度压力方程的系数．压力自由度系数为零，系统矩阵不满足极限上确界(Inf-Sup)条件[9]，导致压力解的发散与震荡．一种解决该问题的方法是为速度自由度和压力自由度选择不同的插值函数，例如对于四面体单元，速度使用二阶插值函数而压力使用线性插值函数．然而高阶有限元单元在计算流体力学中并不是一个好的选择，这些单元在数值计算上并不稳定，同时也会显著增加单元的速度自由度和计算量．鉴于以上问题，本研究使用引入中心速度的MINI四面体单元[10]来稳定压力场．MINI单元不仅满足Inf-Sup条件，且在最终系统矩阵中也并未增加单元节点自由度，仍保持了线性四面体单元的便利性与简洁性．不同于线性四面体单元，四面体MINI单元中心还存在一个额外的速度自由度，因此MINI单元中速度可表示为v=vl+vb，式中：vl为线性速度；vb为气泡速度．使用MINI单元形函数和单元节点速度自由度表示，v=[N1N2N3N4Nb]⋅[v1v2v3v4vb]T;p=[N1N2N3N4][p1p2p3p4]T，式中：N1，N2，N3和N4为线性四面体单元节点形函数；Nb=2-4∑k=14Nk2为气泡速度自由度形函数．使用伽辽金法和MINI单元对式(4)和式(5)进行有限元离散，可得如下速度压力系统矩阵K0L0KbLbLTLbT0vlvbp=000．(6)系统矩阵中各项权重系数如下K=  2k11+k22+k33k21k31k12k11+2k22+k33k32k13k23k11+k22+2k33，kij=η∫V∂NT∂xi∂N∂xjdV；Kb=W1kb21kb31kb12W2kb32kb13kb23W3，kbij=16V5∑k=14∑l=141+δkl∂Nk∂xi∂Nl∂xj；LT=[L1L2L3]，Li=-∫V∂NT∂xiNdV；LbT=[Lb1Lb2Lb3]，Lbi=-45∫V∑k=14∂Nk∂xidV，式中：W1=2kb11+kb22+kb33；W2=kb11+2kb22+kb33；W3=kb11+kb22+2kb33；V为体积．式(6)可做进一步简化消去气泡速度自由度vb，并得到最终单元系统矩阵KLLT-LbKb-1Lbvlp=00．上式可见单元节点自由度仍然只存在线性速度自由度和压力自由度，与线性四面体单元一致，但压力自由度的矩阵系数不再为零，而是-LbKb-1Lb，从而使得MINI单元满足Inf-Sup条件，避免了数值振荡的问题．最终获得的系统矩阵为对称矩阵，可使用共轭梯度法高效求解，从而获得流场的速度与压力．塑料熔体在充填过程中不断向外扩展，直至充满型腔．在此过程中，须要在每一时刻下更新当前流体的充填区域，然后在更新的充填区域之上进行速度场、压力场的求解．本研究采用基于水平集方法[11]的自由表面流模型来进行塑料熔体流动前沿的更新与追踪．3 动态网格失效问题与优化算法任意拉格朗日欧拉法的优点是能求解大变形问题的同时又能精确捕捉移动边界，这是通过调整计算网格来实现的．移动边界的运动会导致塑料熔体网格的位置发生改变，为避免仅移动边界网格而导致塑料熔体内部网格产生畸变，塑料熔体整体计算域的网格均须随着边界网格的移动而相应调整．塑料熔体网格的重构可通过重新划分网格，并将旧网格上的物理量映射至新网格上实现，但每一步均进行网格重新剖分会带来巨大的计算量和时间开销，相比之下，网格变形[12]则是更加可行的方案．无限插值法将移动边界上节点位移等比分配至网格内部节点，从而实现网格移动；反向距离加权插值法中网格内部节点位移等于其相邻节点位移的反向距离加权平均；弹簧模型使用胡克定理计算网格内部节点位移；拉普拉斯模型通过求解拉普拉斯方程来确定内部节点位移，节点位移的各个分量彼此独立无关；线弹性模型使用有限元法求解线弹性材料的变形以获得内部网格节点的位移．上述模型均不具备显式抑制网格变形过程中反转四面体单元产生的能力，因此在网格大变形的情况下，无法避免反转单元的产生．反转单元具有负的单元体积，从而引起数值稳定性问题，导致有限元法的计算失败．由于移动边界上节点位移已知，可将其作为狄利克雷边界条件，使用线弹性网格变形模型求出内部节点的位移x．线弹性网格变形模型具有如下能量方程E(x)=12∫Vσ:εdV，式中：E(x)为应变能；σ为应力张量；ε为应变张量；V为单元体积．线弹性模型的本构方程为σ=2με+λtr(ε)I，式中：μ和λ为拉梅常数；tr(ε)为应变的迹；I为二阶单位张量．由于网格变形中并不需要真实的物理材料，因此取μ=1和λ=1．进一步，已知反转单元导致负体积的四面体单元，由连续介质力学可知，单元变形前后体积的映射关系可以表述为dV=det(F)dV0，式中：dV为单元变形后的体积；dV0为变形前的体积；det(F)为单元变形梯度的行列式．变形梯度张量定义如下Fij=∂xi/∂Xj，式中：xi为当前构形；Xj为初始构形．假设单元变形前为非反转单元，欲使得变形后单元也为非反转单元，则必须det(F)0，进一步可得变形梯度张量的对称分量的行列式det(S)0，且二阶对称张量S的行列式与其特征值存在以下关系det(S)=∏i=1nλi，式中λi为S的第i个特征值．从上式可以得出满足det(S)0的一个充分条件是二阶对称张量S所有的特征值λi0．基于上述条件，本研究引入以下反转单元罚函数项[13]以避免变形过程中负体积单元的产生，Epen(x)=∫V∑i=1np(λi)dV．仅当特征值λ0时，定义p(λ)=λ2使得反转单元的罚函数项非零，否则罚函数项为零．因此引入反转单元罚函数的线弹性网格变形模型的能量方程如下式所示E(x)=12∫Vσ:ε+τ∑i=1np(λi)dV．(7)上式可通过牛顿-拉夫森方法求解，其海森矩阵具备稀疏正定性质，因此可使用预条件共轭梯度法进行求解．τ为罚函数权重系数，迭代中逐步增加τ值，直至罚函数能量为零，网格中所有反转单元均被修正．求解得到新的网格坐标之后，即可计算出当前时间步下的网格速度v^，进而获得式(1)～(3)中的对流速度c=v-v^．4 模拟方法的稳定性验证4.1　大变形下网格优化方法的稳定性验证本研究设计了网格压缩变形测试，实验模型为一个高度为10 mm、直径为4 mm的圆柱体．将该圆柱模型进行四面体网格剖分，网格尺寸为0.2 mm．测试中将圆柱体的底端固定，然后将圆柱体的顶端向下压缩9 mm，压缩比达90%．该圆柱体压缩实验的网格如图3所示．为验证反转单元罚函数项的作用，进行了两组网格变形实验，方案一将式(7)中的罚函数项去除，方案二则使用完整的式(7)进行计算．10.13245/j.hust.230620.F003图3网格压缩测试算例变形后的圆柱网格形状如图4所示，图4(a)为使用实验方案一获得的网格，可见变形后网格出现了严重拓扑错误，且同时产生了5 628个负体积单元，如图中红色区域所示．右图4(b)为使用方案二获得的网格，变形后网格完全符合实验预期，且不存在任何错误的四面体单元．10.13245/j.hust.230620.F004图4网格变形结果对比图5为压缩实验中罚函数权重系数τ与罚函数能量E的曲线．当罚函数权重系数τ取值较小时(如1×10-5)，系统的罚函数值为243 kJ，这意味着未经优化的四面体网格存在大量错误，印证了图4(a)中完全错误的圆柱体压缩网格结果．随着罚函数权重系数τ的增大，罚函数能量趋近于0，最终网格不存在错误单元．实验结果证明在网格极端压缩的情况下，本研究提出的网格优化算法仍然能够稳定地实现网格变形，避免错误单元的产生．10.13245/j.hust.230620.F005图5压缩实验罚函数权重系数τ与罚函数能量E的曲线4.2　速度压力耦合方案的稳定性验证网格变形优化算法保证了网格的正确性，但网格压缩仍会导致有限元单元的质量下降，引起充填模拟的数值稳定性和精度问题．塑料熔体压缩过程中，单元厚度不断变小，导致单元形态恶化．当单元厚度压缩至零或反转时，会导致数值模拟失败．压缩过程中四面体网格质量可使用纵横比进行评价，其定义为单元最大边长与单元最小高度的比值．网格厚度压缩，导致纵横比过大，单元形态变差．当单元厚度压缩至零时，纵横比为无穷大．为验证本研究充填模拟算法在不同网格质量下的数值稳定性，设计了不同纵横比下的三维薄壁平板充填模拟实验．通过对平板模型厚度方向增加四面体单元层数来模拟单元纵横比恶化的情况，并通过对不同网格下压力曲线对比，验证充填模拟算法的稳定性．实验平板模型如图6所示，材料黏度为200 Pa∙s，密度为1 000 kg/m3，充填时间为2 s．10.13245/j.hust.230620.F006图6平板充填模型(mm)表1给出了不同网格方案下的网格纵横比，其中方案1中网格平均纵横比为20.6，四面体网格形态较好；方案3中网格平均纵横比接近方案1的4倍，且最大纵横比为277.8，网格被极度压缩．10.13245/j.hust.230620.T001表1网格方案网格方案平均纵横比最大纵横比方案120.690.5方案241.2167.2方案379.5277.8不同网格方案下压力曲线如图7(a)所示，三组方案下压力均值与方差如图7(b)所示．根据图中结果分析可知，本文提出的塑料熔体充填模拟算法在三种网格方案下均能成功进行求解，且压力结果受到单元质量的影响不显著，其最大压力标准差为1 MPa，所有时间步下平均标准差为0.34 MPa，从而证明了该充填模拟方法在网格畸变情况下的数值稳定性．10.13245/j.hust.230620.F007图7不同网格方案下压力稳定性分析5 实验验证在此通过对碳纤维复合材料注射压缩-热压成形工艺实验进行数值模拟，从而验证本研究提出的模拟方法的精确性和有效性．5.1　实验设计图8为实验用的注塑机、模具和制品．注塑机型号为JSW J110ADC-180，热塑性压缩材料为聚丙烯．实验过程中熔体温度为220 ℃，复合板材预热温度为180 ℃，模具温度为40 ℃；压缩距离为25 mm，压缩速度为25 mm/s；保压时间为2 s，冷却时间为50 s；初始坯料质量为25 g．复合板材使用T300碳纤维及聚丙烯热塑性基体，在此采用文献[14]的材料模型，其模型参数如下：C11=230 GPa，γlim=50°；当纤维剪切角γγlim时，C33=0.054 MPa；当纤维剪切角γ≥γlim时，C33=1.8 MPa．10.13245/j.hust.230620.F008图8注射压缩-热压实验5.2　数值模拟与验证使用提出的方法，通过自行编写的程序对热塑性塑料熔体成形过程进行模拟，并使用开源后处理软件ParaView对模拟结果进行展示．碳纤维复合板材的热压成形使用Abaqus进行模拟[4]，软件在每一个时间步下动态调用Abaqus以实现塑料熔体与复合板材的双向耦合计算．图9展示了在该工艺模拟过程中任意拉格朗日欧拉网格的变形情况，从图9中可看出：任意拉格朗日欧拉网格在成形过程中被显著压缩，此过程中通常会产生大量扭曲和错误的有限元单元，进而造成数值计算失败．10.13245/j.hust.230620.F009图9网格变形结果由于反转四面体单元罚函数能量只在网格全部正确的情况下为零，因此罚函数能量可用作网格质量评价标准，罚函数能量越小，网格质量越好．图10为不同时间步(t)下，线弹性网格变形方法和本文网格变形方法对应的罚函数能量E．从图10可见：在网格变形初期两种算法的罚函数能量均较小，随着网格变形程度的逐渐增大，在1.0 s时刻线弹性网格变形算法的罚函数能量是本文网格变形方法的1×104倍，从而证明在网格极度压缩的情况下，线弹性网格变形算法失效，而本文网格算法仍能保持正确的网格形态．10.13245/j.hust.230620.F010图10不同阶段网格变形罚函数能量图11为在不同压缩行程下塑料熔体充填实验结果与模拟结果的对比[4]，可见随着压缩行程的增加，塑料逐渐扩展并充满整个型腔，且同时伴有复合板材的变形．10.13245/j.hust.230620.F011图11充填结果对比由于该制品为圆盘形的充填模式，本研究使用图11所标注流动直径D10，D15和D25来进行模拟的流动前沿与实际的流动前沿对比，其中D10，D15和D25分别代表压缩行程为10，15和25 mm时塑料熔体的流动直径，从而验证模拟结果的精度．本研究通过对实际制品与模拟结果上的流动直径进行测量，结果如表2所示[4]，可见模拟结果的流动前沿与实验结果的流动前沿符合良好，最大尺寸误差为3.3%．10.13245/j.hust.230620.T002表2不同压缩行程下流动前沿的结果对比压缩行程/mm模拟值/mm实验值/mm误差/%1087.190.03.215114.9118.73.225124.0120.03.36 结语本研究针对塑料熔体充填区域显著改变引起的充填模拟稳定性问题，提出了稳定的速度压力耦合计算方案和引入反转单元罚函数的网格优化方法，实现了移动网格下塑料熔体充填的稳定高效模拟．稳定性验证发现在网格极端变形和不同单元纵横比下，本研究提出的模拟方法均能进行稳定求解．设计了碳纤维复合材料注射压缩-热压成形工艺实验，并根据成形工艺参数进行数值模拟，发现移动边界下塑料熔体流动前沿的模拟误差为3.3%，表明本文提出的移动边界条件下塑料熔体充填的稳定模拟算法能够稳定精确地对成形工艺过程进行模拟，对实际生产有指导意义．