基于可壓縮兩相流模型的注射成型填充模擬

王振海，賈偉，董勤喜，于行

[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]對比?？梢园l現，在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

3 結論

(1) 基于有限體積法以及可壓縮兩相流模型，開發用于注塑模擬的三維數值計算模型。采用高分辨率有限體積法求解VOF自由界面追蹤方程，并使用PISO算法計算壓力和速度場以求解三維熔體流動與傳熱控制方程。使用Cross-WLF黏度模型表征聚合物熔體的流變行為。

(2) 對開發的三維數值計算模型進行了數值驗證，與帶有圓柱形嵌件型腔中熔體填充過程的實驗結果相比，數值結果顯示出非常好的一致性。成功地模擬了速度、溫度、剪切速率和壓力在不同截面上的分布，表明與2.5維模擬相比，本三維模擬程序可以提供更準確、更詳細的流動特性信息。

(3) 對不均勻厚度型腔進行了填充模擬，將其模擬結果與商業軟件Moldflow MPI_3D以及HsCAE3D的模擬結果進行了對比，驗證了模擬結果的正確性，表明筆者開發的用于注塑模擬的三維數值計算模型能比較準確地模擬三維注射成型填充過程。