A quasi-1D model for fast contraction flows of polymer melts

J. Non-Newtonian Fluid Mech. 118 (2004) 157–173

A quasi-1D model for fast contraction flows of polymer melts Jae-Hyeuk Jeong, Arkady I. Leonov∗ Department of Polymer Engineering, The University of Akron, Akron, OH 44325-0301, USA Received 18 July 2003; received in revised form 10 February 2004

Abstract This paper develops an isothermal quasi-1D model for fast (high Deborah number) contraction flows of polymers from a reservoir to a die of circular or rectangular cross-sections. Two components of composite flow models, inhomogeneous elongation and a modified unsteady shearing, employed in different regions of the full flow domain are connected in succession with the use of asymptotic matching boundary conditions. The model, based on a stable and descriptive set of viscoelastic constitutive equations, demonstrates a reasonable good agreement with experimental and direct numerical data with no adjustable parameters. It employs easy PC numerical calculations. The predictions of the model are improved for contraction flows with higher values of Deborah/Weissenberg number and contraction ratio. © 2004 Elsevier B.V. All rights reserved. Keywords: Contraction flow; Deborah number (De); Jet approach; Developing flow

1. Introduction We develop a quasi-1D modeling of fast isothermal contraction flows of polymer melts in the axi-symmetric composite channel (Fig. 1), which consists of a channel with larger cross-sectional area (reservoir) connected to a channel with smaller cross-sectional area (die). The cross-sectional geometries of reservoir and die are considered either circular or rectangular with high ratios of the cross-sectional areas. These geometries are widely utilized in high-speed polymer processing operations, such as extrusion and injection molding, and in the testing of polymer fluids using capillary/slit die rheometers. This type of flow also represents a benchmark problem for the numerical studies of polymer flows at high Deborah (De) numbers. The complexity of these analyses, especially for high De number flows, originates from non-linear character of viscoelastic constitutive equations (CEs) and from the effect of sudden change in the geometry of the composite channels. In this paper, we will not distinguish between Deborah (De) and Weissenberg (We) numbers, and define the Deborah number as: De = θ¯ v¯ / lc . Here v¯ and lc are the average velocity and a characteristic cross-sectional dimension in the ∗ Corresponding author. Tel.: +1-330-972-5138; fax: +1-330-258-2339. E-mail addresses: [email protected] (J.-H. Jeong), [email protected] (A.I. Leonov).

0377-0257/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2004.03.005

die, respectively, and θ¯ is the relaxation time averaged over the linear relaxation spectrum. It is convenient to subdivide schematically the entire complex flow domain in five characteristic geometric regions: (i) the reservoir pre-entrance, (ii) reservoir entrance, (iii) the die entrance, (iv) the die land region, and (v) the die exit region. Additionally, there is the region of free surface flow of swelling extrudate, which is not analyzed in this paper. We assume that a fully-developed simple shearing flow of Poiseuille type occurs in the region (i) located far away upstream from the die entrance. In the reservoir region (ii), which is close to the die entrance, the previous Poiseuille flow dramatically changes to a converging type of flow suited for entering the polymer stream in the die. In this region, many experiments for polymer melts and all direct numerical simulations observe a secondary flow (vortices) at corners of the reservoir entrance. In the developing flow in the die entrance region (iii), additional restructuring happens from the upstream converging flow back to a simple shearing flow downstream in the die. The developing flow is accompanied by the release/dissipation of the stored elastic energy accumulated in the previous contraction flow. If the die is long enough, the die land region (iv) may exist as a continuation of the developing flow region. Here the developed shearing flow of the Poiseuille type occurs whose characteristic is independent of longitudinal coordinate of the die. At the very end of the die exit region (v), the die flow restructures once again to the free-surface extrudate flow. Although this

158

J.-H. Jeong, A.I. Leonov / J. Non-Newtonian Fluid Mech. 118 (2004) 157–173

Pressure

Pressure L

2RR

lb la

2LR

*T1 T

*

Lz lc

* T2

Lz 2RD

2LD (A)

(B)

Fig. 1. Sketches of the slit (A) and circular (B) rheometers: 1, reservoir; 2, die; T*, the pressure transducer locations in the slit/circular rheometers.

type of re-structuring is dynamically important, it happens in a relatively short region. Therefore, flow in this region can be approximately treated by the previous (developing or developed) flow continuing up to the die end with a sudden change in the flow type at the die end cross-section. Many experimental and computational works have been performed for the flows of polymer melts and solutions in the flow regions (ii)–(iv). Flow visualization experiments, such as laser Doppler velocimetry (LDV) [1,2], particle tracer methods [3,4], laser speckle velocimetry [5,6], and birefringence methods [7], have been commonly employed to analyze the development of flow patterns near the die entrance. The contraction polymer flow also attracted attention as a benchmark problem in computational viscoelastic fluid dynamics. This is because this flow involves geometrical complexity along with a non-trivial set of equations describing the momentum balance, continuity and non-linear viscoelasticity of polymeric liquids. After an extensive review of current finite element methods (FEM), a mixed FEM was proposed with the demonstration of advantages in simulating the planar contraction problem [8]. In the review of the pre-entrance/entrance contraction flow, White et al. [9] pointed out that the choice of CEs also influences the numerical stability and evolutionary nature of CEs. Some numerical studies [10,11] of fast viscoelastic flows near sharp entrance corners revealed computational divergence with mesh refinement using the upper convected Maxwell model (UCM). This is partially attributed to the occurrence of corner singularities. Notable research efforts [12,13] with different CEs and several benchmark problems, including the contraction flow, demonstrated that the CE of “Leonov type” is the most stable and the corner singularities do not affect computations for this CE up to the value of Deborah number approaching 700 [13]. Recent publications of direct computer simulations [14,15] of contraction flows for LDPE, which used simple integral, modified models, present some new data. These simulations, however, did not present detail computational data and comparisons with

experimental data, needed for comparison of our computations with both direct simulations and experimental data. Along with the instability problems, the authors of Ref. [9] emphasized the importance of extensional flow in entry flow. They revealed that the common failure of numerical prediction of entry flow comes from inabilities of most CEs to describe well both the extensional and shear flows. Cogswell [16] first realized the dominancy of extensional component in fast entrance flows of polymer melts. This idea has been discussed many times and finally accepted in many experimental papers cited above. Several authors attempted to develop a simplified empirical flow model for the entrance converging flow, using the Cogswell idea. Binding [17] employed empirical equations for simple shearing and simple elongation. Carrot et al. [18] and Guillet et al. [19] involved employed a specific CE of simple integral type in computations. They calculated the stress–strain rate relations in the “Protean” (stream function-longitudinal) coordinate system known well in the theory of convective diffusion [20]. In their formulation, they also employed an arbitrarily chosen “minimum principle” with an empirical distribution of the longitudinal velocity. In spite of the model arbitrariness and unstable character of original simple integral CE, the authors reported successful agreement between the model calculations and experimental data [18,19]. The present work develops a model for high De number contraction flow with abrupt change in entrance geometry. It starts with a set of general CEs and simplifies the flow complexity using two different 1D flow components in different geometrical regions. In the pre-reservoir region (i), the flow is modeled as developed simple shearing flow. Then, following the Cogswell idea [16], the main converging stream of polymer liquid in the reservoir region (ii) is modeled as 1D inhomogeneous extensional flow, similar to that employed in the fiber spinning. The viscoelastic circulatory flow in the reservoir near the die entrance was also taken into account, because it may affect the main extensional flow in the region (ii). In the region (iii) where the entrance die flow is developing, we use once again the CEs for 1D unsteady simple shearing with evolution of entrance elastic strains, unusual in simple shearing, which originate from the reservoir jet flow. The time derivative is simplistically treated in this flow model as space convective derivative as it has been used before in many papers in Newtonian fluid dynamics (e.g. see Refs. [21,22]). Additionally, more detailed momentum balance equation is used in the region (iii) to stabilize numerical calculations. Finally, the characteristics of developed flow in the region (iv) are calculated as limiting regime of flow in very long tubes. The transitions between flow types between different flow regions, except for those in (iii) and (iv), are modeled by specific matching boundary conditions. Unlike many direct [8,10–13] and very simplified [17–19] flow simulations, the present flow model involves the detailed description of the polymer melt flows in simple shearing and simple elongation, keeping these descriptions with

J.-H. Jeong, A.I. Leonov / J. Non-Newtonian Fluid Mech. 118 (2004) 157–173

all the rheological complexity (see also the brief discussion in the next section). Nevertheless, the model treats the transient flows in a very rough and simplistic manner based mostly on intuitive assumptions and experimental findings. Although some of these simplifications could be asymptotically justified, no attempt of that type is demonstrated in the present paper. Therefore, in the following, we will use these “modeling assumptions” in a very declarative manner, judging their validity only by comparison of model calculations with the results of direct computations and/or experimental data. We show that our flow model, designed to save computational cost, demonstrates a reasonably good agreement with existing available experimental results and direct numerical simulations. We should also mention that our model calculations demonstrate better comparison with experimental data and/or direct calculations for higher De number flows and/or for higher contraction ratios. This is perhaps because in flows with higher values of these parameters, the transitions between the above basic flow sub-regions are getting sharper. It is also worth mentioning that no single adjustable parameter is involved in the present flow modeling.

