# Numerical investigation and visualization of polymer behavior within a capillary die.

INTRODUCTION

The development of constitutive models to describe polymer behavior has made significant advances over the last two decades and entangled polymer dynamics can now be accurately predicted over a diverse range of flow regimes [1]. In recent years research has shifted direction in order to develop methodologies that allow use of these equations in practical industrial situations, where high shear rates and polydisperse polymers are common [2, 3]. However, models that accurately describe polymer behavior in these conditions are computationally very expensive with little commercial industrial application. The need for constitutive equations that capture the correct dynamics but that are computationally tractable has attracted much research [4-9]. A recently developed approach is the decoupling of chain orientation from chain stretch; an approach used by the recent pure reptation-based models of Ianniruberto and Marrucci [10, 11] and tested by Wap-perom et al. [12]. None the less, the latest derivations of these models have had little scrutiny to evaluate their performance and their applicability for use in real industrial flow situations. One established model that uses the approach of separating the two timescales and has been used both computationally and experimentally is the Remmelgas-Harrison-Leal (RHL) model [13].

This article publishes results for the RHL model in predicting the behavior of a polymer melt in a capillary rheometer. In particular, visual insights into the type of flow, described by the flow-type parameter, in the region of the barrel and the capillary die are presented. This article describes the RHL model and the necessary equations and CFD techniques that were used to perform analysis of the capillary rheometer.

THE RHL MODEL

The RHL model [13] is a vector-based representation of the DEMG [14] reptation model. The RHL model assumes that the flow-induced orientation and stretch of the polymer chains occurs at timescales that are so diverse that they can effectively be decoupled from one another. Using this approach, the RHL model produces both an independent statistical orientation distribution describing the dynamics of chain orientation and an independent stretch equation describing the dynamics of chain stretch. The orientation of the polymer chains in the RHL model is identical to an equation describing the dynamics of a suspension of rigid Brownian vectors. The second moment of the orientation distribution for a unit vector, n, is then:

[delta]<nn>/[delta]t = -2[nabla]u : <nn> <nn> - [1/[[[tau].sub.d]<[R.sup.2]>]][<nn> - [1/3]I] (1)

where, [delta]/[delta]t is the incompressible codeformational time derivative and u is the velocity gradient tensor:

[delta]<nn>/[delta]t [equivalent to] [[delta]<nn>/[delta]t] + u x [nabla]<nn> - [nabla] x (<nn>u) - [nabla][u.sup.T] x <nn> (2)

Included in Eq. 1 is a time scale analogous to the disengagement time in the DEMG model, [[tau].sub.d]. The length of the vector is included in the orientation by the effective disengagement time, [[tau].sub.d]<[R.sup.2]>. As Brownian forces are present, the orientation of the vector not only depends on the flow but also upon its length [15].

The RHL model assumes that the relaxation of stretch, based upon the longest Rouse time scale [[tau].sub.r], and the disengagement time occur at widely varied time scales. The scalar stretch function, R is equivalent to the end to end length of the polymer scaled by the equilibrium chain length. The equation of the evolution of length of the vector, R, is then given as the mean square value <[R.sup.2]>:

[D<[R.sup.2]>]/Dt = 2[nabla]u : <nn><[R.sup.2]> - [1/[[tau].sub.r]][f(<[R.sup.2]>)<[R.sup.2]> - 1] (3)

where, the extensibility of the vector, R, is modeled by a nonlinear spring function, f(<[R.sup.2]>) which is a Pade approximation to the inverse Langevin function [16]:

f(<[R.sup.2]>) = [1 - [1/3][[R.sup.2]/[[n.sub.t]/[N.sub.e]]]]/[1 - [[R.sup.2]/[[n.sub.t]/[N.sub.e]]]] (4)

and D/Dt is the convected time derivative:

[D<[R.sup.2]>]/Dt [equivalent to] [[[delta]<[R.sup.2]>]/[delta]t] + u x [nabla]<[R.sup.2]> (5)

As with the DEMG model, the maximum stretch of the polymer chain, <[R.sup.2]>[.sub.max.sup.1/2] is modeled as [square root of ([n.sub.t]/[N.sub.e])]. Where, [N.sub.e] is the number of entanglements per chain and [n.sub.t] is the number of Kuhn steps per chain. Equations 1-4 define the orientation and the deformation of a vector that describes the polymer chain. To complete the constitutive model, an expression is needed for the polymer contribution to the stress tensor. The stress tensor, [sigma], is assumed to be purely elastic and factorizable in the form:

[sigma] = [[G.sub.N.sup.0]/[[tau].sub.d]]f(<[R.sup.2]>)<[R.sup.2]><nn> (6)

The two characteristic time scales are considered to be related by the number of entanglements per chain [13]:

[[tau].sub.d] = 3[N.sub.e][[tau].sub.r] (7)

Therefore, the model only has three independent parameters in order to successfully describe a known polymer melt:

* the Rouse time, [[tau].sub.r].

* the number of entanglements, [N.sub.e].

* the polymer contribution to stress, [G.sub.N.sup.0].

NUMERICAL SOLUTION

In this section we present the inclusion of the polymer model in the continuum equations describing the flow of the polymer melt, and discuss their solution using the finite volume CFD technique.

Governing Equations

The governing equations of motion for an incompressible flow are given by the nondimensional Navier-Stokes Equations [17]:

[nabla] x u = 0 (8)

[Du/Dt] - [1/Re][[nabla].sup.2]u = -[nabla]p (9)

where the characteristic Reynolds number is given by:

Re = [[rho][u.sub.c][l.sub.c]]/[[eta].sub.0] (10)

where, [u.sub.c] is a characteristic velocity scale, [l.sub.c] is a characteristic length scale, and [[eta].sub.0] is the polymer zero shear viscosity. The RHL model polymer contribution to the stress, Eq. 6, is included in the Cauchy equation of motion so that Eq. 9 becomes:

