Numerical simulation of three-dimensional viscoelastic flow using the open boundary condition method in coextrusion process

J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

Numerical simulation of three-dimensional viscoelastic flow using the open boundary condition method in coextrusion process Ki Byung Sunwoo a , Seung Joon Park a , Seong Jae Lee b , Kyung Hyun Ahn a , Seung Jong Lee a,∗ a

b

School of Chemical Engineering, Seoul National University, Seoul 151-744, South Korea Department of Polymer Engineering, The University of Suwon, Suwon 445-743, South Korea Received 6 November 2000; received in revised form 19 March 2001

Abstract Three-dimensional numerical simulation of coextrusion process of two immiscible polymers through a rectangular channel has been performed using the finite element method. The upper convected Maxwell (UCM) model and the Phan-Thien and Tanner (PTT) model were considered as viscoelastic constitutive equations. The elastic viscous stress splitting (EVSS) method was adopted to treat the viscoelastic stresses, and the streamline upwinding (SU) method was applied to avoid the failure of convergence at high elasticity. The problem arising from the ambiguous outlet boundary condition that has previously been used in the three-dimensional simulation of a viscoelastic coextrusion process could be avoided by introducing the open boundary condition (OBC) method. The abrupt change or deviation of contact line position near the outlet that was observed when the fully developed outlet boundary condition was applied could be clearly removed by using the OBC method. The effects of viscoelastic properties, such as the shear viscosity ratio, the elasticity, the second normal stress difference, and the extensional viscosity on the interface distortion, the interface curvature, and the degree of encapsulation along the downstream direction have been investigated. The shear viscosity ratio between the polymer melts was the controlling factor of the interface position and the encapsulation phenomena. The interface distortion seems to increase as the elasticity ratio increases under constant shear viscosity, even though it is not so large. The degree of encapsulation seems to increase with increasing the ratio of the second normal stress differences. The extensional viscosity had minor effect on the encapsulation phenomena. The second normal stress difference was found to have a great influence on the increasing of the degree of encapsulation along the downstream direction as compared to the effect of the first normal stress difference. © 2001 Elsevier Science B.V. All rights reserved. Keywords: Coextrusion process; Three-dimensional numerical simulation; Finite element method; Viscoelastic fluid; Open boundary condition method; Viscosity ratio; Second normal stress difference

∗

Corresponding author. Fax: +82-2-888-7295.

0377-0257/01/$ – see front matter © 2001 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 7 - 0 2 5 7 ( 0 1 ) 0 0 1 1 5 - X

126

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

1. Introduction In recent years, multicomponent coextrusion has gained wide recognition as an approach to achieving unique product performance by combining the properties of different materials with lower expenses. Cylindrical dies are mainly used in the production of conjugate and sheath–core fibers and annular dies in the coextruded wire or cable or blown film industry. Multilayer polymer sheets or films are normally produced by the coextrusion of two or more polymers from flat dies. This is accomplished using either of two possible die configurations. In the first, the die is made up of manifolds for each polymer layer, and the layers are joined together just before they exit the die. In the second, one die is used for all layers and a coextrusion feedblock is used before the entrance of the die. The purpose of the feedblock is to combine and distribute various layers such that they have a uniform thickness across the exit of the feedblock. Therefore, the design of both feedblock and die requires the understanding of three-dimensional flow behavior in the coextrusion process. Experimental investigation of a stratified flow in a side-by-side coextrusion has so far identified the viscosity difference between the polymer melts to be the controlling factor of the encapsulation phenomena, i.e. the less viscous melt wraps around the more viscous melt [1–4]. Karagiannis et al. [5] showed that significant interface deformation occurs during the flow in the square die at relatively low viscosity ratio and that the interface curvature increases gradually. The effect occurs near the die walls and may result in undesirable product properties, which increases manufacturing cost due to the need to trim the film or sheet edges. Southern and Ballman [3] investigated the relative importance of viscosity and elasticity effect on the interface shape. They showed that a significant difference in elasticity of the two components produces a small change of interface shape relatively to that observed with a similar viscosity ratio shift, and found that melt elasticity difference causes uneven and ragged interface shape. For the theoretical investigation, White et al. [6] applied the theory of nonlinear viscoelastic fluids to consider the stress and velocity profiles and the interface shape in stratified flows, and showed that the fluid with the greater absolute second normal stress difference tends to protrude into the other fluid. In contrast, Khan and Han [7] showed that the more elastic component tends to wrap around the other one using the Coleman and Noll’s second-order fluid model. Song [8] showed by numerical stability analysis that the stable sheath–core configuration can be obtained when the core consists of higher viscosity and elasticity, while the sheath consists of lower viscosity and elasticity. The effect of polymer viscoelasticity on the uniformity of the layer thickness in the multilayer coextruded structures was investigated with the same materials in each layer to minimize the viscosity effect. The development of secondary motion due to the non-zero second normal stress difference was investigated by observing the flow of multilayered polymer systems [9,10]. Analytical studies were carried out on the basis of a Criminale–Ericksen–Filbey fluid model by Dodson et al. [11] and Townsend et al. [12]. A comparison between experiments and finite element predictions was also made by Debbaut et al. [13] and Debbaut and Dooley [14]. Debbaut et al. studied the straight and tapered channel flows with a square cross-section. Debbaut and Dooley focused on the straight channels with a rectangular or teardrop shaped cross-section, and showed that a reduction of the second normal stress difference leads to less pronounced secondary motions. The modeling and simulation of coextrusion process requires simultaneous consideration of mass and momentum principles as well as constitutive relations. The geometry of the flow domain enters through the boundary conditions. The complex die geometry and the nonlinear polymer melt property result in a mathematical problem that must be solved numerically. Mavridis et al. [15] implemented a

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