2. Constitutive equation In the following, we use a multi-mode viscoelastic CE of differential type. It has been employed [23,24] for detailed and consistent description of all experimentally available rheological data for five polymeric melts, while satisfying all the stability constraints [25]. For each relaxation mode with relaxation time θ k and elastic modulus Gk , the CE operates with a modal elastic Finger tensor c whose evolution for incompressible case is k described as
 ∇ 2θk c + b c2 + 13 c (I2k − I1k ) − δ = 0, k

k

k

∇

c ≡ c˙ − c · ∇v − (∇v)T · c,

I2k =

trc−1 , k

I1k = trc , k

I3k = detc = 1; b = b(I1 , I2 ).

(1)

k

Here ∇v is the velocity gradient tensor, δ is the unit tensor, and b is a positive scaling factor for the relaxation time, which is generally a function of invariants, I1k and I2k with parameters independent of the mode number. This function has been specified for various polymeric melts and elastomers [23]. The extra stress tensor σ and the total stress σ e are defined as σ = −pδ + σ ; e   ∂Wk σ =2 σ = c · , e k k ∂c k

k

k

3Gk (T) Wk ≡ Gk (T)wk (I1k ) = 2(n + 1)



I1k 3

n+1

 −1 .

(2)

159

Table 1 Dimensionless material parameters defined in Eqs. (1) and (2) Polymer melts

Dimensionless parameters

Polyisobutylene (Vistanex) Isayev and Upadhyay [7] Polyisobutylene (LM-MH) LDPE (LDPE Lupolen 1810 H from BASF Ludwigshafen; charge number 912 133 036) [28]

n = 0.1 and b(I1 , I2 ; m) = 1 n = 0.1 and b(I1 , I2 ; m) = 1 n = 0.03 and b(I1 , I2 ; β, ν)  = exp[−β (I − 3)] sinh[ν(I − 3)] , where + ν(I − 3) + 1 β = 0.35, ν = 0.002

Here p is the isotropic pressure and σ is the extra stress k tensor for kth relaxation mode. Here we have simplified the general form of modal elastic potential Wk proposed in [23]. The potential in (2) employs the Hookean modulus Gk (T) that weakly depends on temperature, and the non-dimensional function, wk , which is characterized by only one, mode-independent numerical parameter n. The general choice of the function b(I1 , I2 ) considered in papers [23,24] and parameter n should be specified depending on rheological behavior of a polymer. The formulations and values of these non-linear terms for some polymer melts are presented in Table 1. The additional rheological parameters to be specified are the linear relaxation spectrum (Gk , θ k }. It should be mentioned that in the typical case n > 0 in (2), the CEs (1)–(2) are not reducible to the equations written in terms of only extra stress tensor. This is possible only in the neo-Hookean case when n = 0. However, this case violates the stability constraints and worsens description of high Deborah number properties for molten polymers, especially for elongation flows. Nevertheless only the case n = 0 has been employed before in all the numerical simulations of flows of molten polymers. Additionally, only few CEs for polymer melts exist that demonstrate the good descriptive properties still satisfying all the stability constraints [23,24]. Therefore, in the following text we will use only CEs (1)–(2) keeping in mind that the model might also be used for these other few stable and descriptive viscoelastic CEs.

3. Modeling the entrance contraction flow Our schematic interpretation of complex flow, outlined in Section 1, assumes that two basic regions of polymer flow exist in the reservoir. When the reservoir is long enough, there is a developed steady shear flow of polymer melts in the region (i). We call the flow in this region the far field entrance flow. Additionally, there is the flow region (ii) near the die entrance where the main stream of flow converges to the die. This flow sometimes is accompanied by weak entrance vortices at the flow periphery, near the entrance corner. We call the flow in this region the near field entrance flow. Let x1 be the axial coordinate directed downstream along the symmetry axis with the origin (x1 = 0) at the die

160

J.-H. Jeong, A.I. Leonov / J. Non-Newtonian Fluid Mech. 118 (2004) 157–173

orifice. Then the far and near field flows are located in the respective axial domains (x1 < −l) for the region (i) and (0 > x1 > −l) for the region (ii); the coordinate (x1 = −l), being the parameter searched for. Although the change of the flow characteristics from the far field to the near field entrance flow needs a (presumably small) transitional axial length δ (l), we asymptotically will use the value x1 = −l as the boundary in this model where sudden change in flow type happens. We initially approximate the contraction flow in the near field entrance region (ii) as an inhomogeneous extensional (or “jet”) flow that commonly occurs either in polymer fiber spinning or in polymer sheet processing operations. We assume that such an approach is approximately valid when the Deborah number of the flow is high enough. In order to compare our calculations with the literature birefringence data obtained in moderate De flows, we then develop a more detailed approach involving into the analysis the secondary shear flow, and include the viscoelastic secondary flow in dynamics of the jet. We further demonstrate that for moderate De flows, the contribution of the secondary flow might be important for initial development of the jet (see details of these calculations below) where the local Deborah numbers are too low. Two types of reservoir geometries are analyzed in the paper: the slit type reservoir (Fig. 1A) with cross-sectional sizes 2LR (thickness) and L (width) where LR  L, and the reservoir of circular geometry with the reservoir radius RR . To avoid unnecessary duplications and simplify the notations, a simultaneous analysis of flows in both types of geometries is presented below with asterisked formulae attributed to the circular geometry. The same notations {x1 , x2 , x3 } are used below for the Cartesian and cylindrical coordinate systems to describe the respective flows in the slit (−LR < x2 < LR ) and circular (0 < x2 < RR ) geometries. The index 1 stands for axial (z); 2, for transversal (y or r∗ ); and 3, for neutral (x or ϕ∗ ) directions. 3.1. Flow in the far field entrance region (i): {x1 < −l, LR > x2 > −LR or 0 < x2 < RR } According to our assumption, the steady simple shear flow of Poiseuille type happens in this region. In this flow, the matrices of velocity gradient, ∇v, the extra stress, σ , and e the elastic strain, c , tensors are presented as 

k

  0 0 0 σ11,e    ∇v = γ˙  1 0 0  , σ =  σ12,e e 0 0 0 0   c11,k c12,k 0   c =  c12,k c22,k 0  . k 0 0 1

σ12,e

0

σ22,e 0

σ33,e



 0 ,

(3)

Substituting (3) into (1) written for steady state yields the known solution [23,24]

√

√ 2ζk 2 , , c22,k = √ c11,k = √ ζk + 1 ζk + 1  ζk 2 − 1 c12,k = , c33,k = 1, ζk + 1    2θk γ˙ 2 ˙ = 1+ ; ζk (γ) b    I1,k n dv1 Gk cij,k ; γ˙ = σij,e = ; 3 dx2 k  I1,k = I2,k = c11,k + c22,k + 1 = 2(ζk + 1) + 1.

(4)

Here γ˙ is the shear rate and v1 is the only non-zero velocity component in this flow. The momentum balance equations result in σ12,e (x2 ) =

dP x2 ; dx1

σ12,e (x2 )∗ =

(5)

1 dP x2 . 2 dx1

(5*)

Here P is the pressure and σ 12,e is the shear stress. The expression for flow rate Q is  LR  LR Q = 2L v1 dx2 = −2L x2 γ˙ dx2 ; (6) 0

Q∗ = 2π

0



RR 0



RR

v1 x2 dx2 = −π 0

x22 γ˙ dx2 .

(6*)

Here the constant value of Q is considered below as given. When the Q value is specified, the shear rate, γ(x ˙ 2 ), can be found from Eqs. (4)–(6). It enables us to compute the complete profiles of rheological variables for the steady shear flow up to the distance l from the entrance. In our simple numerical procedures, we employed the Newton–Raphson method to find the gap-wise shear rate distribution for any given value of flow rate Q, while satisfying Eqs. (4)–(6). The numerical integration of Eqs. (6) and (6*) was performed using the trapezoidal method. We should note that the location l of the borderline between the far field and the near field reservoir entrance flows has not been found yet. 3.2. Flow in the near field entrance region (ii) {0 > x1 > −l, LR > x2 > −LR or 0 < x2 < RR } 3.2.1. Jet flow approach We assume that when the flow Deborah number of the flow is high enough, the contraction flow in the near field entrance region (ii) can be modeled as an inhomogeneous extensional “jet” flow with flat profiles of longitudinal velocity and extensional stress. This assumption is directly employed from the Cogswell idea [16] and the flow visualization experiments [7,13,19], where the measured core stress distribution for high De number flows becomes more jet-like near the die. This simplified rough modeling ignores vortices at