[Du/Dt] - [1/Re][[nabla].sup.2]u = -[nabla]p + [sigma] (11)

A finite volume-based computational fluid dynamics package called openFoam [18] was chosen to solve the Equations 1-11 and model the conditions within the capillary die using the RHL polymer model. The source code of openFoam is available with open architecture and thus can be specifically adapted to solve a specific problem. The optimized openFoam CFD code is an efficient method of solving the characteristic polymer equations as well as providing a framework for producing the necessary preprocessing and postprocessing necessary to model the geometries under investigation in this study. The resulting numerical code has been named "PolyFoam."

Spatial Discretization

The results of the RHL model were determined from the computational modeling of a capillary rheometer. The problem domain was simplified to an axi-symmetric case with the same geometry as the long die of a capillary rheometer as shown in Fig. 1.

An axi-symmetric mesh was generated in order to reduce the number of control volumes that would be needed to model the full three-dimensional domain of the capillary. This is accomplished in the three-dimensional openFoam solver using a wedge geometry that is one-cell thick. As the greatest area of change would be at the wall of the capillary, the mesh is graded more finely in this area. It was found that the solution converges with a resolution of 20 control volumes across the radius of the die. Therefore, for the structured meshes used in the study of the capillary die, a resolution of 20 cells across the die radius was used.

[FIGURE 1 OMITTED]

Boundary and Initial Conditions

The computational mesh consists of faces that coincide with the boundaries of the physical domain. Solution of the Navier-Stokes equations requires the velocity and the pressure at the boundaries to be explicitly prescribed [19]. To correctly simulate the conditions within the capillary die, material conditions at these boundaries also need to be specified:

* Inlet Boundary. To model the piston forcing the melt through the barrel a constant velocity profile was prescribed at the inlet boundary.

* Outlet Boundary. The outlet boundary is specified as being at atmospheric pressure.

* Symmetry Plane Boundary. Because of the symmetry of the capillary die a symmetry plane boundary condition can be implemented along the central axis of the die.

* No-Slip Walls. Using the no-slip wall assumption a fixed velocity value of zero was used at the boundary face.

* Wedge Boundary. For 2D cases where no solution is required on parallel planes normal to the 3rd dimension, the flux leaving one parallel boundary is modeled as entering normal to the other parallel boundary. For a 2D axi-symmetric case, where the boundaries are not parallel, a wedge boundary condition is used to ensure the correct modeling of the flux between boundaries.

As well as the physical boundary conditions there are also the initial conditions of the material within the capillary die. Before the piston starts to force the polymer melt through the die, the polymer chains are in their equilibrium state. For this to be the case, the initial condition of the orientation distribution is isotropic with no preferred orientation, as defined in Equation 12.

<nn> = [1/3]I (12)

The equilibrium unstretched state of the polymer is defined as:

<[R.sup.2]> = [^.[beta]] (13)

where [^.[beta]] is the solution to the equation:

f(<[R.sup.2]>)<[R.sup.2]> = 1 (14)

The value of <[R.sup.2]> is the length of the vector used to describe stretch, scaled by its equilibrium value. Therefore the value of <[R.sup.2]> should be unity as its initial equilibrium value. However, the spring law function causes the value of <[R.sup.2]> to equal [^.[beta]] at unity. For the parameter values used in this study, [^.[beta]] [approximately equal to] 0.9. This is a common problem and noticed in other polymer constitutive equations using a nonlinear spring law [20]. The RHL model parameters used were typical for a HDPE with [[tau].sub.r] = 1.2 x [10.sup.-3] s, [N.sub.e] = 190, and [G.sub.N.sup.0] = 1 x [10.sup.5] Pa.

Discretization Schemes and Solution Algorithm

A central differencing scheme was deemed most appropriate for use in this study. This scheme interpolates between two cell centers to find an average value of the flux at the cell face. All interpolation schemes are an approximation, with the central differencing approach being considered second order accurate. A common discretization scheme used in polymer computational rheology is upwind differencing [20]. With this scheme, the interpolation of the face values is found by taking the value of the upwind cell centre. Upwind differencing is used due to the ease of implementation and stability of convergence. However, it is first order accurate and certain effects such as chain stretch overshoots can be lost from the solution. Upwind differencing is also considered to introduce a numerical viscosity [21], which stabilizes the solution at the expense of altering the rheology.

The PISO (Pressure Implicit Splitting of Operators) [22] solution algorithm was to solve the Navier-Stokes equations for transient incompressible flow. The PISO algorithm was suitably modified to include the equations needed for the orientation and stretch of the polymer constitutive equations, along with polymer contribution to the stress tensor. The PISO algorithm was chosen as it can be used for time-dependent flows. Another solution algorithm that can be used for steady state problems is the SIMPLE algorithm (Semi-Implicit Method for Pressure Linked Equations). Comparison of the SIMPLE and PISO algorithms has been conducted by Jang et al. [23]. Although the transient behavior of the flow has not been studied, it is hoped that the code will provide a useful tool for further investigations. The SIMPLE algorithm might have been able to provide a quicker solution but is generally used for steady-state examples. Therefore, to facilitate future work the PISO algorithm was implemented.

RESULTS

Several different numerical parameters can be derived from the computational data to give useful insight into the flow characteristics. Comparisons with classical ideas of flow within the capillary are compared with information taken from literature and experimental results.

Flow Type Parameter

The first set of steady state predictions from the RHL model is of the flow type parameter. The velocity gradient for two-dimensional planar flow, in an appropriate coordinate system, can be written as:

[FIGURE 2 OMITTED]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

For a velocity gradient tensor in arbitrary coordinate system, the flow parameter, [lambda], can be found by:

[lambda] = [|D| - |[omega]|]/[|D| + |[omega]|] (16)