127

double-node finite element technique for the stratified multiphase flow simulation to treat the discontinuity of pressure at the interface. Agassant and Fortin [16], who simulated the coextrusion process of Carreau fluids in two-dimensional formulation, pointed out that the interface location is strongly dependent upon the flow rate and the flow rate ratio, but they could not describe the encapsulation phenomena. Karagiannis et al. [5] and Gifford [17] studied isothermal generalized Newtonian coextrusion process using the three-dimensional analysis. Takase et al. [18] performed a three-dimensional numerical simulation with a viscoelastic model to explain the encapsulation phenomena. They showed that the encapsulation phenomena are affected not only by viscous properties but also by elastic or nonlinear properties. However, the fully developed boundary condition they imposed at the outlet of the analysis region was not correct. One important aspect which must be considered when solving the boundary value problem numericals is the treatment of boundary conditions. In some cases, the boundary can be an actual boundary of the spatial domains under consideration, such as the no-slip condition at walls, while in other cases the boundary may be an artificial one which is chosen for the calculation purpose, such as the outflow velocity profile at the outlet of the die. The choice of good boundary conditions that have physical significance and the way to combine this condition with the numerical scheme employed in the interior is an important subject of research [19,20]. In the numerical simulation of a viscoelastic flow inside a die, velocity profiles are required along the outflow boundary. In the isothermal problem, numerical simulation is normally carried out with the long enough downstream domains, especially when the flow is fast. If the long domain is used, the velocity is fully developed at the outflow boundary and the outflow boundary condition can be matched with the solution calculated in the domain of interest [21]. However, in the nonisothermal case for example, because the velocity profiles are strongly dependent upon temperature and the length required for fully developed profile is very long, it is almost impossible to impose outflow velocity profiles matched with the solutions calculated in the upstream. In order to avoid this difficulty, Wapperom et al. [22] used the solution calculated at the outflow in one-dimensional problem as the outflow velocity boundary condition. They calculated the outflow velocity profile with the previously calculated temperature profile. However, the method they used can only be applied to the decoupled method, which calculates the extra stresses, velocity, and temperature separately. Therefore, if a coupled method is used, this method cannot be applied to the numerical simulation of nonisothermal viscoelastic flow problems. In the case of coextrusion of viscoelastic fluids inside a die, velocity profile should also be imposed at the outflow boundary. However, the velocity profile is not known a priori as in the nonisothermal flow problem. This makes the imposition of proper boundary condition difficult in the coextrusion process of viscoelastic fluids. The imposition of outflow boundary condition problem can be avoided by using the open boundary condition (OBC) method, which is often called as ‘no boundary condition’ or ‘free boundary condition’ method; i.e. not imposing the condition along the outflow boundary. Papanastasiou et al. [23] applied the OBC method to the backward facing step problem at high Reynolds number. OBC method allows the solution calculated in the upstream to pass through the outflow boundary without undergoing significant distortion or without influencing the interior solution. In the OBC method, a line (or surface) integral along the outflow boundary is calculated by unknown outflow nodal values, which are obtained as a part of solutions. Park and Lee [24,25] suggested using the OBC method in the two-dimensional numerical simulation of isothermal and/or nonisothermal viscoelastic flow. The authors [26] performed three-dimensional numerical simulation of nonisothermal coextrusion process with the generalized Newtonian fluid model using the OBC method.

128

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

This work focuses on the effects of the viscoelastic properties on the encapsulation phenomena in a coextrusion process of two immiscible fluids through a rectangular channel. The OBC method will be implemented to remove the difficulty of imposing the outlet boundary condition. The effect of viscosity ratio and elasticity on the interface shape will be investigated with the upper convected Maxwell model. In addition, by using the Phan-Thien and Tanner model [27], the effect of the shear viscosity ratio, the elasticity, the extensional viscosity, and the second normal stress difference will also be considered. 2. Problem description and viscoelastic constitutive models The bicomponent stratified flow of immiscible fluids involves the merging of two fluid streams, the flow inside the die and the flow of the bicomponent system out of the die. In this study, the problem consists of the merging flow of the two fluid streams and the interface shape development in the resulting bicomponent stratified flow inside the die. The three-dimensional continuity and momentum equations for the steady state flow of two viscoelastic fluids (I and II) are ∇ · vk = 0,

