李 勇,鄭建榮

(華東理工大學(xué)機械與動力工程學(xué)院,上海200237)

聚合物熔體矩形流道擠出脹大的數(shù)值模擬

李 勇,鄭建榮

(華東理工大學(xué)機械與動力工程學(xué)院,上海200237)

運用Picard迭代格式的有限元方法,采用Wagner積分型本構(gòu)方程對黏彈流體擠出脹大進行三維模擬分析。每次迭代根據(jù)新的自由面邊界位置重新劃分網(wǎng)格,由前一次迭代得到的速度場,算出單元高斯點流線,沿流線積分計算應(yīng)力,把應(yīng)力作為擬體力,建立非線性方程組迭代格式。對不同寬度和長度的矩形流道的擠出脹大進行模擬計算,分析了流道寬度和長度對脹大率的影響。結(jié)果表明,隨著寬厚比的增加,厚度脹大率隨之增加;隨著流道長度的增加,脹大率逐漸下降。該結(jié)果與二維狹縫流道的數(shù)值模擬和實驗結(jié)果比較表明,用該方法對黏彈流體擠出脹大流動進行三維模擬是可行的和準(zhǔn)確的。

聚合物;熔體;矩形流道;擠出脹大;數(shù)值模擬

0 前言

擠出脹大是黏彈流體在擠出過程中的固有特性,是擠出成型中最被關(guān)心的問題之一。擠出脹大現(xiàn)象給塑料擠出機頭的設(shè)計造成很大困難。擠出脹大的數(shù)值模擬對擠出工藝的控制和提高機頭設(shè)計的品質(zhì)具有重要意義。

聚合物擠出脹大的數(shù)值模擬存在兩大困難,即本構(gòu)方程的選擇和自由面的形狀。由于Maxwell本構(gòu)方程相對簡單,黏彈流體擠出脹大問題的數(shù)值模擬是從研究Maxwell流體的脹大問題開始的,Chang等[1]首先進行了黏彈Maxwell流體擠出脹大的數(shù)值模擬研究。另外,Debae等[2]采用改進的混合有限元法(EVSS/SUPG)模擬,得到收斂解的最高Weissenberg數(shù)是3.75。

但Maxwell方程本構(gòu)關(guān)系無法描述真實黏彈流體的流變特性,積分型本構(gòu)方程包含了松弛時間譜,被認(rèn)為是對復(fù)雜流動綜合預(yù)測能力最強的本構(gòu)關(guān)系,特別是對于擠出脹大問題。而且不會像Maxwell方程那樣遇到高Weissenberg數(shù)的問題。采用積分型本構(gòu)方程關(guān)鍵在于應(yīng)變歷史和應(yīng)力的計算。Luo等[3-4]采用流線有限元法首先對K-BKZ類模型的熔體擠出脹大進行模擬,并在本構(gòu)方程中又引入拉伸系數(shù)譜,準(zhǔn)確擬合了熔體的單軸拉伸黏度。Goublomme等[5]同時用Wagner模型和PSM模型計算了高密度聚乙烯(PEHD)熔體的擠出脹大流動。計算中發(fā)現(xiàn),用Wagner模型更能夠準(zhǔn)確地預(yù)測實驗結(jié)果,而且脹大比隨著模型中表征剪切流動的參數(shù)和表征拉伸流動的參數(shù)的增大而單調(diào)降低。Kiriakidis等[6]采用基于常規(guī)有限單元的應(yīng)變歷史跟蹤方式計算黏彈應(yīng)力,模擬了PE-HD(PSM模型表征)和線形低密度聚乙烯(PE-LLD)熔體經(jīng)過狹縫口模的擠出脹大,結(jié)果表明,口模的長寬比顯著影響?zhàn)椓黧w擠出脹大的結(jié)果,PE-HD熔體的擠出脹大程度隨著流率的增加而增大,計算結(jié)果與實驗值十分吻合,肯定了K-B KZ類模型能夠準(zhǔn)確地預(yù)測黏彈流體在口模內(nèi)的應(yīng)力松弛過程。目前,擠出脹大的研究主要集中在二維流動,對于機頭流道中流動的三維數(shù)值模擬有一定研究[7-9],但對于擠出脹大的三維數(shù)值模擬,由于自由面的確定和網(wǎng)格劃分的復(fù)雜性,處理起來較為困難,研究較少。

在前期的研究中[10],曾采用應(yīng)力增量迭代方法和流線有限元方法,對三維擠出脹大進行了初步研究,本文采用Wanger積分型本構(gòu)方程,運用Picard迭代并根據(jù)矩形流道特點使用自由面流線迭代和內(nèi)部網(wǎng)格重新劃分相結(jié)合的方法對矩形流道的擠出脹大進行了細(xì)致分析,得到了流動的速度場、流線和擠出脹大的三維形狀。

