Capillary flow of low-density polyethylene.
Capillary rheomeiry is extensively used in both industry and academia lo assess the rheological behavior of polymer melts at high shear rates before testing their processabilily in full industrial scale . When such a How is used and the raw data are collected, a number of important corrections should be applied before the rheo- logical data can be compared with corresponding data from a rotational rheometer [2, 3].
First, capillary flow involves flow through a contraction of a certain angle, where there is a large pressure drop associated with such flow, known as entrance pressure or Bagley correction [1, 4]. This pressure is required in order to calculate the true shear stress, and also, frequently, the apparent extensional rheology of molten polymers, a method well practiced in industry [5-9]. Therefore, it is important to understand the origin of this excess pressure and consequently to be able to predict it.
Many studies have previously attempted to examine the origin of entrance pressure and its prediction for low-density polyethylenes (LDPEs). Feigl and Ottinger  have simulated the axisymmetric contraction flow of an LDPE melt using a multimode Rivlin-Sawyers integral constitutive model, which is a special case of the Kaye-Bernstein, Kearsley, Zapas (K-BKZ) model. However, their numerical results well under-predicted the available experimental pressure data for entry flows. Barakos and Mitsoulis [11, 12] reported similar (hidings with the same K-BKZ model for the Bagley correction in capillary flow of an LDPE melt tested by a working group of the International Union of Pure and Applied Chemistry (IUPAC-LDPE). Beraudo et al. 113] used a multimode Phan-Thien/Tanner (PTT) constitutive relation and found that numerical predictions significantly under-estimated the experimental findings. Guillet et al.  also studied the entrance pressure losses for an LDPE melt both experimentally and numerically using the K-BKZ model. Again significant under-estimation was reported. Similar results were obtained by Hatzikiriakos and Mitsoulis , Mitsoulis et al. , and Mitsoulis and Hatzikiriakos  for different polymer melts with the K-BKZ model. All of these works fitted Theological data for the melts and performed isothermal, viscoelastic simulations with no pressure-dependence of the viscosity. Thus, it had become apparent that pressure underestimation was a problem, especially for the more elastic LDPE melts, which needed some extra attention and more careful work to be resolved.
Along these lines, a recent study by the principal authors  examined three different polyethylenes, namely LDPE, a linear low-density polyethylene (LLDPE), and a high-density polyethylene (HDPE) in flow through tapered dies. The novelty was the pressure-dependence of the viscosity, which was assumed a constant. In isothermal simulations, the problem of predicting the Bagley correction for LDPE was solved satisfactorily at extremely high shear rates (up to 1000 [s.sup.-1]) for the first time  by taking into account the effect of pressure on the viscosity on top of viscoelasticity. Thus, it became apparent that both viscoelasticity and pressure-dependence of the viscosity were significant in obtaining high entry pressures in agreement with experiments. It should be noted that such high entry pressures are important for other materials as well, such as pastes, as have been reported by many in the literature [19-22].
Since the effect of pressure on the viscosity was found to be so important, it was decided to re-examine more carefully experimental data on LDPE. It is well known that using capillary data from dies of various length-to-diameter (L/D) ratios, the dependence of viscosity on pressure can be assessed [2, 23-28]. As a first approximation, the following expression (Barus equation) can be used to determine the parameter, [[beta].sub.p], known as the pressure coefficient of viscosity:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1)
where [eta] is the viscosity at pressure p, and [[eta].sub.[p.sup.0]]is the viscosity at ambient pressure. Using a pressurized Sliding-plate rheometer, Dealy and coworkers [26, 29, 30) have determined the coefficient [[beta].sub.p] for various systems, including LDPE, poly-a-methystyrene-co-acrylonitrile (PaMSAN), and LLDPE. For the LDPE (of main interest in the present work), a value between 1.3 x [10.sup.-8] Pa and 4.9 x [10.sup.-8] Pa has been reported by various authors [24, 28, 31. 32]. This coefficient has also been reported to be a function of the shear rate, with [[beta].sub.p] of LDPE decreasing significantly with shear rate [24, 33]. It is. therefore, a partial objective of the present work to carefully re-examine the pressure-dependence of the viscosity for the LDPE melt at hand.
Another issue is the effect of slip of the material along the die walls. Unlike Newtonian fluids, polymer melts slip over solid surfaces when the wall shear stress exceeds a critical value [34-38]. In particular, slip effects have been reported in the capillary flow of molten polyethylenes [18, 34, 35, 37]. Thus, in a comprehensive study of any melt, such as the one undertaken here, possible slip effects should be studied. This can be done by using at least three capillary dies having different diameters and the same L/D in order to keep the effect of pressure constant. If the flow curve shows a diameter dependency similar to that implied by the explicit relations derived by Mooney , slip is present and the Mooney method can be used to determine the slip velocity as a function of the wall shear stress. It is. therefore, another partial objective of the present work to examine whether the LDPE melt at hand slips at the wall.
Another issue is the effect of the non-isothermality of How. Any analysis for the determination of slip and the dependence of viscosity on pressure are complicated by viscous dissipation (viscous heating) effects. Viscous heating in the capillary extrusion of polymer melts has been the subject of several reviews and studies over the past decades [2, 3, 40-46]. Rosenbaum and Hatzikiriakos  have considered the analysis of the combined effects of viscous heating and wall slip using numerical analysis. They have reported that in cases where viscous heating and slip effects are signilicant, the traditional methods used to estimate the slip velocity (for example Mooneyes technique) often fail to give accurate or even physically meaningful results. They have proposed, a data analysis procedure based on a mathematical model for the non-isothermal capillary flow of polymer melts coupled with heat transfer to analyze the capillary data instead. Pressure-dependence of viscosity and compressibility effects have been neglected in the analysis by Rosenbaum and Hatzikiriakos . More recently, Laun [2, 3] considered the superimposed effect of dissipative heating and a pressure-dependent viscosity in circular die flows at high shear-rates. The problem was treated theoretically, and pertinent relations for capillary rheometry were derived. The applicability of the data treatment was demonstrated for the die flow of LDPE. However, compressibility and slip effects were not considered.
It is, therefore, the main objective of this work to study in a comprehensive way the capillary flow of an LDPE melt both experimentally and numerically by considering all possible capillary effects, combined and separately. Namely, the effects of viscoelasticity, compressibility, pressure-dependence of viscosity, and viscous healing (non-isothermal flow) will be considered in order to assess the significance of each effect. It will be shown that viscoelasticity and the pressure-dependence of viscosity are primarily the reasons for high values of the entry pressure losses, especially at elevated shear rates, for LDPE melts,
A typical LDPE melt was used in this work in order to address the effects of pressure, temperature, and compressibility on its capillary flow. The viscosity of LDPE melts is known to be more sensitive to temperature and pressure compared to that of other polyethylenes, such as HDPE and/or LLDPE, and it is the main reason for being selected in this study. It is also known that branched polymers exhibit extensional strain-hardening effects, and therefore increased entrance pressure values are obtained compared to those obtained for their linear counterparts [15-17]. This is expected as the entry pressure has been theoretically associated with extensional viscosity [5-9]. The properties of the particular LDPE used in this study are summarized in Table 1.
TABLE 1. Properties of the LDPE resin used in this study. Sample Zero-shear-rate Resin Melt index Density viscosity ID producer (g/10 min) (g/[cm.sub.3]) (Pa s) (190[degrees]C) (25[degrees]C) 160[degrees]C LDPE DOW 0.47 0.919 617,230 Chemicals
Several rheological experiments were carried out in order to characterize the polymer in the molten state. First, an Anton Paar MCR-501 device was used in the parallel-plate geometry to determine its linear viscoelastic moduli over a wide range of temperatures from 130 to 230 [degree]C. Master curves were obtained using the time-temperature superposition (TTS) principle, and the results are presented at the reference temperature of 160 [degree]C. These values are reported in Table 1 and are found to compare well with those obtained from linear viscoelastic measurements.
The polymer was also rheologically characterized in simple extension using the SER-2 Universal Testing Platform from Xpansion Instruments [48, 49]. Details of sample preparation and operating principles of this rheometer to obtain reliable results are given in Delgadillo et al. [50, 51]. Uniaxial extension tests were performed at the reference temperature of 160 [degree]C.
Capillary Flow Testing
An Instron capillary rheometer (constant piston speed) was used to determine the entrance pressure (known also as Bagley method)  and the viscosity as a function of the wall shear stress, crw, and apparent shear rate, [[gamma].sub.A] = 32Q/[pi][D.sup.3] at 160 [degree]C, where Q is the volumetric How rate and D is the capillary diameter. A series of orifice dies (D = 0.0508 cm and L/D = 0) were also used to determine directly the entrance pressure (Bagley correction) as a function of the apparent shear rate for various contraction angles (15, 30, 45, 60, and 90[degree]).
In addition, several dies having the same diameter (D = 0.0787 cm) and various LID ratios (L/D = 2, 5, 16, 33) have been used to determine the effect of pressure on viscosity. These data were also used to determine the Bagley correction and compared with the one determined from the orifice die. The agreement found was satisfactory. In addition, capillary dies having the same L/D ratio of 20 and different diameters (D = 0.0508, 0.0762, and 0.152 cm) were used to examine slip effects for the case of LDPE by means of the Mooney technique .
GOVERNING EQUATIONS AND RHEOLOGICAL MODELING
We consider the conservation equations of mass, momentum, and energy for weakly compressible fluids under non-isothermal, creeping, steady flow conditions. These are written as [52, 53]:
[bar.u] [nabla][rho] + [rho]([nabla][bar.u]) = 0, (2)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)
Where [rho] is the density, [bar.u] is the velocity vector, p is the pressure, [??] is the extra stress tensor, T is the temperature, Cp is the heat capacity, and k is the thermal conductivity. For a compressible fluid, pressure and density are connected as a first approximation through a simple linear thermodynamic equation of state :
[rho] = [[rho].sub.0](1 + [[beta].sub.c]p), (5)
where [[beta].sub.c] is the isothermal compressibility with the density to be [[rho].sub.0] at a reference pressure [p.sub.0] (= 0).
The viscous stresses are given for inelastic non-Newtonian compressible fluids by the relation [52, 53]:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)
where [eta](|period[gamma]|) is the apparent non-Newtonian viscosity, which is a function of the magnitude |period[gamma]| of the rate-of-strain tensor [??] = [nabla][bar.u] + [nabla][[bar.u].sup.T], which is given by:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)
where [II.sub.[gamma]] is the second invariant of [??]
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (8)
The tensor [??] in Eq. 6 is the unit tensor.
To evaluate the role of viscoelasticity in the prediction of Bagley correction, it is instructive to consider first purely viscous models in the simulations. Namely, the Carreau-Yasuda model was used to fit the shear viscosity data of the LDPE melt. The Carreau-Yasuda model is written as :
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (9)
where [[eta].sub.0, c] is the Carreau zero-shear-rate viscosity, [lambda] is a time constant, [n.sub.C] is the Carreau power-law index, and a is the Yasuda exponent (= 2 for the simple Carreau model). The fitted viscosity of the LDPE melt by Eq. 9 is plotted in Fig. 1, while the parameters of the model are listed in Table 2. We observe that the LDPH melt is very shear-thinning for shear rates above 1 [s.sup.-1] due to the presence of long-chain branching . The Carreau-Yasuda model fits the data well over the range of experiment results.
TABLE 2. Parameters for the LDPE resin obeying the Carreau-Yasuda model (Hq. 9) and the power-law model (Eq. 10) at 160 [degrees]C. Parameter Value [[eta].sub.0] 617,230 Pa s [lambda] 17.727 s [n.sub.c] 0.221 [alpha] 0.215 K 13,780 Pa [s.sup.n] n 0.4
For easy checks and simple analytical formulas, the high shear-rate range of the viscosity data can be lilted to the power-law model for the viscosity [l]:
[eta] = K[|period[gamma]|.sup.n - 1], (10)
where K is the consistency index and n is the power-law index. These values are also given in Table 2.
For the capillary flow simulations, the effect of pressure on viscosity should be taken into account as this becomes evident below. This effect is quite significant for LDPE compared to LLDPE and less significant for HDPB [31, 32]. This effect can be taken into account by multiplying the constitutive relation with a pressure-shift factor, [a.sub.p], defined by the Barus equation, that is [27, 32]:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (11)
where [[beta].sub.p] is the pressure coefficient and p is the absolute pressure, as discussed above.
Viscoelasticily is included in the present work via an appropriate rheological model for the stresses. This is a K-BKZ equation proposed by Papanastasiou et al.  and modified by Luo and Tanner . This is written as:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)
Where t is the current time, [[lambda].sub.k] and [a.sub.k] are the relaxation times and relaxation modulus coefficients, N is the number of relaxation modes, [alpha] and [beta] are material constants, and [I.sub.c], [I.sub.[c.sup.-1]] are the first invariants of the Cauchy-Green tensor [C.sub.t], and its inverse [C.sub.t.sup.-1] the Finger strain tensor. The material constant 0 is given by
[N.sub.2]/ [N.sub.1] = [theta]/(1 - [theta]) (13)
where [N.sub.1], and [N.sub.2] are the first and second normal stress differences, respectively. It is noted that [theta] is not zero for polymer melts, which possess a non-zero second normal stress difference. lis usual range is between -0.1 and -0.2 in accordance with experimental findings .
As discussed above, experiments were performed in the parallel-plate and extensional rheometers for the LDPE melt to rheologically characterize it. Figure 2 plots the master dynamic moduli C' and G" for LDPE at the reference temperature of 160 [degrees]C. The model predictions obtained by fitting the experimental data to Eq. 12 with a spectrum of relaxation times, [[lambda].sub.k] and coefficients, [a.sub.k], determined by a non-linear regression package , are also plotted. The parameters found from the tilting procedure are listed in Table 3. The relaxation spectrum is used to find the average relaxation time, [bar.[lambda]], and zero-shear-rate viscosity, [[eta].sub.0] according to the formulas:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (14)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)
TABLE 3. Relaxation spectra and material constants lor the LDPE resin obeying the K-BKZ model (Eq, 12) at 160 [degrees]C. [alpha] = 7.336, [theta] = 0 k [[lambda].sub.k] (s) [a.sub.k](Pa) [[beta].sub.k] 1 1.28 x [10.sup.-4] 4.50 x [10.sup.5] 1 2 2.51 x [10.sup.-3] 98,795 1 3 2.06 x [10.sup.-2] 48,899 0.18 4 1.62 x [10.sup.-1] 22.089 0.45 5 1.224 8842 0.049 6 6.733 3397 0.026 7 43 948.5 0.024 8 248 287 0.014
The values of these parameters are [bar.[lambda]] = 130 s, [[eta].sub.0] = 150,607 Pa s, indicating a highly elastic melt with a long average relaxation time.
Figure 3 plots a number of calculated and experimental material functions for the LDPE melt at the reference temperature of 160 [degree]C. Namely, data for the shear viscosity, [eta]s, the elongational viscosity, [eta]E, and the first normal stress difference, [N.sub.1] are plotted as functions of corresponding rates (shear or extensional). The parameter [beta] that controls the calculated elongational viscosity was fitted by using the extensional behavior of the melt. Figure 4 shows the extensional behavior of the LDPE at several Hencky strain rates at 160 [degree]C and the model predictions of Eq. 12 using multiple [beta]-values listed in Table 3. It can be seen that the overall rheological representation of all material functions is excellent.
For the non-isothermal calculations, it is necessary to derive a non-isothermal constitutive equation from the isothermal one. This is done by applying the time-temperature shifting concept as explained by Luo and Tanner . This concept is based on the relative difference between the observer's time scale and the material's internal time scale. Some intuition can be obtained considering that as temperature raises so does the amount of molecular motion occurring in one unit of observer's lime. Using [zeta] for the time measured by the material's own internal clock, the following relation holds between [zeta] and the observer's time t:
d[zeta] = dt/[a.sub.T](T) (16)
where the denominator is the time-shifting factor . The above relationship expresses the Morland-Lee hypothesis [52, 58] in a differential form and can be used to obtain the following integral relation between the particle's elapsed time and the observing period:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (17)
Thus as a fluid particle is tracked along a streamline segment [DELTA][l.sub.i], the particle's time corresponding to the residence time [DELTA]t' is given by
[DELTA][zeta] = [DELTA][l.sub.i]/[V.sub.i][a.sub.T]([T.sub.i]) (18)
where [V.sub.i], is the particle velocity.
The above Eq. 16 has been used by Luo and Tanner  to obtain the non-isothermal form of the K-BKZ equation (Eq. 12). This was done by replacing the observer's time t in Eq. 12 with the particle's lime [zeta] given by Eq. 17. The resulting version of the constitutive equation becomes
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (19)
For LDPE melts, the following temperature-shifting function has been found to be adequate [52, 59]:
[a.sub.T](T) = [eta]/[[eta].sub.0] = exp[(E/[R.sub.g])(1/T - 1/[T.sub.0])] (20)
In the above, [[eta].sub.0] is a reference viscosity at [T.sub.0], E is the activation energy constant, [R.sub.g] is the ideal gas constant, and [T.sub.0] is a reference temperature (in K). The activation energy constant can be determined from zero-shear-rate viscosity data ,
Another method can be applied to derive the non-isothermal constitutive equation from the isothermal one. This method is again based on the time-temperature superposition temperature and simply consists of shifting the relaxation times [[lambda].sub.k] from the temperature history within the material's time scale . The equation used to shift the relaxation times in the material's history is given by [60, 61]:
[[lambda].sub.k](T'(t')) = [[lambda].sub.k]([T.sub.0])[a.sub.T](T'(t')). (21)
Both methods have been applied before for the process of film blowing  with virtually the same results. Test runs here confirmed this finding and the current results were derived with the [zeta]-approach.
The viscoelaslic stresses calculated by the above constitutive equation (Eq. 19) enter in the energy equation (Eq. 4) as a contribution to the viscous dissipation term. This assumes that all the elastic stored energy in the material is converted into heat and serves to raise the temperature of the melt. This treatment has been the standard practice in non-isothermal simulations so far (see, e.g., Barakos and Mitsoulis ) and it is followed in the present work as well. Peters and Baaijens  have devised a different approach to handle the splitting between the viscous contribution to energy dissipation and elastic stored energy. However, such a development is beyond the scope of this study.
The thermal properties of the melt have been gathered from various sources. The thermal conductivity was given by Laun  to be 0.17 W/m K, and the heat capacity was equal to 2.25 x [10.sup.3] J/(kg K).
The density is given by the equation :
[rho] (T) = 0.8772 - 0.000597, for T > 150 [degrees]C. (22)
where [rho] is in g/[cm.sup.3] and T is in [degrees]C. Its value at [T.sub.0] = 160 [degree]C is [rho] = 0.7828 g/[cm.sup.3].
The activation energy for this particular melt was found by applying the time-temperature superposition to be equal to E = 64,100 J/mol [50, 51]. Others have reported similar values . All these properties and their provenance are gathered together in Table 4. In this table, we have also added the value for the coefficient of compressibility, [[beta].sub.c] [64) and the value of a constant pressure-coefficient of viscosity, [[beta].sub.p] .
TABLE 4. Values of the various parameters for the LDPE resin at 160 [degree]C. Parameter Value Reference [[beta.sub.c]] 0.00095 [MPa.sup.-1]  [[beta].sub.p] 0.015 [MPa.sup.-1]  M 0.65 [GPa.sup.[n.sub.p] - 1] This work (0.11 [GPa.sup.[n.sub.p] - 1]) [n.sub.p] 0.75 This work P 0.7828 g/[cm.sup.3]  [C.sub.p] 2.25 J/(g K)  K 0.0017 J/(s cm K)  E 64,100 J/mol [50, 51] [R.sub.g] 8.3143 J/(mol K)  [T.sub.0] 160 [degrees]C(433 K) [50, 51] Na = ([bar.[eta]]E[U.sup.2])/(k[R.sub.g][T.sub.0.sup.2]) (24)
The various thermal and flow parameters are combined to give appropriate dimensionless numbers [40. 44]. The relevant ones here are the Peclet number, Pe, and the Nahme-Griffith number, Na. These are defined as
Pe = [rho][C.sub.p]UR/K (23)
where [bar.[eta]] =f(U/R) is a nominal viscosity given by the Carreau-Yasuda model (Eq. 9) at a nominal shear rate of U/R, Ra is the ideal gas constant, and U (= [period[gamma].sub.A] R/4) is the average velocity in the capillary die. The Pe number represents the ratio of heat convection to conduction, and the Na number represents the ratio of viscous dissipation to conduction and indicates the extent of coupling between the momentum and energy equations. A thorough discussion of these effects in non-isothermal polymer melt flow is given by Winter .
With the above properties and a die radius R = 0.0254 cm, the dimensionless thermal numbers are in the ranges: 3.3 < Pe < 658, 0.002 < Na < 3, showing a relatively strong convection (Pe [much greater than] 1), and a moderate coupling between momentum and energy equations (Na ~ l). A value of Na > 1 indicates temperature non-uniformities generated by viscous dissipation, and a strong coupling between momentum and energy equations. For a die radius R = 0.04 cm, the corresponding ranges are: 5.2 < Pe < 1036, 0.003 < Na < 4.17. More details are given in Table 5.
TABLE 5. Range of the dimensionless parameters in the flow of LDPE resin at 160 [degree]C (die radius R = 0.04 cm). Peclet Nahme Apparent Shear Rate, Number, Number, [[gamma].sub.A] ([s.sup.-1]) Pe Na 5 5.2 0.003 11 11.4 0.009 64 66.3 0.10 390 404.1 1.20 1000 1036.1 4.17
The fully developed values for temperature occur for L/D = [varies] and give a maximum temperature rise [DELTA][T.sub.max] at the centreline of the tube according to the formula for a power-law fluid :
[DELTA][T.sub.max] = K/k[((nR)/(3n + 1)).sup.2][[((n + 1/n)[V.sub.max]/R)].sup.n + 1] (25)
According to Eq. 25 and for R = 0.0254 cm, we get a [DELTA][T.sub.max] = 29.2 [degree]C. For R = 0.04 cm, we get a [DELTA][T.sub.max] = 38.3 [degree]C. Thus, a significant temperature rise is expected for very long dies, but for shorter dies as the ones used here, the non-isothermal simulations are necessary to find out the extent to which temperatures rise.
Similarly with the time-temperature superposition principle where the stresses are calculated at a different temperature using the shift factor [a.sub.T], the time-pressure superposition principle can be used to account for the pressure effect on the stresses. In both cases of viscous or viscoelastic models, the new stresses are calculated using the pressure-shift factor [a.sub.p]. For viscous models, Eq. 11 is used to modify the viscosity. For viscoelastic models, such as the K-BKZ model (Eq. 12), the pressure-shift factor modifies the relaxation moduli, [a.sub.k], according to:
[a.sub.k](p(t')) = [a.sub.k]([p.sub.0])[a.sub.p](p(t')) (26)
This is equivalent to multiplying the stresses by [a.sub.p], according to Eq. 12. It should be noted that [a.sub.p] is an exponential function of [[beta].sub.p], which itself may depend on pressure p. Since in a flow field negative pressures may appear, especially around singularities as is the case in contraction flows, special care must be taken numerically to handle these functions for negative numbers. Failure to do so leads to nonsensical results and/or to divergence. It should be noted that in our recent work , the pressure dependence of the viscosity was not handled in this way, but simply the overall pressure results were multiplied by the [a.sub.p]-factor. Thus, this intrinsic handling of the pressure, which assumes different local values, is a major contribution of the present work.
METHOD OF SOLUTION
The solution to the above conservation and constitutive equations is carried out with two codes, one for viscous flows (u-v-p-T-h formulation)  and one for viscoelastic flows [62, 67]. The boundary conditions (BC) for the problem at hand are well known and can be found in our earlier publication . Briefly, we assume no-slip and a constant temperature To at the solid walls; at entry, a fully-developed velocity profile is imposed, corresponding to the flow rate at hand, and a constant temperature [T.sub.0] is assumed; at the outlet, zero surface traction and zero heat flux are assumed; at the centerline, symmetry is assumed. In one case, adiabatic BC are assumed on the solid walls to study the differences with the constant-wall temperature BC.
The novelty in the present work is the addition of the Barus equation (Eq. 11) for the pressure-shift factor [a.sub.p] with either a constant [[beta].sub.p] or an experimentally found function [[beta].sub.p] = [[beta].sub.p] (p), and the calculation of the viscoelastic stresses based on Eq. 26. Furthermore, the non-isothermal viscous and viscoclastic calculations were checked again against analytical solutions for the power-law fluid with viscous heating having either isothermal or adiabatic walls, and against previous numerical solutions for extrudate swell . It is surprising that since the early work by Luo and Tanner , and Barakos and Mitsoulis [621, no other non-isothermal viscoclastic works seem to have appeared in the literature for integral constitutive equations, except for simplified 1-D models for some polymer processes dominated by extension [60, 61]. This despite the fact that polymer processing is almost always occurring under non-isothermal conditions.
The viscous simulations are extremely fast and are used as a first step to study the whole range of parameter values and die designs. The viscoelastic simulations admittedly are much harder to do and they need good initial flow fields to get solutions at such elevated apparent shear rates as 1000 [s.sup.-1]. In our recent work , we explained how it was possible for the first time to do viscoelastic computations up to the highest apparent shear rates with good results. However, those calculations were carried out under isothermal conditions and with no implicit internal calculations for pressure-dependence. The final result for the pressures was multiplied by [[alpha].sub.p] with a constant value of [[beta].sub.p]. In the current simulations, these limitations are lifted, and full viscoelastic non-isothermal simulations are carried out with a pressure dependence of the relaxation moduli for the first time.
The iteration method is direct substitution (Picard scheme). The solution strategy starts at a given apparent shear rate from the viscous non-isothermal solution with pressure dependence, and then using this as an initial solution it continues for the non-isothermal viscoelastic solution without pressure dependence. Once a new solution is achieved, more iterations are performed for the pressure-dependence of the relaxation moduli. In our earlier efforts, only a constant [[beta].sub.p] value could be handled, but in the final efforts an experimentally found [[beta].sub.p]-function gave well behaved results for all cases, provided care was taken to handle numerically correctly the negative pressures in the field.
The results for L/D = 0 were obtained as in Ref. 18 by subtracting from the longest dies results obtained for the same dies without the reservoir. Thus, it was found again necessary to use a long die to fully unravel the high viscoelastic stresses of the LDPE melt.
Entrance (End) Pressure
Figure 5 plots the entrance pressure (or end pressure due to L/D = 0) of LDPE at 160 [degree]C as a function of the die length L/D for an extended range of values of the apparent shear rate from 5 [s.sup.-1] to 1000 [s.sup.-1]. This plot is also known as the Bagley plot. The data from an orifice die are also plotted. The data show a curvature with concavity pointing upwards consistent with Eq. I (effect of pressure on viscosity). These results are for a contraction angle 2 [alpha] = 90[degree]. Our recent work  has shown that the entry pressure, although significant, is nearly independent of contraction angle for angles 2 [alpha] [greater than or equal to] 45[degree].
The Flow Curves
Figure 6 plots the How curves of LDPE at 160 [degree]C obtained with capillary dies having lengths of L/D = 2, 5, 16, and 33. The Bagley and Rabinowitsch corrections have been applied to the data. Also plotted are the linear viscoelastic data at 160 [degree]C, which agree well with the data at L/D = 2 and 5. The data clearly shows that there is an effect of pressure on viscosity. The flow curves have been corrected for these effects according to Barus Equation (Eq. I) and the procedure discussed by others [23, 27, 28]. Based on the superposition procedure, the data for the higher L/D's are shifted horizontally to superpose with the LVE data. In this way the shifting parameter [[beta].sub.p.sup.super] (which is a shear-rate-dependent pressure coefficient) can be obtained. According to Couch and Binding  this parameter could be related to Barus pressure coefficient [[beta].sub.p] as:
The calculated coefficients [[beta].sub.p] depend on the shear rate and the L/D (level of pressure). Duvdevani and Klein  and Liang  have also reported significant shear-rate dependence of the coefficient [[beta].sub.p]. Although the pressure dependence of [[beta].sub.p] has been reported in some studies [68, 69], it has been accepted for practical reasons that this parameter is independent of pressure [70-73]. In our case, [[beta].sub.p] shows a power-law dependency on pressure (see Fig. 7), according to the following equation:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (28)
where m = 0.62 [GPa.sup.[n.sub.p]-1] and [n.sub.p] = 0.75, with the pressure p given in GPa. The pressure-corrected data for the LDPE are plotted in Fig. 8. The data superposes well and the capillary data agrees well with the LVE data. The values of [[beta].sub.p] plotted in Fig. 8 agree well with the reported values in the literature [23, 27, 28, 31].
The diameter dependence of the flow curve has also been studied in order to assess the slip behavior of LDPE. The results are shown in Fig. 9. No diameter dependence is evident as the flow curves for three significantly different diameters practically superpose. However, if viscous heating effects are significant, this may mask the diameter dependence as discussed by Rosenbaum and Hatzikiriakos .
It is instructive to perform first calculations with a purely viscous model, so that the effect of viscoelasticity will become evident later. The numerical simulations have been undertaken using the viscous Carreau-Yasuda model (Eq. 9). This constitutive relation is solved together with the conservation equations of mass and momentum either for an incompressible or compressible fluid under isothermal or non-isothermal conditions (conservation of energy equation) with or without the effect of pressure-dependence of the viscosity.
For the finite element mesh arrangement, we have used our long experience with viscous and viscoelastic flows and chosen a grid that progressively adds more elements as one moves towards the singularity at the entrance to the die, while the elements become bigger as one moves away from the singularity. A typical finite element grid is shown in Fig. 10 for L/R = 40. The domain represents an 18.75:1 tapered circular contraction with an entrance angle 2 [alpha] = 90[degree]. The grid consists of 1700 elements, 7161 nodes, and 22,261 unknown degrees of freedom (d.o.f.), while a four-times denser grid is also used, having been created by subdivision of each element into four sub-elements for checking purposes of grid-independent results. This checking consists of reporting the overall pressures in the system from the two meshes and making sure that the differences are less than 1% between the two results. Having fixed the Carreau-Yasuda model parameters and the problem geometry, the only parameter left to vary was the apparent shear rate in the die ([[gamma].sub.A] = 4Q/[pi][R.sup.3]). Simulations were performed for the whole range of experimental apparent shear rates, namely from 5 [s.sup.-1] to as high as 1000 [s.sup.-1].
First, runs were carried out to study each effect separately, namely the effect of compressibility alone, the effect of a pressure-dependent viscosity alone, and the effect of a temperature-dependent viscosity alone, and these are compared with the base case of no effects at all ([[beta].sub.c] = [[beta].sub.p] = [[alpha].sub.T] = 0). These are shown in Figs. 11-13, respectively. We note in Fig. 11 that compressibility plays a negligible role, as the results are virtually the same with or without compressibility. Admittedly, this is expected because the value for LDPE of [[beta].sub.c] = 0.00095 [MPa.sup.-1] is small . On the other hand, a pressure-dependence of the viscosity, with a constant value of [[beta].sub.p]]= 0.015 [MPa.sup.-1] , plays a significant role (see Fig. 12), augmenting the pressure drop substantially, and this effect is stronger as the apparent shear rate increases. Finally, a temperature-dependent viscosity serves to lower somewhat the pressure in the system (see Fig. 13), but again this effect is negligible for small to intermediate apparent shear rates, and becomes more evident at the high end of shear rate values. The non-isothermal runs have been made with isothermal walls at [T.sub.0] = 160 [degrees]C. At the highest apparent shear rate of 1000 [s.sup.-1] we have also assumed adiabatic walls, and this brings down the pressure appreciably (from 19 MPa to 15 MPa). In reality, the walls are neither isothermal nor adiabatic but somewhere in between; however they are much closer to the isothermal case [40, 44]. We note that the effect of a temperature-dependent viscosity is opposite to the effect of a pressure-dependent viscosity, and these two effects are in competition. The overall effect of the three contributions is shown in Fig. 14, where it can be seen that the pressure is higher than the base case for all apparent shear rates, driven mainly by the pressure-dependent viscosity effect.
When the die is longer (L/D = 40), the overall effects are more pronounced as evidenced in Fig. 15. Again, the main effect is the pressure-dependence of viscosity, which--given a longer die--has more length to show its influence. The pressure drop in the die is not linear but concave, although this is more evident at the highest apparent shear rate of 1000 [s.sup.-1].
The effect of temperature is not very strong, as evidenced in Fig. 16. For L/D = 20, the centreline temperature distribution is shown for different apparent shear rates for isothermal walls and for the two cases of [[beta].sub.p] = 0 and [[beta].sub.p] = 0.015 [MPa.sup.-1]. The maximum temperature rise at the centreline does not exceed 3 [degrees]C, although near the wall is 1-2 degrees higher due to a higher viscous dissipation contribution there.
Collecting all the pressure results together for the live dies with L/D = 0.2, 2, 5, 16, 33, and 5 apparent shear rates [[gamma].sub.A] gives in Fig. 17 the well known Bagley plot. The experimental results are shown as symbols while the numerical results here and in the subsequent graphs are shown as lines. These simulation results have been obtained with a variable [[beta].sub.p] given by Eq. 28. The results are below the experimental data, and this is mainly due to the inability of a viscous model to capture the end correction (for L/D = 0). A single value of [[beta].sub.p] = 0.015 [MPa.sup.-1] gives somewhat higher values, especially for the highest LID and [[gamma].sub.A] values, but it still cannot account for the L/D = 0 end correction. As it was found out in Ret. 18, a purely viscous model for LDPH gives end corrections one order of magnitude lower than the experimental ones. Ii should be noted that underestimates, such as these, have been found in all previous viscous works, even though the effect of pressure was not taken into account in all of them. So the effect of pressure on the viscosity alone is not enough to correctly predict the pressure drop. It is then at this point that we turn our attention to the viscoelastic results.
It is again instructive to compare individual distributions of pressure and temperature between the viscous and the viscoelastic models. This is done in Fig. 18, where we show the pressure results from the Carreau-Yasuda and the K-BKZ models for the case of L/D = 33 with D = 0.0787 cm. The viscoelastic pressures are markedly higher than the viscous ones, and this becomes more apparent for elevated shear rates. The small kinks at z/R = 0 correspond to the singularity at the entry of the die, which is more severe for viscoelastic fluids than for viscous ones. The straight lines in the die region are a consequence of the counter-effects of temperature and pressure, as is well known in rheometry of polymer melts.
The corresponding results for the centreline temperature distributions are shown in Fig. 19 for the same conditions. Again the temperatures from the viscoelastic runs are higher and they reach in the highest conditions a centreline value of 176 [degree]C, an increase of 16 [degree]C (cf. the ultimate increase of 38.3 [degree]C for L/D = [varies]). The higher temperatures are a consequence of the higher viscoelastic stresses, which give rise to more viscous dissipation, especially near the walls.
Collecting again all the pressure results together for the five dies with L/D = 0, 2, 5, 16, 33 and different apparent shear rates [[gamma].sub.A] gives in Fig. 20 the corresponding viscoelastic Bagley plot. We now see an excellent agreement between experiments and simulations. To our knowledge this is the first time that such a good agreement is obtained for LDPE. It is a combination of temperature and pressure dependence that is responsible for that. The pressure dependence with the variable [[beta].sub.p]-value given by Eq. 28 was essential in capturing this behavior. A constant [[beta].sub.p]-value not only gave wrong results but also led to divergence in the most severe cases; it seems to lead to inadmissibly high pressures because of the exponential function in Eq. l.
A commercial LDPE has been studied in entry flows through tapered dies with different L/D ratios with the purpose of predicting the Bagley correction. The experiments have shown a distinct pressure-dependence of viscosity with a pressure coefficient to be a power-law function of pressure. Full rheological characterization was carried out both with a viscous (Carreau-Yasuda) and a viscoelastic (K-BKZ) model. All necessary material properties data were collected for the simulations.
The viscous model was found to underestimate the extrusion pressures. The viscoelastic model showed a very good agreement with the experimental results, which appears to be the first in the literature for the elastic LDPE melt. The simulations showed that: (a) compressibility is not important in these steady flows; (b) viscous dissipation is important, especially for the more severe conditions (high L/D and apparent shear rates); and (c) the pressure-dependence of viscosity is very important and its correct function has to be found experimentally. This is the first time that all these effects are taken into account in a viscoelastic simulation. As in our recent work , the present viscoelastic results are unique in reaching levels of apparent shear rates as high as 1000 [s.sup.-1] and showed the importance of including long lengths to fully relax the viscoelastic stresses for highly viscoelastic polymer melts, such as this LDPE, in order to capture the end correction.
(1.) J.M. Dealy and K.F. Wissbrun, Melt Rheology and its Role in Plastics Processing--Theory and Applications, Van Nostrand Reinhold, New York (1990).
(2.) H.M. Laun, Rheol. Acta, 42, 295 (2003).
(3.) H.M. Laun, Rheol. Acta, 43, 509 (2004).
(4.) E.B. Bagley, J. Appl. Phys., 28, 193 (1957).
(5.) F.N. Cogswell, Trans. Soc. Rheol., 16, 383 (1972).
(6.) F.N. Cogswell, Polymer Melt Rheology--A Guide to Industrial Practice, Wiley, New York (1981).
(7.) D.M. Binding, J. Non-Newtonian Fluid Mech., 27, 173 (1988).
(8.) D.M. Binding, J. Non-Newtonian Fluid Mech., 41, 27 (1991).
(9.) M. Padmanabhan and C.W. Macosco, Rheol. Acta, 36, 144 (1997).
(10.) K. Feigl and H.C. Ottinger, J. Rheol., 38, 847 (1994).
(11.) G. Barakos and E. Mitsoulis, J. Rheol., 39, 193 (1995).
(12.) G. Barakos and E. Mitsoulis, J. Non-Newtonian Fluid Mech., 58,315 (1995).
(13.) C. Beraudo, T. Coupez, A. Fortin, Y. Demay, B. Vergnes, and J.-F. Agassant, in "Proceedings of XIIth International Congress on Rheology," A. Ait-Kadi, J.M. Dealy, D.F. James, and M.C. Williams, Eds., Laval Universily, Quebec City, Canada, 417 (1996).
(14.) J. Guillet, P. Revenue, Y. Bereaux, and J.-R. Clermont, Rheol. Acta, 35, 494(1996).
(15.) S.G. Hatzikiriakos and E. Mitsoulis, Rheol. Acta, 35, 545 (1996).
(16.) E. Mitsoulis, S.G. Hatzikiriakos, K. Christodoulou, and D. Vlassopoulos, Rheol. Acta, 37, 438 (1998).
(17.) E. Mitsoulis and S.G. Hatzikiriakos, Rheol. Acta, 42, 309 (2003).
(18.) M. Ansari, A. Alabbas, E. Mitsoulis, and S.G. Hatzikiriakos, Int. Polym. Proc., 25, 287 (2010).
(19.) D.J. Horrobin and R.M. Nedderman, Chem. Eng. Sic., 53, 3215 (1998).
(20.) A.G. Gibson," Converging Dies", in Rheological Measurements, A.A. Collyer and D.W. Gregg, Eds., Chapman and Hall, London, 455 (1998).
(21.) A.B. Ariawan, S. Ebnesajjad, and S.G. Hatzikiriakos, Can. J. Chem. Eng., 80, 1153 (2002).
(22.) I. Ochoa and S.G. Hatzikiriakos, Powder Technol, 153, 108 (2005).
(23.) M.A. Couch and D.M. Binding, Polymer, 41, 6323 (2000).
(24.) J.-Z. Liang, Polymer, 42, 3709 (2001).
(25.) R. Cardinaels, P. Van Puyvelde, and P. Moldenaers, Rheol. Acta, 46, 495 (2007).
(26.) H.E. Park, S.T. Lim, H.M. Laun, and J.M. Dealy, Macromo-lecules,41, 1023 (2008).
(27.) Y. Son, J. Polym. Res., 16, 667 (2009).
(28.) J. Aho and S. Syrjala, J. App. Polym. Sci., 117, 1076 (2010).
(29.) F. Koran and J.M. Dealy, J. Rheol., 43, 1279 (1999).
(30.) H.E. Park and J.M. Dealy, Macromolecules, 36, 5438 (2006).
(31.) T. Sedlacek, M. Zatloukal, P. Filip, A. Boltizar, and P. Saha, Polym. Eng. Sci., 44, 1328 (2004).
(32.) E.S. Carrcras, N. El Kissi, J.-M. Piau, F. Toussaint, and S. Nigen, Rheol Acta. 45, 209 (2006).
(33.) I.J. Duvdevani and I. Klein, Soc. Plast. Eng. Trans. J., 23, 41 (1967).
(34.) A.V. Ramamurthy, J. Rheol., 30, 337 (1986).
(35.) D. Kalika and M.M. Denn, J. Rheol., 31, 815 (1987).
(36.) S.G. Hatzikiriakos and J.M. Dealy, J. Rheol., 35, 497 (1991).
(37.) S.G. Hatzikiriakos and J.M. Dealy, J. Rheol., 36, 703 (1992).
(38.) L.A. Archer, "Wall Slip: Measurement and Modelling Issues," in Polymer Processing Instabilities, S.G. Hatzikiriakos and K.B. Migler, Eds., Marcel Dekker, New York,73 (2005).
(39.) M. Mooney, J. Rheol., 2, 210 (1931).
(40.) H.H. Winter, AM Heat Transfer, 13, 205 (1977).
(41.) R.M. Ybarra and R.E. Eckert, AlChE J., 26, 751 (1980).
(42.) S.M. Dinh and R.C. Armstrong, AlChE J., 28, 294 (1982).
(43.) J.F. Milthorpe and R.I. Tanner, Int. J. Bum, Methods. Eng., 24,263(1987).
(44.) E. Mitsoulis, R. Wagner, and F.L. Heng, Polym. Eng, Sci., 28, 291 (1988).
(45.) R.C. Warren, Viscous Hearing. Rheological Measurement, Elsevier, London, 119 (1988).
(46.) Y.S. Ko and A.S. Lodge, Rheol. Acta, 30, 357 (1991).
(47.) E.E. Rosenbaum and S.G. Hatzikiriakos, AlChE J., 43, 598 (1997).
(48.) M.L. Sentmanat and Goodyear Co., U.S. Patent 6,578,413 (2003).
(49.) M.L. Sentmanat, Rheol. Acta, 43, 657 (2004).
(50.) O. Delgadillo-Velazquez, S.G. Hatzikiriakos, and M. Sentmanat, Rheol. Acta, 47, 19 (2008).
(51.) O. Delgadillo-Velazquez, S.G. Hatzikiriakos, and M. Sentmanat, J. Polym. Sci. Part B: Polym Phys., 46, 1669 (2008).
(52.) R.I. Tanner, Engineering Rheology, 2nd ed., Oxford University Press, Oxford (2000).
(53.) E. Mitsoulis and S.G. Hatzikiriakos, J. Non-Newtonian Fluid Mech., 157, 26 (2009).
(54.) J.M. Dealy and R.G. Larson, Structure and Rheology of Molten Polymers-from Structure to Flow Behavior and Back Again, Hanser, Munich (2006).
(55.) A.C. Papanastasiou, L.E. Scriven, and C.W. Macosco, J. Rheol., 27, 387(1983).
(56.) X.-L. Luo and R.I. Tanner, Int. J. Num. Methods. Eng., 25, 9 (1988).
(57.) T. Kajiwara, G. Barakos, and E. Mitsoulis, Int. J. Polym. Anal. Charact., 1, 201 (1995).
(58.) X.-L. Luo and R.I. Tanner, Rheol. Acta, 26, 499 (1987).
(59.) J. Meissncr, Pure Appl. Chem., 42, 551 (1975).
(60.) S.M. Alaie and T.C. Papanastasiou, Int. Polym. Proc, 8, 51 (1993).
(61.) M. Beaulne and E. Mitsoulis, J. Appl. Polym. Sci., 105, 2098 (2007).
(62.) G. Barakos and E. Milsoulis, J. Non-Newtonian Fluid Mech., 62, 55 (1996).
(63.) G.W.M. Peters and F.P.T. Baaijens, J. Non-Newtonian Fluid Mech., 68, 205 (1997).
(64.) S.G. Hatzikiriakos and J.M. Dealy, Polym. Eng. Sci., 34, 493(1994).
(65.) R.B. Bird, R.C. Armstrong, and O. Hassager, Dynamics of Polymeric Liquids, Vol. 1: Fluid Mechanics. 2nd ed., Wiley, New York (1987).
(66.) A. Hannachi and E. Mitsoulis, Adv. Polym. Technol., 12, 217 (1993).
(67.) X.-L. Luo and E. Mitsoulis, Int. J. Num. Methods. Fluids, 11, 1015 (1990).
(68.) V.-H. Karl, Angew. Makromol. Chem., 79, 11 (1979).
(69.) S.E. Kadijk and B.H.A.A. van den Brule, Polym. Eng. Sci., 34, 1535 (1994).
(70.) R.C. Penwell, R.S. Porter, and S. Middleman. J. Polym. Sci., 4-2, 9, 731 (1971).
(71.) H.M. Laun, Rheol. Acta, 22, 171 (1983).
(72.) G. Hay, M.E. Maekay, K.M. Awati, and Y. Park, J. Rheol., 43, 1099 (1999).
(73.) A. Goubert, J. Vermanl, P. Moldenaers, A. Gottfert, and B. Ernst, Appl. Rheol., 11, 26 (2001).
Correspondence to: Evan Mitsoulis: e-mail: firstname.lastname@example.org Contract grant sponsor: Natural Sciences and Engineering Research Council (NSERC) of Canada; contract grant sponsor: NTUA: contract grant number: PEBE 2009-2011.
Published online in Wiley Online Library (wileyonlinelibary.com).
[C] 2012 Society of Plastics Engineers
Mahmoud Ansari, (1) Thanasis Zisis, (2) Savvas G. Hatzikiriakos, (1) Evan Mitsoulis (2)
(1) Department of Chemical and Biological Engineering, The University of British Columbia, Vancouver, BC, Canada
(2) School of Mining Engineering and Metallurgy, National Technical University of Athens, Zografou, Athens, Greece
|Printer friendly Cite/link Email Feedback|
|Author:||Ansari, Mahmoud; Zisis, Thanasis; Hatzikiriakos, Savvas G.; Mitsoulis, Evan|
|Publication:||Polymer Engineering and Science|
|Date:||Mar 1, 2012|
|Previous Article:||The role of filler network in nonlinear viscoelastic behavior of vapor grown carbon nanofiber filled polystyrene: a strain dependent rheological...|
|Next Article:||Curing of epoxy/carbon nanotubes physical networks.|