k = I, II,

ρk vk · ∇ vk = −∇pk + ∇ · τk ,

(1) k = I, II,

(2)

where v is the velocity vector, τ the extra stress tensor, p the pressure, and ρ the density. The boundary conditions at the interface are expressed as follows. Kinematic conditions: n · vI = n · vII = 0,

(3)

t1 · vI = t1 · vII ,

(4)

t2 · vI = t2 · vII ,

(5)

Dynamic conditions: t1 · σI = t1 · σII ,

(6)

t2 · σI = t2 · σII ,

(7)

n2 · σI − n2 · σII = 0,

(8)

where σ is the total stress tensor, n the unit normal vector at the interface, and t the unit tangential vector at the interface. In dynamic conditions, surface tension force is neglected, which is generally acceptable due to the high viscosity of polymer melts. 2.1. Upper convected Maxwell model The upper convected Maxwell (UCM) model as a viscoelastic constitutive equation is the one most frequently used in the viscoelastic simulations, which is written as ∇

τ + λ τ = η(∇ v + ∇ vT ),

(9)

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

129

where λ is the relaxation time and η the shear viscosity. The inverse triangular symbol ∇ applied over the extra stress tensor stands for the upper convected time derivative, which is defined as ∇

τ=

∂τ + v · ∇ τ − ∇ vT · τ − τ · ∇ v. ∂t

(10)

To treat the viscoelastic stresses, the elastic viscous stress splitting (EVSS) method [28] is used, where the extra stress tensor is divided into the elastic and viscous parts, i.e. τ = S + 2ηD,

(11)

where S is the modified extra stress tensor, and D the deformation rate tensor defined by D = 21 (∇ v + ∇ vT ).

(12)

The Galerkin finite element method is used to discretize the above set of differential equations. The velocity, pressure, modified stress, and deformation rate unknowns are approximated in terms of Galerkin basis functions as follows:

v= ψ i vi , p= φi pi , (13) i

i


S= φi Si ,


D= φi Di ,

i

i

(14)

where φ i and ψ i are the trilinear and triquadratic basis functions, respectively. The governing equations, weighted with the basis functions and integrated by parts, result in the following weak forms: ∇

∇

S + λ(S + 2ηD); φ = 0,

(15)

ρ v · ∇ v; ψ + −pI + S + η(∇ v + ∇ vT ); ∇ψ = −pI + S + η(∇ v + ∇ vT ) · n; ψ ,

(16)

∇ · v; φ = 0,   D − 21 (∇ v + ∇ vT ); φ = 0,

(17) (18)

where I is the unit tensor. In Eqs. (15)–(18), symbols ; and ; stand for domain boundary surface integrals, respectively. 2.2. Phan-Thien and Tanner model The Phan-Thien and Tanner (PTT) model [27] as a viscoelastic constitutive equation is the one used frequently to represent the second normal stress difference that cannot be obtained by the UCM model. The PTT model is written as      ελ ξ ∇ ξ 1 + trτ τ + λ 1 − τ + τ = 2η0 D, (19) η0 2 2 where λ is the relaxation time, η0 the zero shear viscosity, ξ the parameter indicating the second normal stress difference property, and ε the parameter indicating the extensional viscosity property. The triangular

130

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

symbol  applied over the extra stress tensor stands for the lower convected time derivative, which is defined as ∂τ  + v · ∇ τ + ∇ v · τ + τ · ∇ vT . τ= (20) ∂t The final forms of governing equations, weighted with the basis functions and integrated by parts, result in the following weak forms:       ελ ξ ∇ ξ S+ trS S + (2ελtrS)D + λ 1 − S+ S η 2 2 0     ξ ∇ +λ 2η0 1 − (21) D + η0 ξ D ; φ = 0, 2 ρ v · ∇ v; ψ + −pI + S + η(∇ v + ∇ vT ); ∇ψ = −pI + S + η(∇ v + ∇ vT ) · n; ψ ,

(22)

∇ · v; φ = 0,   D − 21 (∇ v + ∇ vT ); φ = 0.

(23) (24)

In the meanwhile, Takase et al. [18] used the modified PTT model to include Newtonian shear viscosity component to the total shear viscosity, which is defined as τ = 2η0 s D + E,





1+

ελ trE E + λ η0



 1−



 ξ ∇ ξ E + E = 2η0 (1 − s)D, 2 2

(25) (26)

where s is a material parameter. They applied decoupled method to calculate the velocity and the stress fields separately, while we used the coupled method to solve the velocity and the stress fields simultaneously. 3. Numerical methods 3.1. Upwinding method To avoid the failure of convergence at high elasticity, the streamline upwinding (SU) method is applied, where the upwinding effect is imposed only on the convective terms in the constitutive equation. Thus, SU method is an inconsistent formulation unlike the streamline upwind/Petrov–Galerkin (SUPG) method, where upwinding effect is given to the constitutive equation consistently. Debae et al. [29] compared the results of SU, SUPG, and standard Galerkin method using the EVSS formulation in four benchmark problems using the UCM fluid, i.e. the flow around a sphere, through a wavy tube, through an abrupt contraction, and in circular extrusion. Referring to their results for contraction flow problem, despite the artificial diffusion of SU method, the results of SU method were in good agreement with those of standard Galerkin and SUPG methods. In SU method, the modified weighting function φ m is given by φm = φ + k¯ v · ∇φ.