J.-H. Jeong, A.I. Leonov / J. Non-Newtonian Fluid Mech. 118 (2004) 157–173

the entrance corners and neglects the shear effects near the reservoir walls. We now consider the jet flow in the slit and circular reservoir geometries. In the slit case the jet flow can be approximated as an inhomogeneous planar elongation flow. In the circular case, we use the inhomogeneous simple elongation approximation. These approximations are similar to the fiber spinning/sheet processing cases. Then the corresponding expressions for the matrices of elastic strain tensors, c , vek locity gradient tensor, ∇v and the extra stress tensor, σ , in e the jet are   2   λk 0 0 1 0 0     c =  0 λ−2 ∇v = ε˙  0 −1 0  , 0, k k 0 0 0 0 0 1   0 σ11 0   σ =  0 σ22 0  ; (7) e

0 

0

λ2k

0

 c∗ =  0 k

λ−1 k

0 

0

σ33 

0

 0 ,

λ−1 k

1 0  ∗ ∇v = ε˙  0 −1/2 0 0  0 σ11 0  ∗ σ =  0 σ22 0 e

0

0

 0  0 , −1/2   .

(7*)

σ33

Here ε˙ is the elongation rate. Inserting (7) into (1) and employing the “convective approximation”, d/dt ≈ v¯ 1 d/dx1 , commonly used for the description of these inhomogeneous elongation flows, yields   b 1 1 dλk 2 λk − 2 = ε˙ ; (8) v¯ 1 + λk dx1 4θk λk    1 b 1 1 dλk 1+ = ε˙ . (8*) λ2k − v¯ 1 + λk dx1 6θ λk λk Here v¯ 1 is the axial velocity averaged over the jet cross-section, and ε˙ = d¯v1 /dx1 is the elongation rate where the extensional rate (or the cross-sectional area of flow) is a function of x1 . The expression for the elongation stress (a total axial stress in jet cross-sections), σ ext , in the multi-mode case is presented as   I1,k n σext = G (λ2k − λ−2 I1k = 1 + λ2k + λ−2 k ), k ; 3 k

∗ = σext

 k

 G

I1,k 3

n

(9) (λ2k − λ−1 k ),

I1k = λ2k + 2λ−1 k . (9*)

161

Due to the incompressibility condition held in the jet flow v¯ 1 (x1 ) =

Q , A(x1 )

ε˙ ≡

d¯v1 Q dA =− 2 . dx1 A dx1

(10)

Here A(x1 ) is the cross-sectional area of jet, the flow rate Q being the same as in the region (i). With Eq. (10), the evolution Eqs. (8, 8*) for each relaxation modes take the form   Q dλk b Q dA 1 λ2k − 2 = − 2 + ; (11) Aλk dx1 4θk A dx1 λk Q dλk b + Aλk dx1 6θk

 λ2k

1 − λk



1 1+ λk

 =−

Q dA . A2 dx1 (11*)

The set of Eqs. (11, 11*) for each mode is not closed yet for unknowns A and λk . To make the closure, we involve another assumption that the force acting on the jet is constant: Fext = σext × A = const. This is exactly the same as for the real inhomogeneous elongation flow. This condition, neglecting contribution of near-wall shear stresses in the longitudinal force, is expected to work well enough for high De flows where the jet flow dominates. Using this condition, Eqs. (9) and (9*) is written in the final, form   I1,k n Fext ≡ Aσext = A G (λ2k − λ−2 k ) = const, 3 I1k = ∗ Fext

1 + λ2k

≡

k −2 + λk ;

∗ A∗ σext

=A



(12)  Gk

k ∗ I1k = λ2k +

I1k 3

n 

λ2k

1 − λk

2 . λk

 = const, (12*)

Eqs. (11) and (12) represent a closed set. Appropriate boundary conditions for these equations are (see Fig. 1) x1 = −l :

A = Al = 2LLR ,

A = Al = πR2R , x1 = 0 :

λk = λlk .

A = A0 = 2LLD ,

A = A0 = πR2D ,

λk = λlk ;

λk = λ0k .

(13) (13*)

λk = λ0k ;

(14) (14*)

Here the known geometric parameters L and LD (RD ) are, respectively, the width of slit die the half thickness (radius) of slit planar (circular) die. Note that as soon as the values of elongation elastic stretches λlk are known, the steady non-linear boundary value problem is completely closed. To find the values λlk we will use an additional matching condition (A4) at the unknown boundary x1 = −l, derived in Appendix A.

162

J.-H. Jeong, A.I. Leonov / J. Non-Newtonian Fluid Mech. 118 (2004) 157–173

3.2.2. Jet flow with drag from secondary viscoelastic flows We now present a simplified analysis of the secondary shear flow effect in the modeling of the “jet flow”. When the value of De number for the initial jet is high enough, the secondary flow effect is negligible in the near field entrance region, because it would be suppressed by the dominant jet flow. However, the birefringence data [7] clearly showed that for the flows with moderate De numbers, the initial extensional jet would not be dominant over the secondary shear, and therefore, the shear flow may affect the jet development. We approximately model the complex flow in the near field entrance region as a combination of inhomogeneous elongational (jet), core flow near the axis of symmetry, and a near-wall shearing (secondary) flow outside the jet region. The interaction between these two different types of flow results in the drag imposed on the extensional core flow from the shearing one. It means that the secondary flow affects the normal force and the cross-sectional area of jet, changing the jet dynamics over the region (ii). This modeling could be visualized as an extension of a rubber cylinder (or sheet) moved through another liquid in a confined cylindrical container, with neglecting the shear stress (as compared to the extensional one) in the extended cylinder. It is difficult for this type of modeling to distinguish between the basic jet (elongation) flow and the near-wall secondary (shear) flow. To avoid this uncertainty and keep computational simplicity, we simply define the jet as the near-centerline extensional stream that completely contributes to the flow rate. It means that the secondary flow is circulatory and confined in the close-wedge-shaped domain near the corner. The circulatory flow is overwhelmingly observed in numerous experimental data for high De number polymer flows and its existence was well documented in all direct computations. We will model the secondary shear flow using the lubrication approximation [21], which is commonly applied to a wedge-like shearing between two well-defined solid boundaries. In our case the difficulty is that the jet profile is unknown and searched for. Within the lubrication approximation, the momentum balance equations for the secondary flows are written in a simplified form σ12,e (x2 ) =

x2 dP + c1 ; dx1

∗ (x2 ) = σ12,e

 x  dP c1 2 + . 2 dx1 x2

(15) (15*)

Here the pressure P and unknown parameter c1 are commonly considered as independent of x2 , and σ 12,e is the shear stress. In the circulatory flow, the flow rate is equal to zero  LR Qsec = −2L x2 γ˙ dx2 = 0; (16) LEx

Q∗sec = −π



RR REx

x22 γ˙ dx2 = 0.

(16*)

Here LEx (REx ) is the half-thickness/radius of jet and LR (RR ) is the half-thickness/radius of slit/circular reservoir, and γ˙ ≡ ∂v1,sec /∂x2 is the shear rate in the secondary flow, which can be expressed via shear stress using the formulae (4). The non-slip boundary conditions at the reservoir wall and the continuity of velocity at the jet surface are x2 = LR or RR : x2 = LEx or REx :

v1,sec = 0; v1,sec = v¯ 1,jet =

(17) Q . A(x1 )

(18)

Consider now the effect of drag forces on the jet flow. The local force balance between the extensional jet force and shear drag acting in the x1 direction is d (σext LEx ) = −σ12,surf (x1 ), dx1 d ∗ (σ ∗ πR2 ) = −2πREx σ12,surf (x1 ). dx1 ext Ex

(19) (19*)

Here σ ext is the extensional stress acting on the jet and σ 12,surf is the shear stress at the jet surface due to circulatory drag flow. The finite difference iterative algorithm used for calculations in this case is more complicated than that for the pure jet approach. The brief description of the algorithm on the example of plain flow is as follows. Assume that at a certain value of x1 the jet profile, LEx (x1 ), and the extensional stress, σ ext (x1 ), are known, then the jet velocity v¯ 1,jet (x1 ) = Q/A(x1 ) is known too. Now the problem of viscoelastic shearing of circulatory flow, which is presented by (4), (15)–(18), is well defined. It has been elaborated long ago in extrusion theory that a monotonously increasing function σ12 (γ) ˙ can be found from the given “flow curve”. Then the inverse function, γ˙ = ∂v1,sec /∂x2 = −ψ(σ12 ), is well ˙ is defined in defined too. In our case the function, σ12,e (γ), Eq. (4). Taking into account Eqs. (15)–(18), the gap-wise distribution of the axial velocity and the condition for determining two unknown parameters, dP/dx1 and c1 , are   LR  dP v1,sec = ψ x2 + c1 dx2 , dx1 x2   LR  dP ψ x2 + c1 dx2 , v¯ 1,jet = dx1 LEx    LR dP x2 ψ x2 + c1 dx2 = 0. (20) dx1 LEx Using a numerical procedure outlined below and the formulae (20), the problem of circulatory viscoelastic flow can be solved for any value x1 including the determination of the drag shear stress σ 12,surf (x1 ). Then the numerical calculations of Eqs. (11) and (12) with variable extensional force Fext are performed using a small increment )x1 to obtain the values of the jet profile and extensional force at the level x1 + ∆x1 . It is evident that these calculations show a positive increase in the extensional force, )Fext , caused by the