where, D is the rate of deformation tensor and [omega] is the vorticity tensor. The flow parameter has values between rotational flow where [lambda] = -1, simple shear where [lambda] = 0, and planar extensional flow where [lambda] = 1. The flow type parameter is a measure of extension to vorticity. Simple shear can be represented by superposition of a purely rotational flow and a purely extensional flow.

Although the geometric mesh only modeled half of the capillary die due to symmetry, the full result is shown for clarity. Figure 2 shows the flow type parameter, [lambda], for a capillary wall shear rate of 10 [s.sup.-1].

The flow type parameter gives useful information about the flow regime the fluid is experiencing. The region of strong elongational flow at the entrance to the barrel illustrates the method with which the boundary conditions for the piston have been modeled. Modeling the piston, as experimentally observed, would mean the use of a dynamic mesh to model the polymer melt in the barrel as it forced through the die. This would unnecessarily complicate the program code and increase the computational time. A fixed mesh is therefore prescribed along with a constant velocity profile as produced directly in front of a piston. Although not having direct relevance to the modeling of the region of the die, studying this area produces an interesting interpretation of the flow.

The region in front of the piston shows the behavior of the flow during the transition from a constant velocity profile, as specified by the piston movement, to a plug flow velocity of a shear thinning polymer. It is clear that certain areas within this region either have to accelerate or decelerate to match the developed flow profile. The melt nearest the wall has to reduce in velocity, as there is a no-slip boundary condition at the barrel walls. The necessary reduction in velocity from the constant profile to the fully developed flow profile is non linear. To conserve the equations of momentum and conservation of mass, a velocity gradient is formed that causes rotational flow. The melt in the central region of the die has to accelerate in order to match the developed flow velocity profile. The combination of the flow accelerating and rotating causes strong elongational flow in the central region of the die.

Examining the flow type predictions directly before the entrance of the die shows a region of strongly extensional flow. In this area, the streamlines of the polymer melt have to accelerate and converge in order to enter the die. This causes the extensional flow characteristics. A pressure drop caused by the elongational flow in this region is responsible for the Bagley ends correction being used. The elongational flow region extends far back into the barrel region showing the strong effects the sudden restriction of the die has to the flow. The flow type parameter predictions also show a rotational region in the corners of the sudden constriction. The flow within the die region is shown to be shear in nature and this confirms the simple flow characteristics that are necessary to determine the viscosity of the flow.

The results of the flow type parameter for the higher shear rate of 1000 [s.sup.-1] is shown in Fig. 3. Overall, the higher shear flow type prediction does not show significant variation from the 10 [s.sup.-1] predictions. Regions appear more pronounced as the magnitude of the flow is considerably larger. The region within the die is no longer shear alone and contains elongational characteristics. This effect is a manifestation of transient overshoots in stretch and has been noted both in other computational models [2] and observed in birefringence experiments of other polymer melts [24]. The elongational region in front of the die entrance no longer extends to the point of the sudden constriction, as with the 10 [s.sup.-1] sample. This is caused by more material being drawn into the die from the recirculating areas of the corners of the die entrance.

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

The higher shear rate pressure distribution result of Fig. 4 shows a large pressure drop as experimentally observed by Han [25]. This pressure drop is caused by the melt having to accelerate and rearrange as it enters the die. Han observed a typical 35% variation from the pressure in the barrel to that just inside the die. The predictions of the 1000 [s.sup.-1] example in Fig. 4 shows a similar pressure variation to that experimentally observed by Han. The pressure drop observed by the ends of the capillary can significantly alter the driving pressure and is why the Bagley correction used in capillary rheometry. The reason for this effect can be seen more clearly by examining the plots of Figs. 2 and 3 showing the flow type parameter and the strong extensional flow at the entrance to the die.

Figure 4 also shows the pressure distribution at the entrance to the die for the low shear example. The result is consistent with a Newtonian liquid, showing a constant pressure at the entrance and then a pressure drop along the die as it reaches atmospheric pressure at the exit. At this low shear example, the extensional flow at the entrance to the die and the shear within the die is small enough in magnitude that the polymer melt has time to relax from stress. Therefore, the majority of the melt has a viscosity that is equal to the zero shear viscosity and is the reason why the observed pressure drop is similar to a Newtonian fluid. Although the flow directly in front of the capillary region shows strong elongational flow, as shown by Fig. 2, it is also important to examine the magnitude of the deformation as discussed in the next section.

Magnitude of the Deformation Tensor

The results of the flow type parameter show the elongational flow created in front of the piston could prove an interesting experimental method to study elongational flow characteristics. However, the flow type parameter does not show the magnitude of the different flow types. The radius of the capillary die is small in order to cause high shear within this region. In comparison, the deformation of the material in the larger barrel radius is considerably smaller. To illustrate the large range of the deformation experience by the melt, it was necessary to plot the magnitude of the deformation tensor on a logarithmic scale, as shown by Figs. 5 and 6.

[FIGURE 5 OMITTED]

As expected, the plots show that the region of highest deformation is in the capillary. The purpose of this illustration was to show the magnitude of the elongational flow in the proximity of the piston. Both Figs. 5 and 6 show that the deformation in this region tend to be around 1000 times smaller than the shear at the walls of the capillary. Although experimental investigation into this area is a valid method to study elongational flow characteristics, the magnitude of the elongational flow is not very high compared with other techniques. To obtain a large range of rates, the piston would have to move very quickly causing high pressures within the capillary barrel and large material waste.

Polymer Orientation

The orientation of the polymer is described by a tensor whose eigenvectors and eigenvalues describe the orientation distribution function as an ellipsoid. This 3D representation is difficult to visualize in 2D. Many papers comparing predictions of model behavior choose to represent polymer orientation as the angle to the flow direction of the largest eigenvalue [26]. Using this method, chains that are principally aligned in the flow direction have value of 0[degrees], whilst those aligned perpendicular to the direction of flow take a value of 90[degrees]. For the value where the orientation distribution has no preferred direction, the tensor presentation is a perfect sphere where all eigenvectors are at 45[degrees] to the flow direction. Figures 7 and 8 show a 2D representation of the orientation of the polymer melt using this method.