(27)

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

131

The definition of the factor k¯ in Eq. (27), which depends on the velocity and element size, is given by Marchal and Crochet [30] for the two-dimensional case. Here, we simply extended their method to the three-dimensional case. The flow domain is discretized into 27-node hexahedron elements, and the Galerkin finite element procedure is applied [31]. The set of nonlinear algebraic equations is then finally solved by means of Newton’s method and the frontal elimination technique. 3.2. Discontinuity of pressure and extra stresses In a viscoelastic coextrusion process, the pressure and extra stresses are discontinuous along the interface. This discontinuity is handled by assigning two pressure and two extra stress variables per node along the interface, which is called as a double-node technique [5,15]. In the double-node approach, two pressure and two extra stress variables are assigned at each relevant vertex of the three-dimensional elements along the interface. 3.3. Interface update scheme and remeshing procedure The problem is complicated because of the existence of fluid/fluid interface whose position is not known a priori and must be included in the calculation by Picard-type iteration. Pathline method is used to update the interface position [5,15,17]. The interface is described by the path of particles travelling downstream from the end of the separating plate, because particles at the interface remain in a tangential direction to the interface. Their paths are determined by the following relationship: u v w = = , (28) dx dy dz and the interface is calculated separately from the calculation of velocity field based on the latest information of converged velocity as follows:

v dz, (29) y = w where x, y, and z are the coordinate directions, and u, v, and w the velocity components in x, y, and z directions, respectively. The pathline method is a simple method but gives slow convergence, and thus the procedure generally requires a large number of iterations to achieve a special convergence. After updating the interface position, it is necessary to relocate the inner coordinates of the flow domain of interest. We remeshed the flow domain only in thickness direction using the predetermined ratio. The inner nodes are located on a set of straight spines; their distance to the upper wall is based on a proportionality factor that is constant throughout the deformation of the mesh. The coordinates of the inner nodes are thus related to the displacement of the interface nodes. 3.4. Determination of contact line The existence of three-phase static contact line, where the interface between the two fluids meets the wall, adds a further complication to the three-dimensional analysis. The location of fluid/fluid/solid contact line is not known a priori and must be included in the numerical problem [5,26,32]. In this work,

132

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

the extrapolation method for the determination of contact line is used for simplicity and the influence of extrapolation order can be reduced by the refined mesh around the contact line. This extrapolation technique is the simplest method to obtain the location of the contact line, but does not consider the exact wall effect. 3.5. Open boundary condition method In the finite element method, the boundary integral parts in Eqs. (16) and (22) are generally replaced by the boundary conditions, either essential or natural boundary conditions. The boundary condition should be compatible with the solution to be calculated within the domain of interest. If the accurate solution at the boundary is not known a priori, such as in the outflow boundary of the viscoelastic coextrusion process, it is almost impossible to impose a proper boundary condition. In the OBC method, boundary integral parts are treated as parts of unknown equations that have to be solved. Thus, the imposition of boundary condition along the outflow is not required at all. The components of the stiffness matrix related to the outflow boundary need to be modified to accommodate the boundary integral terms. Boundary nodal values are then obtained as a part of the solution in this method. OBC method keeps the solution to pass through the outflow boundary without undergoing significant distortion and influencing the interior solution [24,25].

4. Numerical results and discussion The die geometry and boundary conditions used in this study are shown in Fig. 1. Along the inlet boundary, the fully developed velocity and extra stress profiles are imposed; at the walls, no-slip condition is imposed; along the symmetric plane, zero traction condition is imposed. The flow rate ratio of two layers is set to 0.5 in this study. Two finite element meshes used in this study are given in Fig. 2. The two meshes are identical except the downstream length. The basic information on the meshes used in this simulation is summarized in Table 1, including the number of elements, nodes, and variables.

Fig. 1. Die geometry and the imposed boundary conditions.

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

133

Fig. 2. Finite element meshes used in this simulation: (a) Mesh 1; (b) Mesh 2.

Table 1 Basic information on the meshes used in the simulation

Number of elements Number of nodes Number of decoupling at cross-section Number of variables

Mesh 1

Mesh 2

672 6786 8×6 34008

960 9594 8×6 47892

As a simple criterion for the encapsulation phenomena, the degree of encapsulation (DE) is defined as DE =

yw − yc × 100 (%), 2L

(30)