J.-H. Jeong, A.I. Leonov / J. Non-Newtonian Fluid Mech. 118 (2004) 157–173

drag of the over surrounding melt (recall that the surrounding melt acts like a frictional medium to the extended jet). The above iterative procedure started at the upstream corner x1 = −l, where LEx → LR (REx → RR ). Unfortunately, at this upstream corner, a singularity in stress behavior occurs for our modeling of secondary flow. Analysis in Appendix B shows that this singularity is integrable, and also gives an analytical description of our calculations near the singular point. It should also be noted that the procedure outlined for taking into account the near-wall circulatory flow does not guarantee the continuous value of stress invariant at the jet boundary. Moreover, this value as calculated using the flow model is always discontinuous. However, when the shear stress at the jet boundary exceeds the value of normal stress in the jet, we see it as a clear indication that the main assumption of the jet dominancy at a local place in the polymer entrance flow, is invalid. The numerical procedure for the unknown parameters {γ(x ˙ 2 ), c1 } is performed using a common matrix solver subroutine with satisfying Eqs. (15, 15*) and (16, 16*) and boundary conditions (17) and (18). These calculations completely solve the secondary flow problem including the shear profile at a given axial position and the drag at the jet surface. After that, the jet extensional force for the next step, x1 + )x1 , can be computed from Eqs. (19) and (19*). It is easy to understand that the circulatory shear flow model will decrease the jet length l1 . This is because the drag increases the extensional jet force, which in turn makes the jet thinner. Therefore, when the initial (reservoir) and final (die) areas are fixed the transition from the initial to final jet thickness is getting shorter. 3.2.3. Matching condition at the boundary between the two entrance flow regions The formulation of the jet flow model has not yet been completed because the initial elastic stretches λlk (k = 1, 2, . . . ) in Eqs. (13) and (13*) are still unknown. To determine them we employ at the unknown boundary, x1 = −l, an energetic matching condition (A.4) derived in Appendix A for each kth non-linear Maxwellian mode. This energetic condition, proposed first in Ref. [26], gives a way to setup the asymptotic matching conditions for the flow parameters between the different characteristic regions. Even if those parameters have discontinuities in the values at an effective “interface”, the evaluation of the global flow properties are still well-posed. In real flows, these parameters are changed continuously from one asymptotic flow region to another, but these calculations are the subject of direct numerical simulations. The matching condition (A4) takes the form  1 jet jet x1 = −l : vsh W sh dA = v1 Wk . (21) AR Ω 1 k jet

Here Wksh and Wk are elastic potentials for each kth non-linear Maxwell mode in the far field shearing and near

163

jet

field jet flows, vsh 1 and v1 are the axial velocities of shear and jet flow, and AR is the reservoir cross-sectional area. Using formulae (4) and (9) for simple shearing and elongation flows, Eq. (21) reduces for any kth relaxation mode to the form x1 = −l :

L Q



LR

n+1 vsh dx2 1 (c11,k + c22,k + 1)

0 l 2 = [(λk ) + (λlk )−2

x1 = −l :

2π Q



RR

+ 1]n+1 ,

x2 vsh 1 (c11,k 0 = [(λlk )2 + 2(λlk )−1 ]n+1 , k

k = 1, 2, . . . ;

(22)

+ c22,k + 1)n+1 dx2 = 1, 2, . . . .

(22*)

At x1 = −l, the jet cross-sectional area coincides with that for the reservoir. Eqs. (22) and (22*) allow determination of the initial jet elastic stretches λlk from the known elastic strain tensor profile for the far field shear flow cij ,k , which have been expressed in Eqs. (4)–(6). Now inserting the values of λlk in Eq. (9) with the reservoir cross-sectional area, the total longitudinal jet force can be easily calculated at x1 = −l. As the result, the evolution in the jet area depending on axial position, A(x1 ), can be expressed via σ ext through the calculated value of Fext . Eqs. (11) and (12) along with conditions (13), (14), and (22) represent a simple model for steady entrance flow calculations. Mathematically, this problem is reduced to a well-defined boundary value problem with the parameter l being the eigenvalue. The above model for reservoir flow also provides appropriate boundary conditions for the calculations of developing flow in the die. Our numerical procedure was performed with trapezoidal integration and the root finding routine for Eq. (22) and (22*). The Newton–Raphson methods is applied for solving Eqs. (11) and (12) until the boundary condition (14) is satisfied. When the Deborah number and the contraction ratio of the flow are not high enough to satisfy our prior assumptions for the jet flow model, the other set of force conditions for the jet model with drag from secondary shear flow may work better. The method of calculation using this more complicated model has been demonstrated in the previous section and in Appendix B. One of the goals of this modeling is calculating the values of elastic strains λ0k at the die entrance, which are used as the entrance boundary conditions in the following calculations of developing die flow. These entrance boundary conditions contain memory of the previous reservoir flow, and highly affect the developing die flow. Although the described jet approach oversimplifies the reservoir near field flow in the vicinity x1 = −l, it is getting more realistic when the flow is closer to the entrance. It is also remarkable that this approach has no fitting parameters.

164

J.-H. Jeong, A.I. Leonov / J. Non-Newtonian Fluid Mech. 118 (2004) 157–173

will be guaranteed if the continuity conditions for each component of elastic strain tensor are satisfied. For the plain situation this demand fits the basic condition of existence of 0 = 1. However, for circular geometry, simple shearing, c33,k 0 0 when c33,k = 1/λk , the existence of simple shearing conditions is not satisfied. The calculations for shearing with the 0 initial case c33,k = 1, can also be performed using the CEs (1)–(2), but they are unnecessarily complex. Therefore, we avoid this complexity, because the model is designed to work for large values of Deborah numbers, and consequently at large values of λ0k in the dominant relaxation modes. To avoid excessive of complexity we will use simplified boundary conditions at the die entrance, x1 = 0 following once again (A4) from Appendix A for matching the extensional jet flow in the reservoir and developing shear flow in

4. Modeling developing flow in the die: {Lz > x1 > 0, LD > x2 > −LD (0 < x2 < RD )}

We now model the developing flow in the die for both the slit and circular geometries. We will treat the developing flow as a version of non-steady viscoelastic Poiseuille flow where the time differentiating operator is replaced by to the axial convective space operator as: d/dt ≈ v¯ 1 ∂/∂x1 . Here v¯ 1 = Q/A0 is the die average longitudinal velocity, and A0 is the cross-sectional die area equal to 2L × LD or π(RD )2 for the slit and circular die geometries, respectively. This approach is similar to one that was successfully used for the calculations of developing viscous flows in long tubes [21]. The evolution equations for each kth mode (the index k is omitted) for elastic strain tensors c shown in Eq. (3) for k both the die geometries are  ∂c12 b ∂v1   v¯ + c12 (c11 + c22 ) = c22 (“12” component),   1 ∂x1 2θ ∂x2 ∂c22 b 2 2 v¯ 1 + (c12 + c22 − 1) = 0 (“22” component),   ∂x 2θ  1  2 c11 c22 − c12 = 1 (incompressibility condition). The shear and longitudinal normal stress components for the extra stress tensor are  n  n   I I σ12 = Gk c12,k ; σ11 = Gk c11,k . 3 3

(23)

the die for each mode x1 = 0 :

v01 = v¯ 1 (v2 = 0);

(24)

0 = 0, c12

0 c33 = 1,

Along with (23) and (24), we further use only the longitudinal component of the momentum balance equation ∂p ∂σ11 ∂σ12 = + ; (25) ∂x1 ∂x1 ∂x2

x1 = 0 :

v01 = v¯ 1 (v2 = 0);

0 c12 = 0,

0 c33 = 1,

k

k

∂σ11 1 ∂(x2 σ12 ) ∂p = + , ∂x1 ∂x1 x2 ∂x2

(25*)

and the continuity equation ∂v2 ∂v1 + = 0; ∂x1 ∂x2

(26)

1 ∂(x2 v2 ) ∂v1 + = 0. ∂x1 x2 ∂x2

(26*)

Eqs. (26) and (26*) result in   LD v1 dx2 = −2L Q = 2L 

0

Q = 2π

LD

x2 γ˙ dx2 = const;

(27)

0

RD



RD

x2 v1 dx2 = −π

0

0

x22 γ˙

dx2 = const.

x2 = ±LD : x2 = RD :

v1 = v2 = 0; v1 = v2 = 0.

(28) (28*)

At the die entrance, the continuous boundary conditions should be taken for the stress tensor in each mode. These