[FIGURE 6 OMITTED]

Figure 7 shows the angle the largest eigenvalue makes with the flow direction for a capillary shear rate of 10 [s.sup.-1]. In the region of high extensional flow the polymer melt aligns with the flow direction in order to converge and accelerate with the other streamlines entering the die. Within the die region, there is no preferred alignment except in the high shear region immediately adjacent to the wall. Therefore for this low shear rate, the polymer melt has time to relax stress, as the magnitude of the flow is not high enough to orientate the chains in a preferred direction. On first examination, the region near the corners of the sudden constriction show counter-intuitive results. The orientation suddenly flips from being perpendicular to the flow, to being aligned with the flow. This is a consequence of the method with which a 3D orientation tensor is displayed in 2D by using the angle the largest eigenvalue makes with the flow direction. In this region a change of the eigenvalue that is greatest in magnitude is observed causing the sudden change in orientation.

The PolyFoam predictions for the orientation at the higher shear rate of 1000 [s.sup.-1] is shown in Fig. 8. The orientation outside of the die shows a very similar pattern to the 10 [s.sup.-1] wall shear rate example. The region of the die, however, shows a considerably different situation. Where the flow was not great enough to show a preferred orientation of the polymer melt within the die of the 10 [s.sup.-1] predictions, the higher wall shear rate predictions show the orientation is aligned with the flow. The high flow rate in the die prevents the chains relaxing completely by reptation dynamics.

[FIGURE 7 OMITTED]

CONCLUSIONS

The implementation of the RHL constitutive model has produced a number of predictions that provide novel insights into the flow regime within a capillary rheometer. The operation of the capillary rheometer is designed to obtain simple shear of the melt within the capillary die. The viscosity can then be calculated using appropriate equations. However, practical rheology necessitates the use of the Bagley ends correction to take account for pressure drops observed at the entrance of the die due to the strong elongational flow in this region. The PolyFoam predictions for the flowtype parameter have revealed the reasons for the pressure drop experienced in the entrance. Predictions of the pressure drop have shown to be similar to that of experimentally measured behavior [25].

[FIGURE 8 OMITTED]

Investigations into the flow regime directly in front of the piston within the capillary rheometer have shown an area of strong elongational flow. The magnitude of the flow is considerably smaller than the shear flow present in the capillary but could provide a useful method to examine the behavior of polymer melts in elongational flow.

REFERENCES

1. L.G. Leal and J.P. Oberhauser, Rheol. J., 12, 1 (2000).

2. T.C.B. McLeish, Chem. Eng. Res. Des. Trans. Inst. Chem. Eng. A, 78, 12 (2000).

3. C. Pattamaprom and R.G. Larson, Rheol. Acta, 40, 516 (2001).

4. M.J. Crochet and K. Walters, Annu. Rev. Fluid Mech., 15, 241 (1983).

5. M.J. Crochet, Rubber Chem. Technol., 62, 426 (1989).

6. R. Keunings, Rheol. Acta, 29, 556 (1990).

7. F.P.T. Baaijens, J. Non-Newtonian Fluid Mech., 79, 361 (1998).

8. R.I. Tanner and S.C. Xue, Korea--Aust. Rheol. J., 14, 143 (2002).

9. R.S. Graham, A.E. Likhtman, C.B. McLeish, and S.T. Milner, J. Rheol., 47, 1171 (2003).

10. G. Ianniruberto and G. Marrucci, J. Non-Newtonian Fluid Mech., 102, 383 (2002).

11. G. Ianniruberto and G. Marrucci, J. Rheol., 45, 1305 (2001).

12. P. Wapperom, R. Keunings, and G. Ianniruberto, J. Rheol., 47, 247 (2002).

13. J. Remmelgas, G. Harrison, and L. G. Leal, J. Non-Newtonian Fluid Mech., 80, 115 (1999).

14. G. Marrucci and N. Grizzuti, Gaz. Chim. Ital., 118, 179 (1988).

15. W.L. Olbricht, J.M. Rallison, and L.G. Leal, J. Non-Newtonian Fluid Mech., 10, 291 (1982).

16. A. Cohen, Rheol. Acta, 30, 270 (1991).

17. R. Aris, Vectors, Tensors, and the Basic Equations of Fluid Mechanics, Dover, New York (1989).

18. openFoam, www.opencfd.co.uk.

19. C. Hirsch, Numerical Computation of Internal and External Flows, Wiley, New York (1991).

20. J. Remmelgas, P. Singh, and L.G. Leal, J. Non-Newtonian Fluid Mech., 88, 31 (1999).

21. H.K. Versteeg and W. Malalasekera, An Introduction to Computational Fluid Dynamics, Longman, England (1995).

22. R.I. Issa, J. Comp. Phys., 62, 40 (1986).

23. D.S. Jang, R. Jetli, and S. Acharya, Numer. Heat Transfer, 19, 209 (1986).

24. R. Ahmed, R.F. Liang, and M.R. Mackley, J. Non-Newtonian Fluid Mech., 59, 129 (1995).

25. C.D. Han, Rheology in Polymer Processing, Academic Press, New York (1976).

26. D.W. Mead, D. Yavich, and L.G. Leal, Rheol. Acta, 34, 360 (1995).

M.A. Beard, G.R. Tabor, S.J.K. Ritchie, K.E. Evans

School of Engineering, Computer Science and Mathematics, University of Exeter, UK

Correspondence to: M.A. Beard; e-mail: m.a.beard@ex.ac.uk