where yw is the height of the interface at side wall, yc the height of the interface at symmetric plane (x = 0), and L the characteristic length. The interface profile near the wall with different extrapolation orders was studied by the authors [26] with the same mesh used here. The interface is nearly the same when the extrapolation order is greater than two. Thus, quadratic extrapolation method is used throughout the following work for simplicity. 4.1. Effect of the outlet boundary condition At the outlet boundary, we used the open boundary condition and compared the results with those from the fully developed boundary condition (FBC), where the normal traction force and the tangential velocities are set to zero. Figs. 3 and 4 show the contact line position along the downstream direction when the UCM and PTT models [27] are used. When the FBC was imposed at the outlet, the result from

134

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

Fig. 3. Contact line position along the downstream direction with different outlet boundary conditions and meshes using the UCM model (ηII /ηI = 2.5, λI = 0, λII = 0.1).

the UCM model showed an abrupt change of the contact line position near the outlet. The result from the PTT model also shows a deviation of the contact line position near the outlet and goes to a constant value immediately. In the result from the PTT model, the contact line position should not go to the steady state value because of the continuously developing contact line position along the downstream direction. The results from the OBC method show no abrupt change or deviation of contact line position near the outlet.

Fig. 4. Contact line position along the downstream direction with different outlet boundary conditions and meshes using the PTT model (η0,II /η0,I = 2.5, λI = λII = 0.1, ξ I = ξ II = 0.44, ε I = ε II = 1.5).

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

135

It confirms that the OBC method is more appropriate for the three-dimensional numerical simulation of viscoelastic coextrusion process. It seems that the problems in adopting the ambiguous outlet boundary condition raised by Takase et al. [18] could be resolved by using the OBC method. In order to examine the dependency of the solution on the downstream length, the results obtained with Mesh 1 and Mesh 2 using the OBC method are also compared in Figs. 3 and 4. The result shows that the solution is independent of the downstream length. 4.2. Upper convected Maxwell fluid Elasticity ratio has been considered as a minor factor in the encapsulation phenomena so far. Uneven and ragged interface shape is known to result from elasticity difference. The outlet interface profiles with different elasticity ratio are shown in Fig. 5, where the elasticity of fluid I (λI ) is set to 0.1 and the one of fluid II (λII ) is set to 0.1, 0.2, 0.4, 0.6, and 0.3, respectively. The viscosity of both fluids is the same. As the elasticity ratio (λII /λI ) increases from 1.0 to 8.0, the interface distortion seems to increase even though the effect is not so large and only a little change of outlet interface position from the height of the separating plate (y/L = 1.0) can be observed. Furthermore, the contact line position at the outlet is nearly the same when the elasticity ratio is greater than 1.0. In Fig. 6, the outlet interface profiles are shown when the elasticity of fluid II (λII ) varies from 0 to 0.8. Fluid I is Newtonian (λI = 0) and the viscosity ratio (ηII /ηI ) is set to 2.5. When the elasticity of fluid II (λII ) is set to 0, the interface position moves downward from the height of the separating plate (y/L = 1.0) and fluid I wraps around fluid II due to the viscosity difference. As the elasticity of fluid II (λII ) increases from 0 to 0.8, there is only a little change of interface position arising from the elasticity effect, and the interface distortion increases even though it is not so large. From Figs. 5 and 6, we can conclude that the viscosity difference between the polymer melts is the controlling factor of the interface position and the encapsulation phenomena at the outlet [1–4]. And

Fig. 5. Outlet interface profile with different elasticity ratios using the UCM model (ηII /ηI = 1.0).

136

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

Fig. 6. Outlet interface profile with different elasticity ratios using the UCM model (ηII /ηI = 2.5).

a significant difference in elasticity of the two components produces only a little change of interface position relatively to that observed with a similar viscosity ratio shift, and the melt elasticity difference causes the interface distortion even though it is not so large [3]. From Fig. 7, it can be observed that the transient overshoot of contact line position along the downstream direction becomes larger as the elasticity of fluid II (λII ) increases, but the contact line position at the

Fig. 7. Contact line position along the downstream direction with different elasticity of fluid II using the UCM model (ηII /ηI = 2.5).

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

137

outlet nearly coincides with each other in spite of different λII value. Karagiannis et al. [5] showed experimentally the gradual increase of interface curvature along the downstream direction, but we could not observe the gradual increase using the UCM model. As the second normal stress difference (N2 ) is zero, while the first normal stress difference (N1 ) is non-zero in the UCM model, it may be inferred that the first normal stress difference alone has little effect on the gradual increase of the degree of encapsulation along the downstream direction. 4.3. Phan-Thien and Tanner fluid The parameter ‘s’ used by Takase et al. [18] implies the addition of Newtonian shear viscosity component to the total shear viscosity. Thus, the resulting constitutive equation is different from the original PTT model [27]. Since the parameter ‘s’ is set to zero in this study, the original PTT model is recovered and this form is used throughout the simulation. We will show the dependency of interface shape and the degree of encapsulation on the parameters: η ratio, λ, ξ , and ε. Above all, the effect of the second normal stress difference (N2 ) will be shown clearly. 4.3.1. Effect of shear viscosity ratio and elasticity When the elasticity (λ) increases from 0.1 to 0.3, the shear viscosity curve also changes, as shown in Fig. 8. Because the maximum (dimensionless) shear rate calculated in the simulation is about 9.0, the (dimensionless) shear viscosity can change with the elasticity (λ) in this range of shear rate, and thus higher shear rate gives lower shear viscosity. As the average shear rate of fluid I with lower shear viscosity is greater than that of fluid II especially near the side wall, the shear viscosity ratio (ηII /ηI ) near the side wall goes larger than 2.5, which is the zero shear viscosity ratio (η0,II /η0,I ). Furthermore, the degree of decrease in shear viscosity increases at higher shear rate region as the elasticity (λ) increases, thus the shear viscosity ratio (ηII /ηI ) increases more as the elasticity (λ) increases, as represented in Fig. 9.