1 0 c22

= (λ0 )2 ;

(29)

0 0 c11 + c22 + 1 = (λ0 )2 +

2 . (29*) λ0

Formula (29*) for the circular case means that there is a small restructuring region near the die entrance cross-section. Asymptotic calculations resulted in the fact that this small restructuring zone is λ0k times less than the effective lengths of developing flow region. Using (29*) in the incompressibility condition and the 0 = 1/c0 , exactly the same as in (29), we can relation c22 11 rewrite the boundary conditions (29*) for flow at the entrance of cylindrical die as x1 = 0 :

v01 = v¯ 1 (v2 = 0);

0 c12 = 0,

0 c33 = 1,

(27*)

The non-slip boundary conditions for the components of velocity are

0 c11 =

J(λ0 ) =

0 c11 =

  1 2 (λ0 )2 + 0 . 2 λ

1 0 c22

= J(λ0 ) +



J 2 (λ0 ) − 1; (29**)

The problem formulated by Eqs. (23)–(27) along with boundary conditions (28, 28*) and (29, 29**) treats the viscoelastic developing flow in the entrance region as a transitional flow between the elongation jet flow in the near field entrance and the developed shear Poiseuille flow, asymptotically achieved as x1 → ∞. As discussed above, it means

J.-H. Jeong, A.I. Leonov / J. Non-Newtonian Fluid Mech. 118 (2004) 157–173

that the developing flow in the die reflects the memory of the reservoir elongation entrance flow. Under given flow rate condition, Q = const, the total pressure drop )Ptot between two pressure transducers, T1 and T2 shown in Fig. 1, is now calculated as )Ptot = )P1 + )P2 .

Table 2 Material parameters and experimental conditions for the contraction flows with geometry shown in Fig. 1A Polymer melts

(30)

Here )P1 is the pressure difference due to the steady shear flow in the reservoir without converging flow region and )P2 is due to the shear flow in the slit die. Those quantities can be easily obtained for )P1 from Eqs. (5) and (5*) and for )P2 from Eqs. (25) and (25*). The calculations of initial elastic stretches λ0k , which are established from the analysis of near field entrance flow, represent the first necessary step in our computations of developing die flow. The next numerical step employs a finite difference numerical procedure with using a matrix solver for finding the shear rate in each cross-sectional grid. From the boundary conditions (28) and (29), the initial values of 0 , c0 } at x = 0 are considered as known, variables {v01 , c12 1 22 and the values {c12 , c22 , )P} at the first step )x1 are found using the initial values. The same procedure is applied for calculations of all variables in cross-sections along the die.

5. Results and discussion All the calculations in this paper are performed for isothermal viscoelastic flows under various constant values of flow rate based on the general CEs, Eqs. (1) and (2) with the specification of scaling relaxation factor b given in Table 1. The following necessary information is needed for comparison of model calculations with experimental data or direct numerical simulations. We firstly need the rheological characterization of the polymer system, which provides us with knowledge of the discrete relaxation spectrum {Gk , θ k }. We secondly need a precise description of the geometry with a sharp entrance transition, used in these research papers. We thirdly need experiments or direct computations, presenting at least some space distributions of stress and/or velocity fields that might be compared with our model calculations. Additionally, the contraction flows of a dilute polymer solution are impossible to use for our modeling.The problem  is that the flow De number (De = θ¯ γ˙ App = k Gk θk2 / k Gk θk γ˙ App ) for these solutions is not high enough to use our modeling. Additionally, a high increase in De number also highly increases the Reynolds number, involving significant inertia effects. We reiterate once again that both the flow Deborah number and contraction ratio are the most important factors in the modeling of viscoelastic contraction flows. 5.1. Calculations and experiments for polyisobutylenes For the slit channel flow, we use the experimental data [7] for a polyisobutylene (Vistanex), where the discrete

165

Polyisobutylene (Vistanex) Isayev and Upadhyay [7]

Parameters at 300 K θ k (sec)

Gk , (KPa)

7.025

4.55

1.553

11.59

Geometry (cm)

L = 2.0, Lz = 3.5, la = 4.0, lb = 2.02, lc = 0.95 Die 1:

0.182

86.81

2LR = 0.8, 2LD = 0.4 Die 2 2LR = 0.822, 2LD = 0.121

relaxation spectrum (Maxwellian modes) has been found. For the circular channel contraction flows, the data for another type of polyisobutylene (LM-MH) were used. We obtained the Maxwell modes shown in Table 3 employing the Pade–Laplace procedure [27] from the linear dynamic data and stress relaxation experiments performed using the disc–disc fixture in RMS 800 rheometer. It should also be noted that our CE describes well the uniaxial elongation data obtained using an elongation rheometer of Meissner type. The Monsanto processability tester (MPT) instrument with a circular channel geometry and controlled by piston speed was also used in our experiments. Along with the values of material characteristics, the values of geometrical parameters noted in Fig. 1 are displayed in Tables 2 and 3. Fig. 3 illustrates the comparison between our calculations and experimental data [7] for Vistanex. Using the birefringence technique, the stress was obtained at the various axial positions in near field entrance region from the stress-optical “law” [7], Σ = )n/1.414e − 9 where the stress invariant is 2 ]1/2 . Here )n is the presented as: Σ = [(σ11 − σ22 )2 + 4σ12 birefringence and σ ij are the stress components. Employing the die 2 (see Table 2), which has a higher contraction ratio (≈7:1) than die 1, the experimental data (filled circles) were compared to the jet model (solid lines) and to the jet with drag model (dashed lines). The model of a jet with drag is based on the extensional jet for near-centerline core flow and the shear flow for vortices, with coinciding the jet and shear velocities at the jet boundary (e.g. see the schematic interpretation of the model in Fig. 2). Therefore, when the stress invariant for the core Table 3 Material parameters and experimental conditions for the high De number capillary flows with geometry shown in Fig. 1B Polymer melts

Polyisobutylene (LM-MH)

Parameters at 323 K θ k (s)

Gk (kPa)

0.006 0.046 0.204 0.815 4.091

144.72 47.38 19.63 7.517 1.656

Geometry (cm) Monsanto processability tester Reservoir radius (RR ) = 2.6162, die radius (RD ) = 0.07544, die length (Lz ) = 3.0176

166

J.-H. Jeong, A.I. Leonov / J. Non-Newtonian Fluid Mech. 118 (2004) 157–173

stress (Pa)*10-5

x1 /LD=-4.36

Far-field Reservoir Entrance Flow

1

0 0

Die Developing flow Die Developed flow

x1 /LD=-2.2

(A)

-5

x1 /LD=-3.34 stress (Pa)*10

Real flow velocity profile Approximated Jet vel. profile

stress (Pa)*10-5

Near-field Reservoir Entrance Flow

1 x2 /LR

1

1

0

0 0

0

1

-5

1

0

jet was expressed as Σ = (σ11 − σ22 ), and the stress for vortices was as Σ = 2|σ12 |, a discontinuity usually exists at the boundary due to the different characteristics of flow across the gap direction. According to the experimental data in Ref. [7], the local De numbers of the shear flows in far field entrance and die regions (De = θ¯ γ˙ App = 6.35 × 3Q/2LL2R ) are equal to 2.02 and 93.16, respectively. Because of the low values of the local De number in the far field shear flow, the first two plots in Fig. 3 show that model calculations (lines) disagree with experimental data (symbols). In this transient region, jet profiles have not been developed. However, close to the die where the local De number highly increases, the jet flow profile builds up, and consequently, the model predicts the details of flow profile well enough for the rest of the region. Therefore, we expect that the model will have better predictions when the initial far field flow has a higher local value of De. This can happen as a result of shortening the transient length δ from the far field shear to the near field jet flow. Having only few experimental points for flow in die 1, we cannot demonstrate that an additional shear drag in the near field entrance could significantly improve the stress distribution. But one can see this trend in considerable improving comparisons between our calculations and experimental data for the die 2 (Fig. 3). Fig. 4A presents the comparison between experimental data and our computations for the pressure difference between two transducers shown in Fig. 1. Here we used Eq. (30) with the given values of flow rate. According to Eq. (30), the near field reservoir region does not directly contribute in the total pressure. However, it indirectly affects

x1 /LD=-1.44 stress (Pa)*10

Fig. 2. Schematic classification of contraction flow regions.

stress (Pa)*10

(B)

-5

x1 /LD=-1.74

Jetwith non-Newtonian secondary drag flow profile

1 x 2/LR

x2 /LR

1

0 0

1 x2/LR

0

1 x2/LR

Fig. 3. Gap-wise profiles of stress invariant for PIB (Vistanex) contraction flow in the die 2 geometry at various axial locations at 27 ◦ C, with the flow rate, Q = 7.16 × 10−8 m3 /s: filled circles, experiments; solid lines, jet calculations without drag; dashed lines, jet calculations with drag.

the pressure drop in the developing die flow via the initial values of elastic strains. Therefore, as far as the model describes well enough the developing flow region with reasonable initial values at the die entrance, we expect a good agreement of our calculations with the data for the plots of pressure drop versus flow rate. Fig. 4A shows that the experimental data were described well by the model with and without account of the drag. Our computations also show that the extensional force with the extra shear drag (dashed line) slightly increases the pressure drop ()P1 ) in the reservoir region as compared to the pure jet approach (solid line). This is a result of shortening the length of extended jet. However, the calculated pressures in the die ()P2 ) with the initial conditions, obtained by both model calculations, have almost the same values because )P2  )P1 . Both model types estimate quite successfully the initial values of elastic strains for the die flow. Also, as mentioned earlier, the increase in the contraction ratio and increasing the flow rate causes the increase in extensional force ()Fext ). This effect is demonstrated in Fig. 4B where one can clearly see the difference between the jet flow profiles without (solid line) and with the drag (dashed line) from the circulatory flow.