The development of constitutive models to describe polymer behavior has made significant advances over the last two decades and entangled polymer dynamics can now be accurately predicted over a diverse range of flow regimes [1]. In recent years research has shifted direction in order to develop methodologies that allow use of these equations in practical industrial situations, where high shear rates and polydisperse polymers are common [2, 3]. However, models that accurately describe polymer behavior in these conditions are computationally very expensive with little commercial industrial application. The need for constitutive equations that capture the correct dynamics but that are computationally tractable has attracted much research [4-9]. A recently developed approach is the decoupling of chain orientation from chain stretch; an approach used by the recent pure reptation-based models of Ianniruberto and Marrucci [10, 11] and tested by Wap-perom et al. [12]. None the less, the latest derivations of these models have had little scrutiny to evaluate their performance and their applicability for use in real industrial flow situations. One established model that uses the approach of separating the two timescales and has been used both computationally and experimentally is the Remmelgas-Harrison-Leal (RHL) model [13].

This article publishes results for the RHL model in predicting the behavior of a polymer melt in a capillary rheometer. In particular, visual insights into the type of flow, described by the flow-type parameter, in the region of the barrel and the capillary die are presented. This article describes the RHL model and the necessary equations and CFD techniques that were used to perform analysis of the capillary rheometer.

THE RHL MODEL

The RHL model [13] is a vector-based representation of the DEMG [14] reptation model. The RHL model assumes that the flow-induced orientation and stretch of the polymer chains occurs at timescales that are so diverse that they can effectively be decoupled from one another. Using this approach, the RHL model produces both an independent statistical orientation distribution describing the dynamics of chain orientation and an independent stretch equation describing the dynamics of chain stretch. The orientation of the polymer chains in the RHL model is identical to an equation describing the dynamics of a suspension of rigid Brownian vectors. The second moment of the orientation distribution for a unit vector, n, is then:

[delta]<nn>/[delta]t = -2[nabla]u : <nn> <nn> - [1/[[[tau].sub.d]<[R.sup.2]>]][<nn> - [1/3]I] (1)

where, [delta]/[delta]t is the incompressible codeformational time derivative and u is the velocity gradient tensor:

[delta]<nn>/[delta]t [equivalent to] [[delta]<nn>/[delta]t] + u x [nabla]<nn> - [nabla] x (<nn>u) - [nabla][u.sup.T] x <nn> (2)

Included in Eq. 1 is a time scale analogous to the disengagement time in the DEMG model, [[tau].sub.d]. The length of the vector is included in the orientation by the effective disengagement time, [[tau].sub.d]<[R.sup.2]>. As Brownian forces are present, the orientation of the vector not only depends on the flow but also upon its length [15].

The RHL model assumes that the relaxation of stretch, based upon the longest Rouse time scale [[tau].sub.r], and the disengagement time occur at widely varied time scales. The scalar stretch function, R is equivalent to the end to end length of the polymer scaled by the equilibrium chain length. The equation of the evolution of length of the vector, R, is then given as the mean square value <[R.sup.2]>:

[D<[R.sup.2]>]/Dt = 2[nabla]u : <nn><[R.sup.2]> - [1/[[tau].sub.r]][f(<[R.sup.2]>)<[R.sup.2]> - 1] (3)

where, the extensibility of the vector, R, is modeled by a nonlinear spring function, f(<[R.sup.2]>) which is a Pade approximation to the inverse Langevin function [16]:

f(<[R.sup.2]>) = [1 - [1/3][[R.sup.2]/[[n.sub.t]/[N.sub.e]]]]/[1 - [[R.sup.2]/[[n.sub.t]/[N.sub.e]]]] (4)

and D/Dt is the convected time derivative:

[D<[R.sup.2]>]/Dt [equivalent to] [[[delta]<[R.sup.2]>]/[delta]t] + u x [nabla]<[R.sup.2]> (5)

As with the DEMG model, the maximum stretch of the polymer chain, <[R.sup.2]>[.sub.max.sup.1/2] is modeled as [square root of ([n.sub.t]/[N.sub.e])]. Where, [N.sub.e] is the number of entanglements per chain and [n.sub.t] is the number of Kuhn steps per chain. Equations 1-4 define the orientation and the deformation of a vector that describes the polymer chain. To complete the constitutive model, an expression is needed for the polymer contribution to the stress tensor. The stress tensor, [sigma], is assumed to be purely elastic and factorizable in the form:

[sigma] = [[G.sub.N.sup.0]/[[tau].sub.d]]f(<[R.sup.2]>)<[R.sup.2]><nn> (6)

The two characteristic time scales are considered to be related by the number of entanglements per chain [13]:

[[tau].sub.d] = 3[N.sub.e][[tau].sub.r] (7)

Therefore, the model only has three independent parameters in order to successfully describe a known polymer melt:

* the Rouse time, [[tau].sub.r].

* the number of entanglements, [N.sub.e].

* the polymer contribution to stress, [G.sub.N.sup.0].

NUMERICAL SOLUTION

In this section we present the inclusion of the polymer model in the continuum equations describing the flow of the polymer melt, and discuss their solution using the finite volume CFD technique.

Governing Equations

The governing equations of motion for an incompressible flow are given by the nondimensional Navier-Stokes Equations [17]:

[nabla] x u = 0 (8)

[Du/Dt] - [1/Re][[nabla].sup.2]u = -[nabla]p (9)

where the characteristic Reynolds number is given by:

Re = [[rho][u.sub.c][l.sub.c]]/[[eta].sub.0] (10)

where, [u.sub.c] is a characteristic velocity scale, [l.sub.c] is a characteristic length scale, and [[eta].sub.0] is the polymer zero shear viscosity. The RHL model polymer contribution to the stress, Eq. 6, is included in the Cauchy equation of motion so that Eq. 9 becomes:

[Du/Dt] - [1/Re][[nabla].sup.2]u = -[nabla]p + [sigma] (11)