Fig. 8. Shear viscosity curves with different elasticity in the PTT model (η0,II /η0,I = 2.5, ξ I = ξ II = 0.44, ε I = ε II = 1.5).

138

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

Fig. 9. Shear viscosity ratios along the outlet interface with different elasticity in the PTT model (η0,II /η0,I = 2.5, ξ I = ξ II = 0.44, ε I = ε II = 1.5).

In Figs. 10 and 11, the outlet interface profile and the degree of encapsulation along the downstream direction are shown with increasing the elasticity (λ) from 0.1 to 0.3, while maintaining the elasticity ratio constant. As the elasticity of both layers increases, the interface curvature of the outlet interface as well as the degree of encapsulation along the downstream direction increases due to the increase of shear viscosity ratio, which can be induced from Figs. 8 and 9. In the simulation of coextrusion using the PTT

Fig. 10. Outlet interface profile with different elasticity using the PTT model (η0,II /η0,I = 2.5, ξ I = ξ II = 0.44, ε I = ε II = 1.5).

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

139

Fig. 11. Degree of encapsulation along the downstream direction with different elasticity using the PTT (η0,II /η0,I = 2.5, ξ I = ξ II = 0.44, ε I = ε II = 1.5).

model, higher elasticity results in higher shear viscosity ratio especially near the side wall (higher shear rate region), thus the effect of elasticity and shear viscosity ratio cannot be treated separately. In addition, the degree of encapsulation of the PTT fluid increases along the downstream direction unlike that of the UCM fluid (ξ I = ξ II = 0, ε I = ε II = 0), as shown in Fig. 11. Therefore, it will be natural to suspect that the second normal stress difference (N2 ) and/or the extensional viscosity have substantial effect on the increase of the degree of encapsulation along the downstream direction, because the second normal stress difference is non-zero and the extensional viscosity is finite in the PTT fluid. 4.3.2. Effect of the second normal stress difference The degree of encapsulation along the downstream direction is shown in Fig. 12, where the parameter ξ indicating the second normal stress difference property is set to 0, 0.22, and 0.44, respectively. When ξ is set to 0, the degree of encapsulation along the downstream direction reaches a steady state value. When ξ is set to 0.22 or 0.44, the degree of encapsulation along the downstream direction keeps on increasing in accordance with the experimental results of Karagiannis et al. [5], and the tendency increases as ξ increases. From these results, it may be concluded that the increase of the degree of encapsulation along the downstream direction is primarily due to the second normal stress difference, because the parameter ξ represents the magnitude of the second normal stress difference (N2 = −(ξ/2)N1 ). If the parameter ξ changes from 0.44 to 0.22, the shear viscosity increases a little only at high shear rate region, as shown in Fig. 13. When ξ value of fluid II (ξ II ) is fixed to 0.44 and that of fluid I (ξ I ) decreased from 0.44 to 0.22, the difference between the viscosity curves becomes smaller compared to the one with the same ξ values at higher shear rate region. So the shear viscosity ratio (ηII /ηI ) decreases, especially near the side wall (higher shear rate region), and in turn this may result in the decrease of the curvature of the outlet interface if we only consider the shear viscosity ratio effect. When ξ I is decreased from 0.44 to 0.22, the second normal stress difference of fluid I decreases to a half, and this will also affect

140

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

Fig. 12. Degree of encapsulation along the downstream direction with different ξ parameters using the PTT model (η0,II /η0,I = 2.5, λI = λII = 0.1, ε I = ε II = 1.5).

the interface shape. The curvature of the outlet interface increases as ξ I decreases from 0.44 to 0.22, as shown in Fig. 14. Therefore, it is concluded that the decrease of the second normal stress difference of fluid I while maintaining that of fluid II constant can increase the curvature of the outlet interface despite of the decreased shear viscosity ratio. Fluid II with the higher second normal stress difference in absolute value tends to protrude into fluid I with the lower second normal stress difference in absolute value, and in turn the increased degree of encapsulation can be observed. This result coincides with the theoretical results of White et al. [6] and Song [8].

Fig. 13. Shear viscosity curve with different ξ parameters in the PTT model (η0,II /η0,I = 2.5, λI = λII = 0.1, ε I = ε II = 1.5).

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