J.-H. Jeong, A.I. Leonov / J. Non-Newtonian Fluid Mech. 118 (2004) 157–173

167

P (Pa)

1e+7

1e+6 Die 1 Die 2 calculation w/o drag calculation w/ drag 1e+5 1e-9

1e-8

1e-7

1e-6

3

(A)

Q (m /s)

1.2

jet flow height/LR

1.0 0.8 0.6 0.4 calculation w/o drag calculation w/ drag

0.2 0.0 -20

-15

(B)

-10

-5

0

x 1 /L D

Fig. 4. Contraction flow of PIB (Vistanex) at 27 ◦ C: solid lines, jet calculations without drag; dashed lines, jet calculations with drag. (A) The plots of pressure differences between two mounted transducers (T1* and T2* in Fig. 1 slit rheometer) vs. flow rate Q: (䊊) experimental data for die 1, (䊐) experimental data for the die 2; (B) non-dimensional plots of the jet profiles, calculated with and without drag, vs. axial location in the reservoir. Flow in the die 2 at flow rate Q = 7.16 × 10−8 m3 /s.

We now consider the comparison of our calculations with experimental data for another polyisobutylene (LM-MH). Fig. 5A–D demonstrates basic rheological properties for this material. The data in Fig. 5A–C were obtained using a prallel plate fixture with an RMS 800 instrument, and the data in Fig. 5D were obtained using an elongational rheometer of Meissner type. Employing the Pade–Laplace method [27], we obtained the value of parameters in discrete relaxation spectrum from the dynamic (A) and the stress relaxation (B) experiments in linear region. From the steady simple shear experiments with different temperatures (Fig. 5C), we found the activation energy using the Arrhenius equation, and confirmed the time–temperature superposition. Using the multi-mode constitutive approach (1)–(2) with these Maxwell modes, we described experimental results in Fig. 6, including the non-linear response in uniaxial elongation (Fig. 5D). In the contraction flow with given circular capillary geometry shown in Table 3 (see also Fig. 1) the contraction ratio is about 35, and for the highest flow rate 6.041 × 10−8 m3 /s, the De numbers in reservoir (De = 1.7 × 4Q/π(RR )3 ) and die (De = 1.7 × 4Q/π(RD )3 ) are equal to is 7.3 × 10−3 and 304, respectively. Because of the high difference between both the contraction ratios and De numbers, this situation is definitely favorable for exposing the effect of jet flow in

the near field reservoir entrance region, as compared to the slit channel case. Fig. 6 shows very good agreement with experimental data for the die flow within the range of De numbers from 6 to over 300. Fig. 7A and B demonstrates the centerline velocity and pressure profiles along the circular channel at De ≈ 250 using the jet flow models with and without drag, respectively. Because there is a big difference between the lengths of the two modeled reservoir jets, we changed the scale for negative (reservoir) values of longitudinal variable x1 in Fig. 7A. Our calculations qualitatively agree with the experimental results obtained by Paskhin [1] who demonstrated that the centerline velocity goes first through a maximum and finally stabilizes. Direct comparison with data [1] was impossible because Paskhin used a reservoir with a conical entrance and also did not report important material characteristics. The comparison between Fig. 7A and B show that the simple jet model predicts a longer length of reservoir jet and shorter length of the developing die flow region than the more realistic model of the reservoir jet with drag. For example Fig. 7A shows that the longer dimensionless length ≈16of the developing flow (evaluated by the axial velocity profile) is needed for the jet with drag flow model as compared to that ≈5 for the jet flow only. The reason for this difference is the longer relaxation of higher values of initial

J.-H. Jeong, A.I. Leonov / J. Non-Newtonian Fluid Mech. 118 (2004) 157–173 1e+5

1e+5

1e+4

1e+4 G ( Pa )

G' & G" ( Pa )

168

1e+3 1e+2 w - G' w - G" calculations

1e+1

1e+2 experiments calculations

1e+1

1e+0

1e+0 0.01

0.1

1

10

100

0.01 0.1

ω ( rad /s )

(B)

elongational viscosity (η E: Pa*s)

(A)

0.0034 viscosity: ln(η) ( Pa*s )

1e+3

0.0032

0.0030

0.0028

8

10 1/Temp. ( K )

(C)

10

100

1e+6

1e+5 -3

-1

el.rate:9.263e s el.rate:1.077e-2 s-1 -2 -1 el.rate:2.035e s

1e+4 0

12 -1

1

time ( sec )

5

10

time (sec)

(D)

Fig. 5. Basic rheological experimental data and their modeling for PIB (LM-MH) for: (A) dynamic test at 50 ◦ C; (B) stress relaxation test at 50 ◦ C; (C) temperature dependence of Newtonian viscosity; (D) uniaxial extension test at 18 ◦ C. Various symbols in the figures indicate experimental data, lines denote calculations. 10000

app.viscosity(Pa*s)

MPT experiment calc ulation w/o drag calc ulation w/ drag

1000 10

100

1000

app. shear rate (1/s)

(A)

pressure(Pa)

1.5 e+7

1.0 e+7

5.0 e+6 MPT experiment cal culat ion w/o drag calculation w/ drag 0.0 e+0 0. 0e+0

(B)

2.0e-8

4.0e-8

6.0e-8

8.0e-8

1.0e-7

1.2e-7

3

flow rate (m /s)

Fig. 6. Isothermal circular capillary flow for PIB (LM-MH) at 50 ◦ C in Monsanto processability tester (geometrical description is given in Table 3): (䊊) experimental data; solid lines, jet calculations without drag; dashed lines, with drag. (A) Plots of apparent viscosity vs. apparent shear rate; (B) plots of transducer T* pressure (Fig. 1B) vs. flow rate.

axial velocity at x2=0 (m/s)

J.-H. Jeong, A.I. Leonov / J. Non-Newtonian Fluid Mech. 118 (2004) 157–173

169

0. 05

calculation w/o drag calculation w/ drag

0. 00

- 40

- 20

0 x 1 /R

(A)

10

20

-4

-2

0

10

x 1 /R

D

20

D

1e+7

calculation w/o drag calculation w/ drag

Press.(Pa)

8e+6 6e+6 4e+6 2e+6 0e+0 0

(B)

5

10 x 1 /R

15

20

25

0

5

10 x 1 /R

D

15

20

25

D

50 ◦ C

Fig. 7. The flow in the circular rheometer for PIB (LM-MH) at under flow rate, Q = 5.0 × 10−8 m3 /s: solid lines, jet calculations without drag; dashed lines, jet calculations with drag. (A) Centerline axial velocity profiles in the reservoir and die regions; (B) pressure profile along the die region.

elastic strains at the die entrance, which is characteristic for the model of jet flow with drag. One can also see from the calculations of both models that the length of developing flow region evaluated from pressure profile is considerably less than that evaluated using the axial velocity profile. It indicates that the stabilization of pressure gradient does not completely describe the developing flow region. 5.2. Experiments and calculations for LDPE The contraction flow in Ref. [28] was investigated using LDPE (Lupolen 1810 H from BASF Ludwigshafen). It should be noted that LDPE needs special modeling consideration [24,29] for the function b in Eqs. (1) and (2). This type of modeling is shown in Table 1. The experiments [28] were performed using 10:1 (reservoir/die) slit contraction controlled by constant plunger velocities at a temperature of 150 ◦ C. At this temperature, the density of the material is ρ = 778 kg/m3 and the zero-shear-rate viscosity is η ≈ 65 kPa s. The average axial velocities in the die were ranged from 4.6 to 49.3 mm/s. No instability was observed in these typical creeping flow experiments. Laser Doppler velocimetry was used to measure 3D velocity field. The detailed material characterization and 10:1 slit dimensions are shown in Tables 1 and 4. In the present paper, instead of directly employing the 14 relaxation modes found in [28], we used