1 控制方程

根據(jù)擠出加工的特點,聚合物熔體等溫、蠕變流動運動方程為:

連續(xù)性方程:

本構(gòu)方程采用Wagner給出的積分型本構(gòu)方程[11],其應(yīng)力由式(3)給出:

式中 T(t)——應(yīng)力張量

p——靜壓力,Pa

τ——偏應(yīng)力張量,Pa

V——速度向量,m/s

I——單位矩陣

m(t-t′)——代表時間依賴性的記憶函數(shù),由材料的黏彈性數(shù)據(jù)決定

C-1——Finger應(yīng)變張量,是Cauchy張量C的逆

I1——C-1的第一不變量

I2——C的第一不變量

h(I1,I2)——非線性的衰減函數(shù)

t——當(dāng)前時間,s

t′——時間的積分變量,s

記憶函數(shù)m(t-t′)是松弛時間的函數(shù),可以近似地表達為不連續(xù)的時間段的迭加:

式中 ai——松弛模量,Pa

λi——松弛時間,s

對于任意變形的流動,要求給出h與I1、I2之間的關(guān)系,定義:

其中α=0.32,有:

對于PE-LD熔體,式(6)中β=0.35,n1=0.31,n2=0.106。

2 數(shù)值方法

2.1 非線性本構(gòu)方程的迭代

熔體本構(gòu)關(guān)系的非線性特性,是高分子流體擠出問題的主要困難之一。主要解決方法分為兩類:第一類是混合法,即把應(yīng)力、速度和壓力作為未知量一起求解;另一類稱分裂格式,即Picard迭代,它僅把速度和壓力作為未知量求解,用前一次迭代的速度場計算非牛頓應(yīng)力,再把應(yīng)力項當(dāng)作擬體力處理。該方法雖然迭代次數(shù)多、但迭代格式實施簡單、成熟,一次求解的未知變量小于混合格式,計算量低,同時易于方便地結(jié)合微分型或積分型本構(gòu)方程以及通過分子模型表達的應(yīng)力關(guān)系進行黏彈流體流動數(shù)值模擬,由于應(yīng)力不能表達成速度的顯函數(shù),所以本文采用后一種方法,以速度矢量(ν)和壓力(p)為未知量,對式(1)、(2)采用Galerkin弱解形式,由于應(yīng)力不能表達成速度的顯函數(shù),所以把應(yīng)力項放在方程的右端,建立迭代格式:

式中 φu——速度的形函數(shù)

φp——壓力的形函數(shù)

?！e分區(qū)域Ω的邊界

ηr——參考黏度,Pa·s

τ(n)——由第n次迭代結(jié)果ν(n)、pn通過本構(gòu)方程計算得到

d——應(yīng)變速率張量,s-1

N——面力,Pa

f——體力,N/m3

2.2 自由邊界的迭代方法

自由邊界的邊界條件為混合邊界條件,即應(yīng)力和法向速度都為零,所以邊界線就是流線。積分型本構(gòu)方程的熔體二維擠出脹大的模擬大都采用流線有限元方法[3-4],即沿著流線劃分網(wǎng)格。它利用自由邊界為流線這一性質(zhì),先假定自由邊界位置,根據(jù)應(yīng)力邊界條件求解,得到速度分布,對速度分量積分,得到流線即新的自由邊界,同時內(nèi)部節(jié)點也沿流線分布。由于每次迭代中計算出的速度結(jié)果并不是最終速度,所以每次迭代都采用速度場計算內(nèi)部節(jié)點流線的意義不大,收斂速度較慢。并且由于三維自由面的形狀相對于二維曲線復(fù)雜得多,當(dāng)速度存在誤差時,如果積分步長過大會造成單元網(wǎng)格畸形。所以,對于三維情況,本文僅對邊界流線進行迭代求解,以計算自由面位置,每次邊界更新后,重新對網(wǎng)格進行自動劃分,確定內(nèi)部節(jié)點,所以內(nèi)部節(jié)點并不沿流線分布。計算中采用了網(wǎng)格自動劃分技術(shù),對矩形流道的單元網(wǎng)格自動劃分。

設(shè)x為流動方向,則新的自由面由式(10)、(11)得到:

2.3 單元內(nèi)部流線的計算