141

Fig. 14. Outlet interface profile with different ξ parameters of fluid I using the PTT model (η0,II /η0,I = 2.5, λI = λII = 0.1, ε I = ε II = 1.5).

4.3.3. Effect of extensional viscosity The ε value of 0, 0.1, and 1.5 corresponds to the infinite, extension thickening, and extension thinning extensional viscosity (ηE ), respectively. The degree of encapsulation along the downstream direction hardly changes in spite of the change of ε value while maintaining the ratio of ε constant, as shown in Fig. 15. From this result, it is confirmed that the extensional viscosity has little effect on the increase of

Fig. 15. Degree of encapsulation along the downstream direction with different ε parameters using the PTT model (η0,II /η0,I = 2.5, λI = λII = 0.1, ξ I = ξ II = 0.44).

142

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

Fig. 16. Outlet interface profile with different ε parameters of fluid I using the PTT model (η0,II /η0,I = 2.5, λI = λII = 0.1, ξ I = ξ II = 0.44).

the degree of encapsulation along the downstream direction, while the second normal stress difference has great influence on the increase of the degree of encapsulation along the downstream direction, as already shown in Fig. 12. When ε value of fluid II (εII ) is fixed to 1.5 (extension thinning) and that of fluid I (ε I ) decreased from 1.5 to 0.1 (extension thickening), the curvature of the interface at the outlet decreases a little bit, as shown in Fig. 16. This result can be explained with the decrease of the extensional viscosity ratio (ηE,II /ηE,I ), and this tendency coincides with the shear viscosity effect, but the sensitivity is very small. These results agree with the previous results of Takase et al. [18].

5. Conclusions In this work, we have carried out three-dimensional numerical simulation of viscoelastic coextrusion process in a square die using the OBC method, and examined the effect of viscoelastic properties on the interface distortion and the encapsulation phenomena. The problem of adopting the ambiguous outlet boundary condition in the three-dimensional numerical simulation of viscoelastic coextrusion process could be avoided by introducing the OBC method. The abrupt change or deviation of contact line position near the outlet that commonly occurs in using the fully developed outlet boundary condition could be clearly removed by using the OBC method. Using the upper convected Maxwell model as a constitutive equation, it was shown that the viscosity difference between the polymer melts is the controlling factor of the interface position and the encapsulation phenomena at the outlet. The interface distortion increases as the elasticity ratio increases, although it is not so large, and the first normal stress difference alone has little effect on the gradual increase of the degree of encapsulation along the downstream direction.

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

143

When the Phan-Thien and Tanner model is used, the effect of shear viscosity ratio and elasticity could not be treated separately in the flow simulation of viscoelastic fluids unlike the UCM fluid. When the elasticity of both layers increases maintaining the elasticity ratio constant, the degree of encapsulation increases resulting from the increased viscosity ratio especially near the side wall (higher shear rate region). The extensional viscosity has little effect on the degree of encapsulation along the downstream direction, but the ratio of extensional viscosities has some effect on the interface curvature. When the ratio of the second normal stress difference (ξ II /ξ I ) increases, the curvature of the outlet interface increases resulting from the increased protruding effect of fluid II with the higher second normal stress difference in absolute value into fluid I with the lower second normal stress difference. The second normal stress difference was found to have a great influence on the increase of the degree of encapsulation along the downstream direction. In this study, we used the single-mode PTT model, which overpredicts the shear sensitivity in the viscosity at high shear rates. Takase et al. [18] added the shear viscosity component to the total shear viscosity, and this addition results in a plateau region of shear viscosity at high shear rate region which cannot be observed in the normal polymer melts. Therefore, it will be necessary to expand the PTT model to the multimode to predict more realistic shear viscosity and the encapsulation phenomena observed in the coextrusion process. Acknowledgements The authors wish to acknowledge the support of Honam Petrochemical Corporation in Korea. This study was also supported by research grants from the Korea Science and Engineering Foundation through the Applied Rheology Center. References [1] C.D. Han, A study of bicomponent coextrusion of molten polymers, J. Appl. Polym. Sci. 17 (1973) 1289. [2] B.L. Lee, J.L. White, An experimental study of rheological properties of polymer melts in laminar shear flow and interface deformation and its mechanism in two-phase stratified flow, Trans. Soc. Rheol. 18 (1974) 467. [3] J.H. Southern, R.L. Ballman, Additional observations on stratified bicomponent flow of polymer melts in a tube, J. Polym. Sci. 13 (1975) 863. [4] C.D. Han, Y.W. Kim, Further observations of the interface shape of conjugate fibers, J. Appl. Polym. Sci. 20 (1976) 2609. [5] A. Karagiannis, A.N. Hrymak, J. Vlachopoulos, Three-dimensional studies on bicomponent extrusion, Rheol. Acta 29 (1990) 71. [6] J.L. White, R.C. Ufford, K.R. Dharod, R.L. Price, Experimental and theoretical study of the extrusion of two-phase molten polymer systems, J. Appl. Polym. Sci. 16 (1972) 1313. [7] A.A. Khan, C.D. Han, On the interface deformation in the stratified two-phase flow of viscoelastic fluids, Trans. Soc. Rheol. 20 (1976) 595. [8] J.K. Song, Numerical study on the interface movement in the sheath–core type coextrusion, Ph.D. Thesis, Seoul National University, 1993. [9] J. Dooley, R. Ramanathan, Coextrusion thickness variation — effects of geometry on layer rearrangement, in: Proceedings of the ANTEC’94 Conference XL, 1994, p. 89. [10] J. Dooley, K.S. Hyun, K. Hughes, An experimental study on the effect of polymer viscoelasticity on layer rearrangement on coextruded structures, Polym. Eng. Sci. 38 (1998) 1060. [11] A.G. Dodson, P. Townsend, K. Walters, Non-Newtonian flow in pipes of non-circular cross-section, Compt. Fluids 2 (1974) 317.