A finite volume-based computational fluid dynamics package called openFoam [18] was chosen to solve the Equations 1-11 and model the conditions within the capillary die using the RHL polymer model. The source code of openFoam is available with open architecture and thus can be specifically adapted to solve a specific problem. The optimized openFoam CFD code is an efficient method of solving the characteristic polymer equations as well as providing a framework for producing the necessary preprocessing and postprocessing necessary to model the geometries under investigation in this study. The resulting numerical code has been named "PolyFoam."

Spatial Discretization

The results of the RHL model were determined from the computational modeling of a capillary rheometer. The problem domain was simplified to an axi-symmetric case with the same geometry as the long die of a capillary rheometer as shown in Fig. 1.

An axi-symmetric mesh was generated in order to reduce the number of control volumes that would be needed to model the full three-dimensional domain of the capillary. This is accomplished in the three-dimensional openFoam solver using a wedge geometry that is one-cell thick. As the greatest area of change would be at the wall of the capillary, the mesh is graded more finely in this area. It was found that the solution converges with a resolution of 20 control volumes across the radius of the die. Therefore, for the structured meshes used in the study of the capillary die, a resolution of 20 cells across the die radius was used.

[FIGURE 1 OMITTED]

Boundary and Initial Conditions

The computational mesh consists of faces that coincide with the boundaries of the physical domain. Solution of the Navier-Stokes equations requires the velocity and the pressure at the boundaries to be explicitly prescribed [19]. To correctly simulate the conditions within the capillary die, material conditions at these boundaries also need to be specified:

* Inlet Boundary. To model the piston forcing the melt through the barrel a constant velocity profile was prescribed at the inlet boundary.

* Outlet Boundary. The outlet boundary is specified as being at atmospheric pressure.

* Symmetry Plane Boundary. Because of the symmetry of the capillary die a symmetry plane boundary condition can be implemented along the central axis of the die.

* No-Slip Walls. Using the no-slip wall assumption a fixed velocity value of zero was used at the boundary face.

* Wedge Boundary. For 2D cases where no solution is required on parallel planes normal to the 3rd dimension, the flux leaving one parallel boundary is modeled as entering normal to the other parallel boundary. For a 2D axi-symmetric case, where the boundaries are not parallel, a wedge boundary condition is used to ensure the correct modeling of the flux between boundaries.

As well as the physical boundary conditions there are also the initial conditions of the material within the capillary die. Before the piston starts to force the polymer melt through the die, the polymer chains are in their equilibrium state. For this to be the case, the initial condition of the orientation distribution is isotropic with no preferred orientation, as defined in Equation 12.

<nn> = [1/3]I (12)

The equilibrium unstretched state of the polymer is defined as:

<[R.sup.2]> = [^.[beta]] (13)

where [^.[beta]] is the solution to the equation:

f(<[R.sup.2]>)<[R.sup.2]> = 1 (14)

The value of <[R.sup.2]> is the length of the vector used to describe stretch, scaled by its equilibrium value. Therefore the value of <[R.sup.2]> should be unity as its initial equilibrium value. However, the spring law function causes the value of <[R.sup.2]> to equal [^.[beta]] at unity. For the parameter values used in this study, [^.[beta]] [approximately equal to] 0.9. This is a common problem and noticed in other polymer constitutive equations using a nonlinear spring law [20]. The RHL model parameters used were typical for a HDPE with [[tau].sub.r] = 1.2 x [10.sup.-3] s, [N.sub.e] = 190, and [G.sub.N.sup.0] = 1 x [10.sup.5] Pa.

Discretization Schemes and Solution Algorithm

A central differencing scheme was deemed most appropriate for use in this study. This scheme interpolates between two cell centers to find an average value of the flux at the cell face. All interpolation schemes are an approximation, with the central differencing approach being considered second order accurate. A common discretization scheme used in polymer computational rheology is upwind differencing [20]. With this scheme, the interpolation of the face values is found by taking the value of the upwind cell centre. Upwind differencing is used due to the ease of implementation and stability of convergence. However, it is first order accurate and certain effects such as chain stretch overshoots can be lost from the solution. Upwind differencing is also considered to introduce a numerical viscosity [21], which stabilizes the solution at the expense of altering the rheology.

The PISO (Pressure Implicit Splitting of Operators) [22] solution algorithm was to solve the Navier-Stokes equations for transient incompressible flow. The PISO algorithm was suitably modified to include the equations needed for the orientation and stretch of the polymer constitutive equations, along with polymer contribution to the stress tensor. The PISO algorithm was chosen as it can be used for time-dependent flows. Another solution algorithm that can be used for steady state problems is the SIMPLE algorithm (Semi-Implicit Method for Pressure Linked Equations). Comparison of the SIMPLE and PISO algorithms has been conducted by Jang et al. [23]. Although the transient behavior of the flow has not been studied, it is hoped that the code will provide a useful tool for further investigations. The SIMPLE algorithm might have been able to provide a quicker solution but is generally used for steady-state examples. Therefore, to facilitate future work the PISO algorithm was implemented.

RESULTS

Several different numerical parameters can be derived from the computational data to give useful insight into the flow characteristics. Comparisons with classical ideas of flow within the capillary are compared with information taken from literature and experimental results.

Flow Type Parameter

The first set of steady state predictions from the RHL model is of the flow type parameter. The velocity gradient for two-dimensional planar flow, in an appropriate coordinate system, can be written as:

[FIGURE 2 OMITTED]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

For a velocity gradient tensor in arbitrary coordinate system, the flow parameter, [lambda], can be found by:

[lambda] = [|D| - |[omega]|]/[|D| + |[omega]|] (16)

where, D is the rate of deformation tensor and [omega] is the vorticity tensor. The flow parameter has values between rotational flow where [lambda] = -1, simple shear where [lambda] = 0, and planar extensional flow where [lambda] = 1. The flow type parameter is a measure of extension to vorticity. Simple shear can be represented by superposition of a purely rotational flow and a purely extensional flow.