eight modes for characterizing LDPE as found by Laun [30]. It was impossible to use the Pade–Laplace method [27] to extract Maxwell modes from the experimental data in linear region, because the stress relaxation for long time durations was unavailable. Nevertheless, using eight modes from [30] and carefully chosen non-dimensional parameters shown in Table 1, we were able to describe reasonably well all basic, shearing and extensional rheological experiments (Fig. 8). With our selection of Maxwell modes, the  average re laxation time was calculated as θ¯ = k Gk θk2 / k Gk θk = 58.7 s. Then the De numbers in the slit die were defined as: Table 4 Material parameters and experimental conditions for the contraction flows with geometry shown in Fig. 1A Polymer melts

LDPE Melt I (Lupolen 1810 H from BASF Ludwigshafen) [28]

Parameters at 423 K [30] θ k (sec)

Gk , (KPa)

1000 100 10 1

0.001 0.18 1.89 9.8

0.1 0.01 0.001 0.0001

26.7 58.6 94.8 129

Geometry (cm)

L = 2.0 Reservoir: 2LR = 2.0, la = lb = 10.0, lc = 0 Die: 2LD = 0.2, Lz = 5.0

1e+5

1e+7

1e+4

1e+5 1e+4

1e+3

N 1 (Pa )

1e+6

1e+3 1e+2

1e+2 0.01

1

10

100 1000

Q = 1.02e-6 (m3/s)

∆ Q = 2.25e-6 (m3/s)

-50

0

50

1e+3

O G’ G” 0.1

1

10

100

1000

ω (1/s)

axial vel. (m/s) at x2=0

1e+4

0.01

x1/LD

(A)

1e+5

(B)

0.10 0.08 0.06 0.04 0.02 0.00 -20

0

1e+6 1e+5

.

O ε. = 0 .1 s -1-1 ε. = 0 .3 s ∆ ε = 1 .0 s -1

1e+4 0.1

1

10

100

time (s)

De = θ¯ γ˙ App = θ¯ × 6Q/(4LL2D ). Here γ˙ App is the apparent shear rate in slit die, Q the flow rate, and L and LD are the width and the half thickness of the die, respectively. Because of unclear definition of experimental flow rates (or De numbers) in Ref. [28], we recalculated the flow rates from the gap-wise axial velocity profile as: Q = 4.00 × 10−7 , 1.02 × 10−6 and 2.25 × 10−6 m3 /s for the low De number flows. Fig. 9 illustrates comparison between our calculations and experimental data for centerline axial velocity for flows with three different Deborah numbers. Our calculations compared to the experimental data presented in Fig. 9A used the pure jet approach and in Fig. 9B the jet with drag. Comparing Fig. 9A and B indicates that the most significant difference in model predictions is the jet length in the reservoir, the more detailed approach (jet with drag) showing a good agreement with experiments (Fig. 9B). Definitely, the above mentioned experimental conditions are favorable for implementing our model, but the local Deborah numbers of reservoir shear flow are not high enough for initial extensional jet to be dominant. Thus, there is considerable difference between the converging lengths for the models of jet with or without drag. But interestingly enough, both the models for the die developing flow agree well with experimental data.

(C)

40

60

40

60

0.10 0.08 0.06 0.04 0.02 0.00 -20

Fig. 8. Basic rheological experimental data and the CE predictions for LDPE: (A) steady simple shear viscosity and the first normal stress difference N1 at 150 ◦ C, (B) dynamic test at 150 ◦ C; (C) uniaxial extension test at 150 ◦ C: various symbols stand for experimental data and lines for calculations.

20 x1/LD

(B) axial vel. (m/s) at x2=0

elongational viscosity (η E:Pa)

O Q = 4.00e-7 (m3/s)

-100

1e+2

(C)

0.10 0.08 0.06 0.04 0.02 0.00

-1

shear rate (s )

(A)

G' & G" (Pa)

0.1

axial vel. (m/s) at x2=0

J.-H. Jeong, A.I. Leonov / J. Non-Newtonian Fluid Mech. 118 (2004) 157–173

visco sity (Pa *s)

170

0

20 x1/LD

Fig. 9. Experimental data for LDPE and calculations for the centerline axial velocity under various given flow rates at 150 ◦ C; the origin of x1 locates at the die entrance. (A) Computations using the jet approach; (B) computations using the jet approach with drag; (C) direct computation and experimental data [28]. Symbols stand for experimental data and lines for calculations.

The comparison between direct numerical simulations and their experiments [28] is presented in Fig. 9C. Comparing Fig. 9B and C, one can see that our calculations are slightly better for flows in the die entrance and in the region of developing die flow. Although both simulations evaluate the non-dimensional axial length of developing flow as 10–20 in x1 /LD units, the direct simulations underestimate the die developing behavior, especially for higher values of flow rate. It seems that the deviation of the direct simulation from experimental data increases with the increase in flow De number. Fig. 10 demonstrates the velocity profiles (dots) at x1 /LD = 40 for each flow De. The results of our computations presented by lines are obtained using the model of jet with drag. Unfortunately, the measurements of velocity profiles were only performed at the land (developed flow) region. Therefore, the detailed developing behavior of die flow, the most interesting in the contraction flow, could not be seen from Fig. 10.

J.-H. Jeong, A.I. Leonov / J. Non-Newtonian Fluid Mech. 118 (2004) 157–173 O Q = 4.00e-7 (m3/s) • Q = 1.02e-6 (m33/s) ∆ Q = 2.25e-6 (m /s)

axial vel. (m/s)

0.10

0.05

0.00 -1

0

1

x2/LD

Fig. 10. Gap-wise axial velocity profile with various flow rates. Symbols stand for experimental data and lines for calculations.

171

tational analyses of non-isothermal contraction flows are relatively easy to carry out using the present model. An example of such an analysis will also be presented. Additionally (vi) other types of viscoelastic flows with geometric and physical complexities, such as non-symmetrical entrance flows and polymer wall slip, can also be treated relatively easily when using the present model. It should be finally mentioned that numerical analyses of high Deborah number polymer flows are constrained by the onset of melt flow instabilities. Two basic instabilities of either entrance or wall slip types, could also be approximately analyzed within the present model.

6. Conclusions Although there are a lot of published data described contraction flows of polymeric fluids, they do not provide necessary information for our modeling. Therefore, we were unable to compare our model calculations with many direct calculations and/or experimental data. Nevertheless, the comparison of our calculations with available literature results demonstrated that our model captures the most essential features of contraction flow with considerably less computational effort than that used in direct computations. Definitely, our models have limitations in cases when there is no dominant extensional flow in the near field entrance region. These cases are the low De number flows and/or the flows with small contraction ratio. Those cases are, however, rarely applied to real polymer processing and are outside the scope of our modeling. As discussed above, we ignored in our modeling realistic smooth transitions from one characteristic flow zone to another, substituting them with “sudden” changes. Therefore, the model cannot describe details of flow in the transitional zones. However, the evaluation of the change in flow characteristics using matching conditions, was successful for the flow regions in succession. It seems that our approach works better for higher De flows. This happens because of the shortening of the transitional regions between the basic flow zones. With due respect to direct numerical computations, there are remarkable advantages in our model predictions, such as (i) better agreement with experimental data at higher De numbers (high-speed processing) and higher contraction ratios, which are well known to be troublesome conditions for existing direct computational methods. Therefore, the present model could be a good candidate for modeling steady polymer contraction flows at high De numbers. Additionally, our model is (ii) stable when employing a complicated but very realistic rheological model, and (iii) effective in treating flow in a geometrically simplified manner. (iv) The model operates with well-defined rheological parameters and does not use any adjustable parameters. Since the model approximately (v) traces entire deformation history of the flow, it is possible to apply it for calculating swelling in short dies. This analysis will be demonstrated elsewhere. Also compu-

Appendix A. Effective matching conditions for different types of steady channel viscoelastic flows We start from the exact local balance of the mechanical energy: ∂ ∂ ∂ (K + W) + [vi (K + W)] + D = (vj σij ). ∂t ∂xi ∂xi

(A.1)

Here K = ρ|v|2 and W are the local kinetic and elastic energies, D is the dissipation, v and σ are the velocity and total stress. Eq. (A.1) immediately follows from the momentum balance and general viscoelastic constitutive equations. In cases of flows of polymer melts, the contribution of kinetic energy K in energy balance can always be neglected. Then Eq. (A.1) can be split in the case of constitutive Eqs. (1) and (2) into independent energy balances for each kth non-linear relaxation modes: ∂ ∂ ∂ [vi (Wk )] + Dk = (vj σij,k ). (A.2) (Wk ) + ∂t ∂xi ∂xi Here Wk , Dk , and σ ij ,k are, respectively, the elastic potential, dissipation, and the full stress tensor in the kth relaxation mode. Consider now a steady flow of an incompressible viscoelastic liquid in a composite cylindrical tube whose arbitrary cross-section(s) is represented by a domain Ω. Then integrating Eq. (A.2) over the tube cross-section with the use of non-slip condition at the wall yields    d d v(Wk )dω = vσ11,k dω − Dk dω. (A.3) dx1 Ω dx1 Ω Ω Here v≡v1 is the longitudinal component of velocity. It should be mentioned that in any steady developed flow, the right-hand side of Eq. (A.3) vanishes. We now assume that in the region x1 < x0 the flow is of a certain type, while in the region, x1 > x0 + δ, it is of another type, with a transient flow in the region x0 < x1 < x0 + δ. Here x1 is the longitudinal axis in flow region. If outside the transition region the flows could be approximately considered as steady, one can neglect the right-hand side in Eq. (A.3). Then integrating (A.3) over the transition area and neglecting the length of the transient zone (δ → 0) will