144

K.B. Sunwoo et al. / J. Non-Newtonian Fluid Mech. 99 (2001) 125–144

[12] P. Townsend, W.M. Waterhouse, K. Walters, Secondary flow in pipes of square cross-section and the measurement of the second normal stress difference, J. Non-Newtonian Fluid Mech. 1 (1976) 107. [13] B. Debbaut, T. Avalosse, J. Dooley, K. Hughes, On the development of secondary motions in straight channels induced by the second normal stress difference: experiments and simulations, J. Non-Newtonian Fluid Mech. 69 (1997) 255. [14] B. Debbaut, J. Dooley, Secondary motions in straight and tapered channels: experiments and three-dimensional finite element simulations with multimode differential viscoelastic model, J. Rheol. 43 (6) (1999) 1525. [15] H. Mavridis, A.N. Hrymak, J. Vlachopoulos, Finite-element simulation of stratified multiphase flows, AIChE J. 33 (1987) 410. [16] J.F. Agassant, A. Fortin, Prediction of stationary interfaces in coextrusion flows, Polym. Eng. Sci. 34 (1994) 1101. [17] W.A. Gifford, A three-dimensional analysis of coextrusion, Polym. Eng. Sci. 37 (1997) 315. [18] M. Takase, S. Kihara, K. Funatsu, Three-dimensional viscoelastic numerical analysis of the encapsulation phenomena in coextrusion, Rheol. Acta 37 (1998) 624. [19] D. Givoli, Numerical Methods for Problems in Infinite Domains, Elsevier, Amsterdam, 1992. [20] R.L. Sani, P.M. Gresho, Resume and remark on the open boundary condition minisymposium, Int. J. Num. Methods Fluids 18 (1994) 983. [21] M. Gupta, C.A. Hieber, K.K. Wang, Viscoelastic modelling of entrance flow using multimode Leonov model, Int. J. Num. Methods Fluids 24 (1997) 493. [22] P. Wapperom, M.A. Hulsen, J. Zanden, A numerical method for steady and nonisothermal viscoelastic fluid flow for high Deborah and Peclet numbers, Rheol. Acta 37 (1998) 73. [23] T.C. Papanastasiou, N. Malamataris, K. Ellwood, A new outflow boundary condition, Int. J. Num. Methods Fluids 14 (1992) 587. [24] S.J. Park, S.J. Lee, On the use of the open boundary condition method in the numerical simulation of nonisothermal viscoelastic flow, J. Non-Newtonian Fluid Mech. 87 (1999) 197. [25] S.J. Park, S.J. Lee, Numerical simulation of coextrusion process of viscoelastic fluids, Rheol. Acta. 13 (2001) 37. [26] K.B. Sunwoo, S.J. Park, S.J. Lee, K.H. Ahn, S.J. Lee, Three-dimensional numerical simulation of nonisothermal coextrusion process with generalized Newtonian fluids, Korea–Australia Rheol. J. 12 (2000) 165. [27] N. Phan-Thien, R.I. Tanner, A new constitutive equation derived from the network theory, J. Non-Newtonian Fluid Mech. 2 (1977) 353. [28] D. Rajagopalan, R.C. Armstrong, R.A. Brown, Finite element methods for calculation of steady viscoelastic flow using constitutive equations with a Newtonian viscosity, J. Non-Newtonian Fluid Mech. 36 (1990) 159. [29] F. Debae, V. Legat, M.I. Crochet, Practical evaluation of four mixed finite element methods for viscoelastic flow, J. Rheol. 38 (2) (1994) 421. [30] J.M. Marchal, M.J. Crochet, A new mixed finite elements for calculating viscoelastic flows, J. Non-Newtonian Fluid Mech. 26 (1987) 77. [31] S.J. Lee, S.J. Lee, Numerical prediction of three-dimensional extrudate swell, Korean J. Rheol. 2 (1990) 35. [32] A. Torres, N. Hrymak, J. Vlachopoulos, J. Dooley, B.T. Hilton, Boundary conditions for contact lines in coextrusion flows, Rheol. Acta 32 (1993) 5–13.