Although the geometric mesh only modeled half of the capillary die due to symmetry, the full result is shown for clarity. Figure 2 shows the flow type parameter, [lambda], for a capillary wall shear rate of 10 [s.sup.-1].

The flow type parameter gives useful information about the flow regime the fluid is experiencing. The region of strong elongational flow at the entrance to the barrel illustrates the method with which the boundary conditions for the piston have been modeled. Modeling the piston, as experimentally observed, would mean the use of a dynamic mesh to model the polymer melt in the barrel as it forced through the die. This would unnecessarily complicate the program code and increase the computational time. A fixed mesh is therefore prescribed along with a constant velocity profile as produced directly in front of a piston. Although not having direct relevance to the modeling of the region of the die, studying this area produces an interesting interpretation of the flow.

The region in front of the piston shows the behavior of the flow during the transition from a constant velocity profile, as specified by the piston movement, to a plug flow velocity of a shear thinning polymer. It is clear that certain areas within this region either have to accelerate or decelerate to match the developed flow profile. The melt nearest the wall has to reduce in velocity, as there is a no-slip boundary condition at the barrel walls. The necessary reduction in velocity from the constant profile to the fully developed flow profile is non linear. To conserve the equations of momentum and conservation of mass, a velocity gradient is formed that causes rotational flow. The melt in the central region of the die has to accelerate in order to match the developed flow velocity profile. The combination of the flow accelerating and rotating causes strong elongational flow in the central region of the die.

Examining the flow type predictions directly before the entrance of the die shows a region of strongly extensional flow. In this area, the streamlines of the polymer melt have to accelerate and converge in order to enter the die. This causes the extensional flow characteristics. A pressure drop caused by the elongational flow in this region is responsible for the Bagley ends correction being used. The elongational flow region extends far back into the barrel region showing the strong effects the sudden restriction of the die has to the flow. The flow type parameter predictions also show a rotational region in the corners of the sudden constriction. The flow within the die region is shown to be shear in nature and this confirms the simple flow characteristics that are necessary to determine the viscosity of the flow.

The results of the flow type parameter for the higher shear rate of 1000 [s.sup.-1] is shown in Fig. 3. Overall, the higher shear flow type prediction does not show significant variation from the 10 [s.sup.-1] predictions. Regions appear more pronounced as the magnitude of the flow is considerably larger. The region within the die is no longer shear alone and contains elongational characteristics. This effect is a manifestation of transient overshoots in stretch and has been noted both in other computational models [2] and observed in birefringence experiments of other polymer melts [24]. The elongational region in front of the die entrance no longer extends to the point of the sudden constriction, as with the 10 [s.sup.-1] sample. This is caused by more material being drawn into the die from the recirculating areas of the corners of the die entrance.

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

The higher shear rate pressure distribution result of Fig. 4 shows a large pressure drop as experimentally observed by Han [25]. This pressure drop is caused by the melt having to accelerate and rearrange as it enters the die. Han observed a typical 35% variation from the pressure in the barrel to that just inside the die. The predictions of the 1000 [s.sup.-1] example in Fig. 4 shows a similar pressure variation to that experimentally observed by Han. The pressure drop observed by the ends of the capillary can significantly alter the driving pressure and is why the Bagley correction used in capillary rheometry. The reason for this effect can be seen more clearly by examining the plots of Figs. 2 and 3 showing the flow type parameter and the strong extensional flow at the entrance to the die.

Figure 4 also shows the pressure distribution at the entrance to the die for the low shear example. The result is consistent with a Newtonian liquid, showing a constant pressure at the entrance and then a pressure drop along the die as it reaches atmospheric pressure at the exit. At this low shear example, the extensional flow at the entrance to the die and the shear within the die is small enough in magnitude that the polymer melt has time to relax from stress. Therefore, the majority of the melt has a viscosity that is equal to the zero shear viscosity and is the reason why the observed pressure drop is similar to a Newtonian fluid. Although the flow directly in front of the capillary region shows strong elongational flow, as shown by Fig. 2, it is also important to examine the magnitude of the deformation as discussed in the next section.

Magnitude of the Deformation Tensor

The results of the flow type parameter show the elongational flow created in front of the piston could prove an interesting experimental method to study elongational flow characteristics. However, the flow type parameter does not show the magnitude of the different flow types. The radius of the capillary die is small in order to cause high shear within this region. In comparison, the deformation of the material in the larger barrel radius is considerably smaller. To illustrate the large range of the deformation experience by the melt, it was necessary to plot the magnitude of the deformation tensor on a logarithmic scale, as shown by Figs. 5 and 6.

[FIGURE 5 OMITTED]

As expected, the plots show that the region of highest deformation is in the capillary. The purpose of this illustration was to show the magnitude of the elongational flow in the proximity of the piston. Both Figs. 5 and 6 show that the deformation in this region tend to be around 1000 times smaller than the shear at the walls of the capillary. Although experimental investigation into this area is a valid method to study elongational flow characteristics, the magnitude of the elongational flow is not very high compared with other techniques. To obtain a large range of rates, the piston would have to move very quickly causing high pressures within the capillary barrel and large material waste.

Polymer Orientation

The orientation of the polymer is described by a tensor whose eigenvectors and eigenvalues describe the orientation distribution function as an ellipsoid. This 3D representation is difficult to visualize in 2D. Many papers comparing predictions of model behavior choose to represent polymer orientation as the angle to the flow direction of the largest eigenvalue [26]. Using this method, chains that are principally aligned in the flow direction have value of 0[degrees], whilst those aligned perpendicular to the direction of flow take a value of 90[degrees]. For the value where the orientation distribution has no preferred direction, the tensor presentation is a perfect sphere where all eigenvectors are at 45[degrees] to the flow direction. Figures 7 and 8 show a 2D representation of the orientation of the polymer melt using this method.

[FIGURE 6 OMITTED]