應(yīng)力的計算有2種方法,Tanner沿流線即單元邊界積分計算應(yīng)力,而范毓?jié)欉M行改進,采用沿通過高斯點的流線進行積分的方法。本文由于在每次迭代時單元內(nèi)部的節(jié)點由網(wǎng)格自動劃分得到,所以節(jié)點和高斯點并不沿流線分布,必須首先得出單元高斯點流動軌跡,再沿軌跡積分計算應(yīng)力。根據(jù)速度場,可計算高斯點流線。穩(wěn)態(tài)流動質(zhì)點的加速度可以通過速度表達為:

根據(jù)式(12)可以一步步地給出當(dāng)前時刻t在高斯點的質(zhì)點的流動軌跡,用張量表示為:

式中 x′i——質(zhì)點在時刻t′的坐標(biāo)

式(13)中上標(biāo)1表示在xi1處的值,而xi1則由式(15)給出:

根據(jù)式(13)、(15)可以一步步地向回求出質(zhì)點的位置,公式具有三次精度,另外時間推移量要盡量小,因為整個積分包括很多時間步,如果步長較大,每一步的誤差會產(chǎn)生較大的累計誤差。

3 出口脹大的計算結(jié)果

本文采用上述三維自由面迭代方法,對低密度聚乙烯(PE-LD)在180℃時通過矩形流道擠出的擠出脹大過程進行了仔細(xì)分析。計算中采用20節(jié)點六面體單元,其中速度為20節(jié)點,采用二次插值,壓力為8節(jié)點,采用一次插值。由于計算量較大,計算中最大單元數(shù)為432個。

其中自由邊界的迭代步長不能過大,自由面經(jīng)過大約500次迭代,可以達到較高的收斂精度。用該方法進行擠出脹大的自由面迭代,雖然由于內(nèi)部節(jié)點不在流線上,在計算應(yīng)力時需要根據(jù)速度場計算經(jīng)過高斯點的流線,增加了每次迭代的計算量,但避免了由于局部流線不均勻而可能造成的流體內(nèi)部單元畸形。本構(gòu)方程中的材料松弛材特性如表1所示。

表1 PE-LD在180℃時的流變參數(shù)Tab.1 Rheological parameters of PE-LD at 180℃

在二維狹縫流動中,出口脹大率定義為在出口外一定位置的流體寬度與出口寬度的比值。對于矩形流道三維流動,定義為流道兩個對稱面處的脹大尺寸與原流道出口尺寸的比值,其中厚度方向的脹大率為B1,寬度方向的脹大率為B2,對于寬厚比大的流道,B1可用以與二維狹縫流道的脹大進行比較。

首先考慮等截面矩形流道的出口脹大情況,三維網(wǎng)格采用自動劃分技術(shù)生成網(wǎng)格。擠出平均速度va=1 mm/s,流道寬度(D)為6 mm,流道厚度用H表示。圖1～圖3分別給出了2種寬厚比(D/H)的流道的有限元網(wǎng)格的脹大圖,以及厚度方向?qū)ΨQ面上的流線圖和速度等值線,其中速度等值線數(shù)值是相對于平均速度的數(shù)值。

圖1 矩形流道擠出形狀Fig.1 Shape of the extrudate flowing out of the rectangular channel

圖2 對稱面上的流線圖Fig.2 Streamline pattern of the symmetrical surface

圖3 對稱面上速度等值線圖Fig.3 Axial velocity contour of the symmetrical surface

對不同寬厚比矩形流道的擠出脹大進行分析,在表2中給出出口12 mm處厚度方向脹大率隨寬厚比的變化情況。從表2可以看出,隨著寬厚比的增加,脹大率B1隨之增加,當(dāng)寬厚比達到6時,厚度的脹大率與二維數(shù)值模擬得到的脹大率一致,與實驗數(shù)據(jù)也很接近[12],而寬度方向的脹大率B2則隨著寬厚比的增加而減少。

表2 不同寬厚比的矩形流道的擠出脹大率Tab.2 Swell ratio of rectangular channel with differentratio of width to thickness

由于黏彈性流體應(yīng)力松弛較為緩慢,熔體在流道的停留時間決定應(yīng)力的松弛程度,所以擠出速度和流道長度對脹大率也有一定的影響。表3給出了不同流道長度的正方形流道擠出脹大率,對于正方形流道,B1與B2相同。可以看出隨著流道長度的增加,脹大率逐漸下降,但當(dāng)流道長度增加到一定程度時,脹大率變化幅度減小。所以流體在流道中的停留時間對脹大率有較大影響。

表3 流道長度對正方形流道擠出脹大率的影響Tab.3 Swell ratio of square channel with different length

對于收斂流道的擠出脹大的分析,考慮2種不同收斂角度的流道,流道總長度相同,為15.6 mm,寬度不變,厚度收斂角分別取作15°和20°。脹大率分別為1.46和1.49,脹大率要大于等截面流道的脹大率1.39,同時由于收斂角為15°的流道的出口流道長度要長些,所以其脹大率小于20°的流道。圖4、圖5給出了15°收斂流道的擠出形狀和流線圖。

