基于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

為了準確求解熔體在型腔內的黏彈性充填，假設在型腔模壁上為無滑移無滲透邊界，法向偏應力張量梯度為零，入口速度為恒定值，入口偏應力張量為零。假設充填過程中，型腔排氣良好，即出口處速度和偏應力張力的法向梯度均為零為零，壓力恒為大氣壓。對于溫度，入口溫度即為熔體的加工溫度，模壁和流動前沿溫度采用第三類邊界條件：

其中，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）本文所發展算法可以準確預測注射成型過程中熔體的流動前沿的形狀和位置、澆口壓力大小、熔接線位置和形狀，并且可以用之計算充填過程中熔體的流動誘導應力。