172

J.-H. Jeong, A.I. Leonov / J. Non-Newtonian Fluid Mech. 118 (2004) 157–173

finally result in the jump-like continuity condition for the “free energy flux” in any kth relaxation mode:   vWk ds|x0 −0 = vWk ds|x0 +0 . (A.4) Ω

Ω

The condition of (A.4) type has been used by Graessley et al. [26] for the evaluation of extrudate swell.

Appendix B. Evaluation of stress singularity in modeling of the secondary flow near the upstream corner To simplify the evaluation of this singularity we consider in the expressions (20) and (20*) the drag from the circulatory Newtonian flow, and express the shear stress at the jet surface as LEx 4η¯v1 /LR σ12,surf = − , ξ= , (B.1) 1−ξ LR η¯v1 ∗ = σ12,surf RR ξ ∗



 (1 − ξ ∗2 )(1 − 3ξ ∗2 ) − 4ξ ∗4 ln(ξ ∗ ) , (1 − ξ ∗2 )[(1 − ξ ∗2 ) + (1 + ξ ∗2 )ln(ξ ∗ )] REx . (B.1*) ξ∗ = RR

Here σ 12,surf is the shear stress at the jet surface, η the Newtonian viscosity, v¯ 1 the averaged axial velocity at x1 = −l, LEx (REx ) is half thickness (radius) of jet at x1 = −l + dx1 , LR (RR ) is the half thickness (radius) of slit (circular) reservoir. Inserting (B.1, B.1*) into (19, 19*) and rearranging the formulae with the aid of (9) and (11), yields    3 dξ ∗ 2 − G k λk + 2 λk dx1 k    Gk 4η¯v1 1 4 λk − 4 , ξ → 1, = 2 + 2θk v¯ 1 LR (1 − ξ ∗ ) λk k (B.2) 

−

k

  2 dξ ∗ Gk λ2k + λk dx1

   Gk η¯v1 1 = 2 2λk + 2 + 12θk v¯ 1 RR (1 − ξ ∗ ) λk k   1 × λ2k − (λk + 1), ξ ∗ → 1. λk

(B.2*)

Integrating Eqs. (20) and (20*) in the limits ξ → 1, ξ ∗ → 1 yields  ξ ≈ 1 − α)x10 ;  −1 Gk (λ2k + 3λ−2 α = 8η¯v1 (L2R k )|λk =λl ) k

(α)x10  1),

k

(B.3)

ξ∗ ≈ 1 −



α∗ )x10 ;  −1 Gk (λ2k + 2λ−1 α∗ = 2η¯v1 (R2R k )|λk =λl ) k

(α∗ )x10  1).

k

(B.3*)

Here λlk are the initial stretching values at the point x1 = −l. Once the initial gap distance of secondary shear flow is determined from (B.3, B.3*) for a small numerical step increase )x10 , the iterative numerical method described in the main text, can be implemented. The numerical procedure has been performed using a root finding subroutine for the unknown gap distance at x1 = −l + )x10 . And the shear stress at the jet surface x1 = −l + )x10 has been computed from Eqs. (B.1) and (B.1*) with the calculated gap distance value, ξ (or ξ ∗ ).

References [1] E.D. Paskhin, Motion of polymer liquids under unstable conditions and in channel terminals, Rheol. Acta 17 (1978) 663. [2] S.J. Muller, E.S. Shaqfeh, R.G. Larson, Experimental studies of the onset of oscillatory instability in viscoelastic Taylor–Couette flow, J. Non-Newtonian Fluid Mech. 46 (1993) 315. [3] B.S. Baumert, S.J. Muller, Flow visualization of the elastic Taylor–Couette instability in Boger fluids, Rheol. Acta 34 (1995) 147. [4] S.A. White, D.G. Baird, Flow visualization and birefringence studies on planar entry flow behavior of polymer melts, J. Non-Newtonian Fluid Mech. 29 (1988) 245. [5] R.J. Binnington, G.J. Troup, D.V. Boger, A low cost laser-speckle photographic technique for velocity measurement in slow flows, J. Non-Newtonian Fluid Mech. 12 (1983) 255. [6] T.D. Dudderar, P.G. Simpkins, Laser speckle photography in a fluid medium, Nature 270 (1977) 45. [7] A.I. Isayev, R.K. Upadhyay, Two-dimensional viscoelastic flows: experimentation and modeling, J. Non-Newtonian Fluid Mech. 19 (1985) 135. [8] K.-J. Bathe, Finite Element Procedures in Engineering Analysis, Prentice Hall, New York, 1982. [9] S.A. White, A.D. Gotsis, D.G. Baird, Review of the entry flow problem: experimental and numerical, J. Non-Newtonian Fluid Mech. 24 (1987) 121. [10] R.A. Brown, R.C. Armstrong, A.N. Beris, P.W. Yeh, Elastic flows, Comput. Math. Appl. Mech. 58 (1986) 201. [11] R. Keunings, On the high Weissenberg number problem, J. Non-Newtonian Fluid Mech. 20 (1986) 209. [12] S. Dupont, M.J. Crochet, Finite element simulation of viscoelastic fluids of the integral type, J. Non-Newtonian Fluid Mech. 17 (1985) 157. [13] R.K. Upadhyay, A.I. Isayevy, Simulation of two-dimensional planar flow of viscoelastic fluid, Rheol. Acta 25 (1986) 80. [14] E. Mitsoulis, Numerical simulation of entry flow of the IUPAC-LDPE melt, J. Non-Newtonian Fluid Mech. 97 (2001) 13. [15] E. Mitsoulis, M. Schwetz, H. Munstedt, Entry flow of LDPE melts in a planar contraction, J. Non-Newtonian Fluid Mech. 111 (2003) 41. [16] F.N. Cogswell, Converging flow of polymer melts in extrusion dies, Polym. Eng. Sci. 12 (1972) 64. [17] D.M. Binding, An approximate analysis for contraction and converging flows, J. Non-Newtonian Fluid Mech. 27 (1988) 173.

J.-H. Jeong, A.I. Leonov / J. Non-Newtonian Fluid Mech. 118 (2004) 157–173 [18] C. Carrot, J. Guillet, R. Fulchiron, Converging flow analysis, entrance pressure drops, and vortex sizes: measurements and calculated values, Polym. Eng. Sci. 41 (2001) 2095. [19] J. Guillet, P. Revenu, Y. Bereaux, J.-R. Clermont, Experimental and numerical study of entry flow of low-density polyethylene melts, Rheol. Acta 35 (1996) 494. [20] V.G. Levich, Physicochemical Hydrodynamics, Prentice Hall, New York, 1962. [21] N.A. Sliozkin, Dynamics of Incompressible Viscous Fluids, Gos. Tech. Theor. Izdat., Moscow, 1955, p. 357. [22] S.M. Richardson, Fluid Mechanics, Hemisphere Publishing Co., 1987, p. 182. [23] M. Simhambhatla, A.I. Leonov, On the rheological modeling of viscoelastic polymer liquids with stable constitutive equations, Rheol. Acta 34 (1995) 259. [24] A.I. Leonov, in: D.A. Siginer, D. De Kee, R.P. Chabra (Eds.), Advances in the Flow and Rheology of Non-Newtonian Fluids, Elsevier, New York, 1999, pp. 519–576.

173

[25] Y. Kwon, A.I. Leonov, Stability constraints in the formulation of viscoelastic constitutive equations, J. Non-Newtonian Fluid Mech. 58 (1995) 25. [26] W.W. Graessley, S.D. Glassock, R.L. Crawley, Die swell in molten polymers, Trans. Soc. Rheol. 14 (1970) 519. [27] M. Simhambhatla, A.I. Leonov, The extended Pade–Laplace method for efficient discretization of linear viscoelastic spectra, Rheol. Acta 32 (1993) 589. [28] K. Feigl, H.C. Ottinger, A numerical study of the flow of a low-density-polyethylene melt in a planar contraction and comparison to experiments, J. Rheol. 40 (1996) 21. [29] M. Zatloukal, J. Vlcek, C. Tzoganakis, P. Sa’ha, in: Proceedings of the 6th European Conference on Rheology: Molecular Modelling, Numerical Simulation, and Micro-structural Modeling, 2002, pp. 247–248. [30] H.M. Laun, Description of the non-linear shear behavior of a low density polyethylene melt by means of an experimentally determined strain dependent memory function, Rheol. Acta 17 (1978) 1.