Figure 7 shows the angle the largest eigenvalue makes with the flow direction for a capillary shear rate of 10 [s.sup.-1]. In the region of high extensional flow the polymer melt aligns with the flow direction in order to converge and accelerate with the other streamlines entering the die. Within the die region, there is no preferred alignment except in the high shear region immediately adjacent to the wall. Therefore for this low shear rate, the polymer melt has time to relax stress, as the magnitude of the flow is not high enough to orientate the chains in a preferred direction. On first examination, the region near the corners of the sudden constriction show counter-intuitive results. The orientation suddenly flips from being perpendicular to the flow, to being aligned with the flow. This is a consequence of the method with which a 3D orientation tensor is displayed in 2D by using the angle the largest eigenvalue makes with the flow direction. In this region a change of the eigenvalue that is greatest in magnitude is observed causing the sudden change in orientation.

The PolyFoam predictions for the orientation at the higher shear rate of 1000 [s.sup.-1] is shown in Fig. 8. The orientation outside of the die shows a very similar pattern to the 10 [s.sup.-1] wall shear rate example. The region of the die, however, shows a considerably different situation. Where the flow was not great enough to show a preferred orientation of the polymer melt within the die of the 10 [s.sup.-1] predictions, the higher wall shear rate predictions show the orientation is aligned with the flow. The high flow rate in the die prevents the chains relaxing completely by reptation dynamics.

[FIGURE 7 OMITTED]

CONCLUSIONS

The implementation of the RHL constitutive model has produced a number of predictions that provide novel insights into the flow regime within a capillary rheometer. The operation of the capillary rheometer is designed to obtain simple shear of the melt within the capillary die. The viscosity can then be calculated using appropriate equations. However, practical rheology necessitates the use of the Bagley ends correction to take account for pressure drops observed at the entrance of the die due to the strong elongational flow in this region. The PolyFoam predictions for the flowtype parameter have revealed the reasons for the pressure drop experienced in the entrance. Predictions of the pressure drop have shown to be similar to that of experimentally measured behavior [25].

[FIGURE 8 OMITTED]

Investigations into the flow regime directly in front of the piston within the capillary rheometer have shown an area of strong elongational flow. The magnitude of the flow is considerably smaller than the shear flow present in the capillary but could provide a useful method to examine the behavior of polymer melts in elongational flow.

REFERENCES

1. L.G. Leal and J.P. Oberhauser, Rheol. J., 12, 1 (2000).

2. T.C.B. McLeish, Chem. Eng. Res. Des. Trans. Inst. Chem. Eng. A, 78, 12 (2000).

3. C. Pattamaprom and R.G. Larson, Rheol. Acta, 40, 516 (2001).

4. M.J. Crochet and K. Walters, Annu. Rev. Fluid Mech., 15, 241 (1983).

5. M.J. Crochet, Rubber Chem. Technol., 62, 426 (1989).

6. R. Keunings, Rheol. Acta, 29, 556 (1990).

7. F.P.T. Baaijens, J. Non-Newtonian Fluid Mech., 79, 361 (1998).

8. R.I. Tanner and S.C. Xue, Korea--Aust. Rheol. J., 14, 143 (2002).

9. R.S. Graham, A.E. Likhtman, C.B. McLeish, and S.T. Milner, J. Rheol., 47, 1171 (2003).

10. G. Ianniruberto and G. Marrucci, J. Non-Newtonian Fluid Mech., 102, 383 (2002).

11. G. Ianniruberto and G. Marrucci, J. Rheol., 45, 1305 (2001).

12. P. Wapperom, R. Keunings, and G. Ianniruberto, J. Rheol., 47, 247 (2002).

13. J. Remmelgas, G. Harrison, and L. G. Leal, J. Non-Newtonian Fluid Mech., 80, 115 (1999).

14. G. Marrucci and N. Grizzuti, Gaz. Chim. Ital., 118, 179 (1988).

15. W.L. Olbricht, J.M. Rallison, and L.G. Leal, J. Non-Newtonian Fluid Mech., 10, 291 (1982).

16. A. Cohen, Rheol. Acta, 30, 270 (1991).

17. R. Aris, Vectors, Tensors, and the Basic Equations of Fluid Mechanics, Dover, New York (1989).

18. openFoam, www.opencfd.co.uk.

19. C. Hirsch, Numerical Computation of Internal and External Flows, Wiley, New York (1991).

20. J. Remmelgas, P. Singh, and L.G. Leal, J. Non-Newtonian Fluid Mech., 88, 31 (1999).

21. H.K. Versteeg and W. Malalasekera, An Introduction to Computational Fluid Dynamics, Longman, England (1995).

22. R.I. Issa, J. Comp. Phys., 62, 40 (1986).

23. D.S. Jang, R. Jetli, and S. Acharya, Numer. Heat Transfer, 19, 209 (1986).

24. R. Ahmed, R.F. Liang, and M.R. Mackley, J. Non-Newtonian Fluid Mech., 59, 129 (1995).

25. C.D. Han, Rheology in Polymer Processing, Academic Press, New York (1976).

26. D.W. Mead, D. Yavich, and L.G. Leal, Rheol. Acta, 34, 360 (1995).

M.A. Beard, G.R. Tabor, S.J.K. Ritchie, K.E. Evans

School of Engineering, Computer Science and Mathematics, University of Exeter, UK

Correspondence to: M.A. Beard; e-mail: m.a.beard@ex.ac.uk

Printer friendly Cite/link Email Feedback | |

Author: | Beard, M.A.; Tabor, G.R.; Ritchie, S.J.K.; Evans, K.E. |
---|---|

Publication: | Polymer Engineering and Science |

Geographic Code: | 1USA |

Date: | Oct 1, 2007 |

Words: | 4310 |

Previous Article: | Effects of raw fiber materials, fiber content, and coupling agent content on selected properties of polyethylene/wood fiber composites. |

Next Article: | Toughening of epoxy resin with amine functional aniline acetaldehyde condensate. |

Topics: |