圖4 15°收斂流道的擠出形狀Fig.4 Shape of the extrudate flowing out of converging channel with the convergence angle of 15°

圖5 15°收斂流道對稱面上的流線圖Fig.5 Streamline pattern of the symmetrical surface of the 15°converging channel

4 結(jié)論

(1)采用自由面迭代結(jié)合網(wǎng)格自動重新劃分的方法,可以避免流線有限元處理擠出脹大的三維模擬時可能造成的網(wǎng)格畸形,從而較為準(zhǔn)確地預(yù)測擠出脹大的三維模擬的自由面形狀,使擠出脹大的研究更接近實際應(yīng)用;

(2)PE-LD熔體在矩形流道中的擠出脹大過程呈現(xiàn)出厚度脹大率隨寬厚比的增加而增加,脹大率隨流道長度的增加而減小的規(guī)律;

(3)收斂流道的脹大率大于等截面流道的脹大率。

[1] Chang P W,Pattern T W,Finlayson B A.Collocation and Galerkin Finite Element Methods for Viscoelastic Fluid Flow I.Description of Method and Problems with Fixed Geometry[J].Computers&Fluids,1979,7(4):267-283.

[2] Debae F,Legat V,Crochet M J.Practical Evaluation of Four Mixed Finite Element Methods for Viscoelastic Flow[J].Journal of Rheology,1994,38(2):421-443.

[3] Luo X L,Tanner R I.A Streamline Element Scheme for Solving Viscoelastic Flow Problems Part II Integral Constitutive Models[J].Journal of Non-Newtonian Fluid Mechanics,1986,22:61-89.

[4] Luo X L,Tanner R I.Finite Element Simulation of Long and Short Circular Die Extrusion Experiments Using Integral Models[J].International Journal forNumerical Methods in Engineering,1988,25:9-22.

[5] Goubfomme A,Draily B,Crochet M J.Numerical Prediction of Extrudate Swell of a High-density Polyethylene[J].Journal of Non-Newtonian Fluid Mechanics,1992,44:171-195.

[6] Kiriakidis D G,Mitsoulis E.Viscoelastic Simulations of Extrudate Swell for a HDPE Melt Through Slit and Capillary Dies[J].Advanced Polymer Technology,1993,12(2):107-117.

[7] 柳和生,涂志剛,熊洪槐.聚合物熔體在L型異型材擠出口模內(nèi)三維粘彈流動數(shù)值模擬[J].塑性工程學(xué)報,2007,14(2):114-117.

[8] 武停啟,江 波,畢 超,等.聚合物熔體全三維非等溫數(shù)值模擬[J].塑料,2006,35(1):73-79.

[9] 李 勇,江體乾.塑料熔體三維流動的數(shù)值模擬[J].中國塑料,2004,18(9):68-72.

[10] 李 勇,江體乾.聚合物熔體三維擠出脹大的數(shù)值模擬[J].力學(xué)學(xué)報,2002,(6):856-862.

[11] Wagner M H,Demarmels A.A Constitutive Analysis of Extensional Flows ofPolyisobutylene[J].Journal of Rheology,1990,34(6):943-958.

[12] Ahmed R,Liang R F,Mackley M R.The Experimental Observation and Numerical Prediction ofPlanarEntry Flow and Die Swell for Molten Polyethylenes[J].Journal ofNon-NewtonianFluidMechanics,1995,59:129-153.

Numerical Simulation of the Extrudate Swell of Polymer Melt Through Rectangular Channels

LI Yong,ZHEN GJianrong

(School of Mechanical and Power Engineering,East China University of Science and Technology,Shanghai 200237,China)

In this paper,the finite element method of Picard iterative format was applied for analyzing extrudate swell of viscoelastic fluid flows in three-dimension using Wager integral constitutive equations.During each iteration,the updated streamline of boundaries was generated according to boundary condition of velocity.The streamline of Gauss point was acquired according to the velocity field of previous iteration.The stress was treated as pseudo-body and the format of iteration was given.Extrudate swells of different width and length of the rectangular die were analyzed.Results of swell ratios influenced by width and length of die were shown.It showed that the swell ratio increased with increasing ratio of width to thickness,and decreased with increasing length of runner.The results were very accurate and were compared with two dimensional numerical simulation and experiments.It proved that this method was efficient and reliable.

polymer;melt;rectangular channel;extrudate swell;numerical simulation

TQ320.66+3

B

1001-9278(2010)02-0078-05

2009-09-06

聯(lián)系人,yong_li@ecust.edu.cn