ASIA-PACIFIC JOURNAL OF CHEMICAL ENGINEERING Asia-Pac. J. Chem. Eng. 2010; 5: 732–742 Published online 13 October 2009 in Wiley Online Library (wileyonlinelibrary.com) DOI:10.1002/apj.400 Research Article Constitutive equations for modelling of polymeric materials and rubbers B. Crawford1 , J. K. Watterson1 , P. L. Spedding1 *, R. I. Gault1 , W. Herron2 and M. Proctor3 1 School of Mechanical and Aerospace Engineering, Queen’s University Belfast, Ashby Building, Stranmillis Road, Belfast, Northern Ireland, UK Dow-Corning Corp., Midland, MI, USA 3 Dow-Corning Ltd, Barry, UK 2 Received 7 January 2009; Revised 18 August 2009; Accepted 18 August 2009 ABSTRACT: Extrusion of polymeric materials can be made more efficient by modelling with computational techniques. Material-dependent data obtained through a constitutive equation is required to close the fluid dynamic equations. The paper outlines constitutive equations developed for non-Newtonian behaviour and details the choice that is required to obtain a constitutive equation that encapsulates actual material behaviour. 2009 Curtin University of Technology and John Wiley & Sons, Ltd. KEYWORDS: constitutive equations; non-Newtonian flow; differential viscoelastic models; integral viscoelastic models INTRODUCTION Production of polymeric materials often involves extrusion. The process is, however, only understood on a perproduct, trial and error basis, which is effective but inefficient. Therefore a more fundamental understanding of extrusion is needed using modelling and computational techniques. Generally a constitutive equation is required to close the system of equations in any fluid dynamics problem. The constitutive equation in question provides essential material-dependent information about the variables that appear in the constitutive equations. For Newtonian flow, the constitutive equation defines the viscosity and hence the relationship between deformation and stress. The absence of an appropriate constitutive equation for any system renders the problem ill-posed making it impossible to calculate the field variables. In non-Newtonian fluid mechanics there is no single constitutive equation that is suitable for all materials in all flow fields. There are many constitutive equations available to describe the nonlinear behaviour of polymeric materials.[1,2] For any situation the challenge is to choose an appropriate constitutive equation that encapsulates the actual behaviour of the material in question. The simplest constitutive equations are the generalized Newtonian fluid (GNF) models, which are only *Correspondence to: P. L. Spedding, School of Mechanical and Aerospace Engineering, Queen’s University Belfast, Ashby Building, Stranmillis Road, Belfast BT9 5AH, Northern Ireland, UK. E-mail: a.p.doherty@qub.ac.uk 2009 Curtin University of Technology and John Wiley & Sons, Ltd. capable of describing nonlinear viscous behaviour. They are used where there is a lack of rheological data or there is interest in shear-thinning phenomena. At the other end of the spectrum, sophisticated nonlinear viscoelastic models attempt to capture all of the observed nonlinear viscoelastic phenomena such as elastic recoil, i.e. fluid memory, normal direction shear stress, stress relaxation and elongational viscosity. There are two main categories of viscoelastic models namely differential and integral models. Differential models are simpler to integrate into numerical threedimensional schemes but are computationally expensive when simulating flows with a broad spectrum of relaxation times. This is because multiple modes are required to fit the rheological properties accurately. Integral models require an accurate calculation of the deformation history of the fluid particles to calculate the stress in the material, thus requiring complex particle tracking algorithms. Despite this limit, these models are more suited for describing a broad range of relaxation times. This is because with integral models most computer time is spent on particle tracking and little time is devoted to the stress calculation. Thus adding additional relaxation modes does not significantly increase computational cost. All constitutive equations must obey four formulation principles to be admissible for calculations. First is the principle of frame-invariance. The predicted behaviour of the material as expressed by the constitutive equation must not depend on the coordinate system used in the calculation. This is automatically satisfied if the constitutive equations are formulated using a consistent Asia-Pacific Journal of Chemical Engineering MODELLING OF POLYMERIC MATERIALS AND RUBBERS tensorial form. The relevant tensor notation is given in texts.[1,2] Second is the principle of material objectivity. All the constitutive equations must be independent of any super-imposed rigid body motions such as rotations, translations and any time-dependent motions. Third is the principle of determinism of stress and local action. It was proposed by Oldroyd[3] and states that the stress in a fluid element must be determined only by the rheological history of that element independently of neighbouring elements. Finally there is the principle of fading memory. Coleman and Noll[4] proposed that the influence of deformation on the stress decreases as time progresses after deformation. Thus in the limit for slow flows the material behaves like a viscous liquid, whilst for rapid flows it behaves like an elastic solid. Owing to the third principle, the concept of fading memory must be associated with fluid particles and not with points in space. Thus in the formulation of constitutive equations of whatever type the Eulerian notion of fixed points in space must be discarded. The main thrust of this and subsequent work will be directed to the modelling of silicone gum and silicone rubbers recognizing that predictive models that apply to these systems may not be suitable for other polymers. The most common constitutive equations will now be reviewed. There are various viscosity functions available to describe the shear behaviour, the simplest of which is the power-law model (n−1) η = K |IID | where K is the consistency index and n is the nonNewtonian power-law index. The relation is only useful for modelling flows where the typical shear rates lie in the power-law region of the shear viscosity curve. A more realistic model, capable of describing the low and high shear rate plateaus in the viscosity curve is the Carreau model[5] 2 (n−1)/2 η = η∞ + (η0 − η∞ ) 1 + λN |IID | The simplest constitutive equation is the Newtonian fluid τ = 2ηD (1) where τ is the extrastress tensor, D is the rate of deformation tensor and η is the viscosity constant which is a function of temperature. The Newtonian fluid cannot describe the shear-thinning behaviour of polymonic materials. However, by generalizing Eqn (1) allowance can be made for a nonlinear shear viscosity by use of the invariants of the rate of deformation tensor. (2) τ = 2η(ID , IID , IIID )D (6) where η0 is zero-shear viscosity, η∞ is infinite shear viscosity which is usually taken to be zero and λN is the natural time of the model which should not be confused with the characteristic relaxation time of the material. The value of λN is taken as the reciprocal of the shear rate at which the viscosity curve deviates from Newtonian to power-law behaviour. The Yasuda model[6] addresses the problem of possible variation in the curvature of the transition from Newtonian plateau to the power-law region by using a parameter ‘a’ in the Carreau model a (n−1)/a η = η∞ + (η0 − η∞ ) 1 + λN |IID | THE GNF (5) (7) Figure 1 presents the shape of the shear viscosity curves for these models. The power-law generalized Newtonian fluid (GNF) model grossly over-predicts the viscosity at low shear rates where many materials exhibit a Newtonian plateau. For incompressible materials in simple shear, ID and IIID are zero, respectively. So Eqn (2) reduces to τ = 2η(IID ) D (3) where η(IID ) is the viscosity function and Eqn (3) is the GNF model. The model can be written in terms of the shear rate γ̇ (4) γ̇ = |IID | 2009 Curtin University of Technology and John Wiley & Sons, Ltd. Shear viscosity curves for some common generalized Newtonian fluid (GNF) models. This figure is available in colour online at www.apjChemEng.com. Figure 1. Asia-Pac. J. Chem. Eng. 2010; 5: 732–742 DOI: 10.1002/apj 733 734 B. CRAWFORD et al. Asia-Pacific Journal of Chemical Engineering For purely viscous constitutive equations the extensional viscosity ηe is related to the shear viscosity η through the Trouton ratio Tr Tr = ηe =3 η (8) In general Eqn (8) does not hold for polymeric materials, which typically exhibit a much higher extensional viscosity. Linear polymers may show tension-thinning extensional behaviour, while branched polymers can show a tension-thickening response. Thus GNF models are usually extremely poor at modelling the extensional viscosity of polymeric materials. For simple axisymmetric extrusion, GNF models typically will under-estimate die swell and the entrance pressure drop and will fail to capture the normal stress differences and any vortex enhancement behaviour that may occur. Furthermore GNF models cannot handle time-dependent phenomena associated with memory effects. For these purposes a viscoelastic constitutive equation is required. DIFFERENTIAL VISCOELASTIC MODELS Linear differential viscoelastic model For differential viscoelastic models the extrastress tensor is related to the velocity gradient by means of the linear Maxwell equation [7] τ +λ dτ = 2ηD dt (9) where λ is the relaxation time of the material. The Maxwell model is the differential form of the general equation of linear viscoelasticity and so is only applicable in situations where the strain is small. As the Maxwell model is not frame invariant, it is not admissible as a constitutive equation. Quasi-linear differential viscoelastic models To make the Maxwell model frame invariant, it must be written in a more general form using a coordinate system that is convected or deformed with the material elements.[3] Thus the time derivative in Eqn (9) must be replaced, thus δτ = 2ηD τ +λ δt (10) with either the base coordinate vectors parallel to material lines or the base vectors normal to material planes. The first system is known as contravariant base 2009 Curtin University of Technology and John Wiley & Sons, Ltd. vectors and yield the upper-convected derivative, which can be written in terms of the stress as Dτ ∇ − (∇v)T · τ − τ · ∇v = Dt τ (11) where the substantial or material time derivative that translates with the material particles is given by Dτ ∂τ = + v · ∇τ Dt ∂t (12) In the second convected coordinate system, the base vectors are normal to material planes and these are termed covariant base vectors, which yield the lowerconvected derivative τ given by Dτ + τ · (∇v)T + ∇v · τ = Dt τ (13) Both the upper and lower-convected derivatives only allow for affine deformations of the polymer strands. Affine motion means that any strand end-to-end vector initially coincident with a macroscopic vector embedded in the continuum, remains parallel to and of equal length to the macroscopic vector during deformation. These two convected derivatives yield two new frame invariant quasi-linear Maxwell models, the upperconvected Maxwell (UCM) model τ +λ ∇ = 2ηD τ (14) and the lower-convected Maxwell (LCM) model τ +λ = 2ηD τ (15) The UCM and LCM models both have a constant viscosity and quadratic first normal stress difference. However, the LCM model gives an unrealistic value for the ratio of the second and first normal stress differences N2 /N1 and so is rarely used. The UCM model has been extensively used to test numerical algorithms but does exhibit an infinite elongational viscosity at finite extension rates and a zero value for N2 . The Oldroyd[3] proposal involved the Jeffreys’ equation[8] δτ δD = 2η D + λτ (16) τ +λ δt δt where λr is the retardation time. This equation appends an additional term to the Maxwell equation as a solvent contribution to the stress. Rendered frame invariant results in the Oldroyd-B fluid ∇ ∇ (17) τ + λ = 2η D + λτ D τ Asia-Pac. J. Chem. Eng. 2010; 5: 732–742 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering MODELLING OF POLYMERIC MATERIALS AND RUBBERS and Schowalter[10] defined the Gordon and Schowalter (GS) convected derivative, thus and the Oldroyd-A fluid ∇ τ + λ = 2η D + λτ τ D The Oldroyd-A fluid is disregarded for the same reasons as the LCM model. The Oldroyd-B fluid is equivalent to the following equation τ = τp + τs (19) where τp is the polymer stress and τs the solvent stress. The polymer stress satisfies the UCM equation, τp + λ ∇ = 2ηp D τp (20) whilst the solvent stress satisfies the Newtonian relation τs = 2ηs D ληs η (22) The Oldroyd-B fluid thus reduces to the UCM model when the retardation time is zero. Adding a purely viscous term to a constitutive equation has important implications for the mathematical behaviour of the model. Apelian et al .[9] later developed the Modified UCM model by making the relaxation time a function of trτ . NONLINEAR DIFFERENTIAL VISCOELASTIC MODELS The former quasi-linear models do not exhibit rheological properties typical of actual polymeric materials. This can only be achieved using fully nonlinear viscoelastic models where the assumption of affine motion is dropped. There are three main ways to develop nonlinear models. These are: (23) where the parameter ξ is the slippage parameter and has the range of values 0 ≤ ξ ≤ 2. The GS derivative introduces a nonaffine motion of the polymer strands that controls the rate at which stress builds up in the material. Equation (23) reduces to the upper-convected derivative when ξ = 0 and to the lower-convected derivative when ξ = 2. The corotational derivative is obtained when ξ = 1. Johnson and Segalman[11] proposed the use of the GS convected derivative in the UCM model rather than the upper-convected derivative yielding the JohnsonSegalman (JS) model τ +λ (21) The retardation time is related to the total shear viscosity η = ηp + ηs and the relaxation time λr τr = ξ ξ∇ = 1− + τ 2 τ 2τ (18) τ = 2ηD (24) This model has many advantages over the UCM model particularly in shear flows. It is able to predict a shearthinning viscosity, a nonquadratic first normal stress difference and a negative second normal stress difference. The ξ parameter controls the shear behaviour, while the ratio of the normal stresses is given by −ξ N2 = N1 2 (25) The JS model does not show any improvement in the prediction of extensional behaviour over the UCM model as an unbounded extensional viscosity still appears at finite extension rates. Also the JS model exhibits excessive shear-thinning with a maximum in the shear stress curve beyond which the shear stress decreases with increasing shear rate. Such physically unreasonable behaviour can be cured by adding a purely viscous term to the model with the extrastress tensor τ comprising a viscoelastic element τ1 and a purely viscous contribution τ2 τ = τ1 + τ2 (26) and Using nonaffine convected derivatives A nonlinear model can be developed by using a convected derivative that allows nonaffine motion of the polymer strands. As both this upper-convected derivative of Eqn (11) and the lower-convected derivative of Eqn (13) are admissible, then a linear combination of both must also be admissible. Using this method Gordon 2009 Curtin University of Technology and John Wiley & Sons, Ltd. τ1 + λ = 2η1 D τ1 τ2 = 2η2 D (27) (28) η If η2 ≥ 81 then the excessive shear-thinning is avoided 1 with the shear stress always increasing with shear rate. However, there is then a plateau in the shear viscosity curve as the shear rate approaches infinity. Asia-Pac. J. Chem. Eng. 2010; 5: 732–742 DOI: 10.1002/apj 735 736 B. CRAWFORD et al. Asia-Pacific Journal of Chemical Engineering Another model that uses the GS derivatives is the Phan-Thien-Tanner model (PTT),[12,13] which adds a functional dependence on the τ1 term. f (τ1 )τ1 + λ τ1 = 2η1 D (29) The function f (τ1 ) correlates the rate of creation and destruction of network junctions in either of two forms, the PTT(a) exponential model f (τ1 ) = exp ελ η1 trτ1 (30) trτ1 (31) or the PTT(b) form f (τ1 ) = 1 + ελ η1 For both forms of the equation, there are two additional fitting parameters ξτ and ετ that define the nonlinear behaviour. The first parameter ξτ controls the shearthinning behaviour, whilst ετ controls the extensional behaviour and blunts the elongational singularity that otherwise would be present. For the case when ετ = 0, the PTT model reduces to the JS model. Equation (30) predicts a maximum in the elongational viscosity before reaching a tension-thinning region at higher extension rates, while Eqn (31) predicts a tension-thickening region followed by a plateau at higher extension rates. Numerical and experimental investigations by Guillet et al .[14] on low density polyethylene Glomsaker et al .[15] on polyvinylchloride and Stachmuylders[16] on styrene-butadiene polymers showed that the PTT(a) Eqn (30) is more coherent with experimental data, giving a good fit for a wide variety of deformations and so is preferred in numerical simulations. The PTT(a) model does have drawbacks. It predicts spurious oscillations in the start up of shear flow when ξ = 0. These oscillations disappear when ξ = 0, but then N2 is predicted to be zero. It inherits the excessive shear-thinning behaviour shown by the JS model and so requires the addition of a purely viscous term to the extrastress tensor in the same manner as Eqns (26) and (28). Inclusion of nonlinear stress terms A nonlinear differential model also can be developed by using a nonlinear stress term rather than the linear stress term proposed by Oldroyd.[3] The Giesekus model [17] is an example that contains a term that is quadratic in stress ∇ αλ = 2η1 D (32) τ1 τ1 + λ I+ η1 τ1 2009 Curtin University of Technology and John Wiley & Sons, Ltd. The Giesekus model describes how the relaxation time of a molecule is altered when the surrounding molecules are orientated. The parameter α is termed the mobility parameter with values 0 ≤ α ≤ 1 and is associated with anisotropic drag on the constituent polymer molecules. When α = 0 the UCM model is obtained. The α parameter controls the shear thinning, the elongational viscosity and normal stress differences. At low shear rates α N2 =− N1 2 (33) When α > 0.5 the model may exhibit instabilities in shear flow and thus it is necessary to add a Newtonian component to the extrastress tensor in a similar way to that used for the JS and PTT models. The Giesekus model is able to describe accurately shearthinning behaviour and normal stress differences but generally does not perform well in the prediction of extensional flows. As it only has one additional fitting parameter, the Giesekus model is not as flexible as the PTT model which often leads to a tradeoff when fitting the shear and extensional behaviour. When α = 0.5 the Giesekus equation duplicates the model of Leonov[18] in viscometric flows. Brooks and Hughes[19] introduced streamline upwind (SU) methods particularly the Petrov-Galerkin (SUPG) method which introduce an anisotropic diffusivity acting only in the direction of the streamlines. King et al .[20] extended the approach with the Explicitly Elliptical Momentum Equation method, while Rajagopalan et al .[21] gave a more useful extension with the elastic viscous split stress (EVSS) method. Inclusion of invariants The simplest way of obtaining a nonlinear model is to make the quasi-linear model dependent upon the invariants of the deformation tensor in a manner similar to that used with the GNF model. This approach is used in the empirical White-Metzner (WM) model.[22] The UCM model is thus modified as τ1 + λ(IID ) ∇ = 2η1 (IID )D τ1 (34) where the function η1 (IID ) is selected on the basis of the shear viscosity curve and λ (IID ) chosen on the basis of the first normal stress difference curve. These functions commonly take the form of a powerlaw (Eqn (5)) or Carreau law (Eqn (6)). The WM model is relatively simple but while viscosity and the first normal stress difference curves are well handled the model still predicts an infinite elongational viscosity at finite extension rates and N2 to be zero. Barnes and Roberts[23] proposed a form for both η ( IID ) and λ Asia-Pac. J. Chem. Eng. 2010; 5: 732–742 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering MODELLING OF POLYMERIC MATERIALS AND RUBBERS Table 1. Functions in Eqn. 35 for differential viscoelastic models. Fluid model UCM Oldroyd-B JS PTT Giesekus WM f (τ1 ) I Comment No viscous term, i.e. τ2 = 0 0 I UCM model if η2 = 0 0<ξ ≤2 I UCM model if ξ =0 trτ I JS model if 0 < ξ ≤ 2 exp ελ 1 η1 ε =0 Leonov model τ 0 I + αλ 1 η1 when α = 0.5 0 I λ(IID ), η1 (IID ) ξ 0 (IID ) that enables the extensional viscosity to remain finite with the WM model. Table 1 summarises the more common nonlinear differential models in the form of Eqns (35a and b). f (τ1 )τ1 + λ τ1 = 2η1 D τ2 = 2η2 D (35a) (35b) Single mode differential models cannot accurately describe material behaviour over a wide range of shear rates. They, however, can be generalized in the form of a multimode model. Here the extrastress is taken as the sum of the Newtonian plus a number of individual mode stresses τi N τi + 2η2 D INTEGRAL VISCOELASTIC MODELS Linear integral viscoelastic model The basis of the integral viscoelastic model is Boltzmann’s superposition principle,[24] which states that the effects of mechanical history are linearly additive. That is, the influence of strains at different times in the past can be added linearly to give the stress at the present time. Consider the change in stress at the present time t due to a change in strain at the past time t . dτ = Gdγ (t ) (37) where G is the relaxation modulus dτ = G dγ (t ) dt = G γ̇ (t ) dt dt (38) The present shear stress is the integral of Eqn (38) over all past times from t = −∞ up to the present + t. t G(t − t )γ̇ (t ) dt (39) τ= −∞ Multimode differential viscoelastic models τ= For example with a multimode model with N relaxation modes in two-dimensional flow the number of unknown is 3N + 3. (36) i where N is the number of modes. Each τi is given by an equation such as Eqn (35a) most commonly using the GS or PTT expression for f (τi ). Thus each mode has a relaxation time λi , modulus gi and partial viscosity ηi . In terms of the Maxwell model of Eqn (9) this is equivalent to connecting N Maxwell elements in parallel. The values of gi and λi are obtained from the linear viscoelastic rheological behaviour. Small amplitude oscillatory shear deformations are used to obtain the discrete relaxation spectrum. Generally 1–1.5 modes per decade of experimental data of either shear rate or frequency are required to obtain a good fit to the rheological data. A constrained optimisation procedure is required to find λi , gi and ηi values that best fit rheological data. The main problem with multimode differential models is their significant increase in computer cost. 2009 Curtin University of Technology and John Wiley & Sons, Ltd. This is Boltzmann’s general equation of linear viscoelasticity. Thus the stress is an integral, over all past time, of the relaxation modulus (which is a material property times the deformation rate which is a flow property). The Eqn (39) enables information to be extracted on the stress relaxation modulus by measuring the stress and the strain in a test with a specific time pattern. This idea is central to all linear viscoelastic rheological tests where the type of deformation and the time pattern are not important provided the strain is small. Fluid memory can be found by integration of Eqn (39) by parts and measure γ (t ) t τ= m(t − t )γ (t ) dt (40) −∞ with the memory function of linear viscoelasticity being defined by m(t − t ) = − dG(t − t ) dt (41) The memory function is positive as the relaxation modulus G(t) decreases with time such that G∞ → 0. Viscoelastic fluids have a fading memory such that deformations experienced by a fluid element in the recent past contribute more to the current stress Asia-Pac. J. Chem. Eng. 2010; 5: 732–742 DOI: 10.1002/apj 737 738 B. CRAWFORD et al. Asia-Pacific Journal of Chemical Engineering than earlier deformations. In experiments, the stress relaxation modulus shows an exponentially decaying behaviour so can be modelled as G(t) = G0 e −t/λ (42) where λ is the relaxation time of the material. So Eqn (39) gives t τ = G0 e −(t−t )/λ γ̇ (t ) dt (43) −∞ Equation (43) is the integral form of the Maxwell model Eqn (9). Just as was the case with the differential models, a single mode does not give an acceptable model of the data. To improve the situation a number N of relaxation times are used so G(t) can be expressed as the sum of exponentially decaying functions thus G(t) = N gi e −t/λi (44) i =1 the general linear viscoelastic model is obtained from Eqn (43) τ= t N −∞ i =1 gi e −(t−t )/λi γ̇ (t ) dt τ= t N −∞ i =1 gi e −(t−t )/λi 2D(t ) dt The most realistic models of the integral type belong to the KBKZ familyof relations developed by Kaye[26] and Bernstein et al .[27] From classical rubber elasticity theory, the strain energy W is W = where C−1 is the Finger-Cauchy deformation tensor. Generalizing the dependence of W on the strain gives The extrastress tensor τ will require to be decomposed into a viscoelastic and purely viscous contribution, in a similar fashion to Eqns (26) and (28). The memory function m (t − t ) is expressed as a sum of exponential decays thus i =1 C−1 (50) (51) For an incompressible fluid I3 = 1 so W = W (I1 , I2 ) (48) 2009 Curtin University of Technology and John Wiley & Sons, Ltd. (52) The stress can now be expressed as τ = ∂W ∂ε (53) where ε is the Hencky strain τ =2 e −(t−t )/λi W =W (46) −∞ λi (49) I1 = trC−1 1 −1 2 −2 trC − trC I2 = = trC 2 I3 = det C−1 An admissible integral viscoelastic constitutive equation can be obtained by using the Finger–Cauchy finite deformation tensors.[1,2] . This is analogous to the use of convected derivatives in the development of admissible differential models. There exists an integral equivalent to the UCM model of Eqn (14) called the Lodge’s rubber-like liquid [25] given by t τ1 = m(t − t ) C−1 (t ) dt (47) N gi 1 G trC−1 2 The strain energy W is a scalar quantity and must be made independent of the frame of reference by being written as a function of the invariants of C−1 thus Quasi-linear integral viscoelastic models m(t − t ) = Nonlinear integral viscoelastic models (45) The eqn can be extended to three-dimensions thus The Lodge model reduces to the general linear viscoelastic model in the limit of small displacement gradients. At high strains in the nonlinear region the Lodge model does not predict the correct material behaviour by giving a similar response to the UCM model. The Oldroyd differential models also have an equivalent integral form. ∂W −1 ∂W C −2 C ∂I1 ∂I2 (54) For viscoelastic fluids the elastic energy and stress are allowed to relax and W can be assumed to be the historic integral of a function of I1 , I2 and t − t in the form t u(I1 , I2 , t − t ) dt (55) W = −∞ Asia-Pac. J. Chem. Eng. 2010; 5: 732–742 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering MODELLING OF POLYMERIC MATERIALS AND RUBBERS with the functions u(I1 , I2 , t − t ) being called the potential function. The corresponding stress tensor is given by t ∂u(I1 , I2 , t − t ) −1 C (t ) 2 τ= ∂I1 −∞ ∂u(I1 , I2 , t − t ) −2 C(t ) dt (56) ∂I2 Equation (56) represents the KBKZ class of constitutive equations. The Lodge equation is a special case of Eqn (56) obtained when the potential function is u(I1 , I2 , t − t ) = 1 m(t − t )I1 2 (57) The Eqn (56) is difficult to handle as the potential function u(I1 , I2 , t − t ) must be determined and it depends on two strain invariants and time. Simplification can be achieved by considering the shearing deformations where the shear stress and first normal stress differences are given by t ∂u ∂u 2 − γ (t, t ) dt (58A) τ12 = ∂I ∂I 1 2 −∞ t ∂u ∂u N1 = 2 − γ 2 (t, t ) dt (58B) ∂I ∂I 1 2 −∞ Therefore a kernel function is required ∂u ∂u φ(γ , t − t ) ≡ 2 − ∂I1 ∂I2 (59) which is found by performing step strain experiments by imposing a shear strain γ at time zero to obtain the shear stress given as 0 ∞ φ (γ , t − t ) dt = γ φ (γ , s) ds τ12 = γ −∞ t (60) where s is the elapsed time. The function Ø (γ , s) can be found by taking the time derivative of the relaxing shear stress. By knowing Ø (γ , s), τ12 and N1 can be predicted by the KBKZ equation for any shearing deformation history. The KBKZ equation can be further simplified by introducing the concept of time-strain separability. In nonlinear step strain experiments performed at different strain levels plots of G (γ , t) against t are all parallel suggesting that the stress relation modulus is factorable or separable into time-dependent and strain-dependent terms. G(γ , t) = G(t)h(γ ) (61) This implies that the kernel function also is factorable. φ(γ , t − t ) = m(t − t )h(γ ) (62) 2009 Curtin University of Technology and John Wiley & Sons, Ltd. is In simple shear flows, the separable KBKZ equation t m(t − t )h(γ )γ (t, t ) dt (63) τ12 = −∞ where h(γ ) ≤ 1 is the damping function. At small strains as γ → 0, G(γ , t) must reduce to the linear relaxation modulus G(t). It is possible to find h(γ ) graphically by establishing the amount of vertical shift required to superpose each nonlinear relaxation curve onto the linear relaxation curve. Separability has been shown to hold for deformations other than shear, so the potential function can be factorised as u(I1 , I2 , t − t ) = m(t − t )u(I1 , I2 ) (64) The KBKZ equation takes on the general form t ∂u(I1 , I2 ) −1 τ= m(t − t ) 2 C (t ) ∂I1 −∞ ∂u(I1 , I2 ) −2 C(t ) dt (65) ∂I2 known as the separable or factorable KBKZ equation. In Eqn (65) time effects described through the memory function are separated from the strain effects associated with I1 and I2 . Simple linear viscoelastic experiments allow the memory function m(t − t ) to be determined leaving only nonlinear experimental data to find the u(I1 , I2 ) term. However, the u(I1 , I2 ) term depends implicitly on time through the invariant terms and the Finger–Cauchy tensors. Equation (65) is a special case of the more general superposition integral equation of Rivlin and Sawyers.[28] t τ= m(t − t ) φ1 (I1 , I2 ) C−1 (t ) −∞ +φ2 (I1 , I2 ) C (t ) dt (66) In this equation the functions Ø1 and Ø2 are not necessarily derivatives of a potential function as was the case with Eqn (59). In deformations other than shear the KBKZ equations require extensive experimental data to describe the arbitrary deformation histories. Wagner[29] proposed a simplification that involves dropping the Cauchy term from Eqn (66) to give t m(t − t )h(I1 , I2 ) C−1 (t ) dt (67) τ1 = −∞ The memory function is used in the form of Eqn (48) and the damping function h depends on the first and second invariants of the Finger deformation tensor, C−1 , Asia-Pac. J. Chem. Eng. 2010; 5: 732–742 DOI: 10.1002/apj 739 740 B. CRAWFORD et al. Asia-Pacific Journal of Chemical Engineering from Eqn (51). The original Wagner equation used an exponential function for h (68) h(I2 ) = exp −n(I2 − 3)1/2 where n is a material coefficient. However by omitting the Cauchy term eqn (67) predicts N2 = 0 in shear. Also by using a reversible damping function an incorrect response is given in reversing deformations as well as exhibiting excessive elastic recoil in start-up experiments of uniaxial elongation. Wagner[30] corrected these deficiencies by using an irreversible damping function H(I1 , I2 ) to give t m(t − t ) H (I1 , I2 ) C−1 (t ) dt (69) τ1 = Unfortunately like the Wagner equations N2 is predicted as zero in purely viscometric flows. Luo and Tanner[33] proposed improvements to the PSM model. They suggested to use a spectrum of β values, associated with each relaxation mode rather than a single β value, thus improving the prediction in extensional flows. Secondly, the Finger tensor C−1 was replaced by (1)/(1 − θ ))[C−1 (t ) + θ C(t )] where θ is a parameter controlling the ratio of normal stress differences. θ N2 = N1 1−θ −∞ The damping function H was made irreversible by the following. As the strain changed between time t and t the lowest value of h that existed over the time interval was selected such that the damping function could only decrease. H (t, t ) = min h(t, t ) t=t t=t (70) I = βI1 + (1 − β)I2 (71) with the damping function H (I ) = e −n(I −3) 1/2 (72) So the advanced Wagner II equation became N gi −(t−t )/λi −n(I −3)1/2 −1 τ1 = e e C (t ) dt −∞ i =1 λi (73) With this equation, apart from the linear viscoelastic data only two adjustable parameters are required to fit the shear (n) and elongational (β) behaviour. The PSM model [32] used the general form of the Wagner Eqn (69) with a rational form of the damping function α (74) H (I ) = α + (I − 3) t where α is a material parameter and I is the general strain invariant of Eqn (71). t N gi α C−1 (t ) dt α + (I − 3) −∞ i =1 (75) Again there are only two adjustable parameters, α to fit the shear and β for the elongational behaviour. τ1 = λi This ratio had an important impact on swelling behaviour. Thus the Luo-Tanner model is t N θ gi −(t−t )/λi e τ1 = 1 − θ −∞ i =1 λi × α [C−1 (t ) (α − 3) + βi I1 + (1 − βi )I2 + θ C (t )] dt To better reconcile data in shear and extension, Wagner et al .[31] used a generalized strain invariant e −(t−t )/λi 2009 Curtin University of Technology and John Wiley & Sons, Ltd. (76) (77) Goublomme and Crochet [34] modified the Wagner Eqn (69) to incorporate the effect of the second normal stress difference. t N θ gi −(t−t )/λi −n(I −3)1/2 −1 e e [C (t ) τ1 = 1 − θ −∞ i =1 λi + θ C(t )] dt (78) In the same manner, they included the effect of the second normal stress differences in the PSM model. t N θ gi −(t−t )/λi e τ1 = 1 − θ −∞ i =1 λi × α [C−1 (t ) (α − 3) + βI1 + (1 − β)I2 + θ C(t )] dt (79) Equation (79) is identical to Eqn (77), but only permits a single value of the β parameter. One of the main drawbacks with the PSM damping function is that it is unable to simultaneously fit a strain-hardening response in planar elongation and a shear-thinning viscosity. For example Olley et al .[35,36] reported that the PSM model predicted diminishing vortex behaviour at increasing flow rates contrary to observation. They proposed a damping function h(S ) = α α + βS + (1 − β)S 2 (80) Asia-Pac. J. Chem. Eng. 2010; 5: 732–742 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering where S 2 = I1 − 3 MODELLING OF POLYMERIC MATERIALS AND RUBBERS (81) The damping function has proved to be useful for planar strain-hardening materials flowing in geometries where there is significant elongational flow. It can predict vortex growth in planar contraction flows but tended to over-predict N1 values. NUMERICAL EFFECTS The Weissenberg number describes the degree of elasticity in the flow and the importance of nonlinearities. We = λc Vd /Rd = λc γw (82) With λc as the material characteristic relaxation time, Vd is the mean flow velocity, Rd tube radius and γw developed wall shear rate. Elastic effects increase with We making numerical simulations highly nonlinear and increasingly difficult to solve until a critical limit Wec is reached above which simulation is impossible. Keunings[37,38] outlined the difficulties caused by Wec but Lipscomb et al .[39] suggested that problem could be alleviated by both relaxing the nonslip boundary condition and the viscoelastic contribution in the vicinity of singularities. Later workers verified these measures[40,41] . Further differential viscoelastic models exhibited a hyperbolic character that leads to instabilities at low We values leading Marchal and Crochet[42] to use 4 × 4 subelements for the viscoelastic stresses in conjunction with SU methods to delay the appearance of wiggling in the stress fields. The method however proved to be computationally expensive when extended to multimode calculations. Debae et al .[43] comprehensively evaluated these types of simulations.[19 – 21] The SU 4 × 4 method gave good simulations but was computationally expensive. In contrast, the EVSS[21] simulation was less demanding for the same accuracy of results. With integration simulations Debae[43] reported that the Galerkin method was unable to achieve high We values but otherwise its results were accurate. The SUPG method by contrast was unpredictable and failed to handle singularities. The SU 4 × 4 while being stable required mesh refinement studies to validate its performance. Thus from a numerical point of view emphasis remains in formulating/choosing simulations which can deliver mesh convergent results at realistic higher values of We . CONCLUSION As all GNF models are inexpensive in computer time they are used to simulate the purely viscous behaviour 2009 Curtin University of Technology and John Wiley & Sons, Ltd. of the material. The Yasuda model of Eqn (7) usually gave the best fit to the material behaviour offering more flexibility than the Carreau model of Eqn (6) in the transitional region between Newtonian and powerlaw behaviour. This region is especially important when modelling polymeric materials as pointed out by Swallow.[44] Of the differential viscoelastic models, the quasilinear type models prove to be special cases of the more general nonlinear models, the PTT models, the Giesekus model and the WM model. Of the PTT models the PTT(a) equations of Eqns (29) and (30) predict a maximum in the elongational viscosity before reaching a tension-thinning region at higher extension rates that was more coherent with experimental data than the PTT(b) form of Eqns (29) and (31) and therefore the PTT(a) equations were more suitable when modelling polymeric materials. The PTT(a) model has two fitting parameters ξ and ε controlling the shearthinning and extensional behaviour, respectively. The Giesekus model of Eqn (32) had the advantage of having only one fitting parameter α that controlled the rheological behaviour of the material. The model can accurately describe shear-thinning behaviour and normal stress differences, but cannot fit simultaneously the shear-thinning shear viscosity and tension-thinning extensional viscosity. The WM model of Eqn (34) is a relatively simple relation that can describe accurately the shear behaviour over a wide range of shear rates, but exhibits strange behaviour in predicting free surface deformations and die swell and an infinite uniaxial viscosity at finite extension rates. On balance, the best performance appears to be with the PTT(a) model. The quazi-linear integral viscoelastic models do not predict the correct material behaviour and gave a response similar to the UCM equation. With the nonlinear integral viscoelastic models the KBKZ Eqn (56) proved to be difficult to handle without simplification. That by Wagner as Eqns (69) and (73) only require two adjustable parameters n to fit the shear and β to fit the elongational behaviour. Subsequently the use of an additional damping function led to the PSM model of Eqn (75), which seems to be the preferred general form of the KBKZ models, as it is able to handle die swell characteristics. The Luo–Tanner modification of Eqn (77) employed a spectrum of β values to improve the prediction of extensional flows and the normal stress difference ratio. The Goublomme–Crocket modification of Eqn (79) allows prediction of the second normal stress difference. All these models require the use of fitting parameters to obtain acceptable prediction of material behaviour. These parameters must first be determined for each material by detailed experimentation. The resulting constitutive equations can then be used to predict the exact result expected from the material in question. Later work will follow the determination of the constitutive Asia-Pac. J. Chem. Eng. 2010; 5: 732–742 DOI: 10.1002/apj 741 742 B. CRAWFORD et al. Asia-Pacific Journal of Chemical Engineering equations and their use in computational dynamics for two silicone systems. REFERENCES [1] R.B. Bird, R.C. Armstrong, O. Hassager. Dynamics of Polymeric Liquids, Vol. 1, 2nd edn. Wiley: New York, 1987. [2] R.G. Larson. Constitutive Equations for Polymer Melts and Solutions, Butterworths: Boston, 1988. [3] J.G. Oldroyd. Proc. Roy. Soc., 1950; A200, 523. [4] B.D. Coleman, W. Noll. Rev. Mod. Phys., 1961; 33, 239. [5] P.J. Carreau, PhD Thesis, Department of Chemical Engineering, University of Wisconsin: 1968. [6] K. Yasuda, PhD Thesis, Massachusetts Institute of Technology, Cambridge, 1979. [7] J.C. Maxwell. Phil. Trans. Roy. Soc., 1867; A157, 49. [8] H. Jeffreys, The Earth, Cambridge University Press: London, 1929. [9] M.R. Apelian, R.C. Armstrong, R.A. Brown. J. NonNewtonian Fluid. Mech., 1988; 27, 299–321. [10] R.J. Gordon, W.R. Schowalter. Trans. Soc. Rheol., 1972; 16, 79. [11] M.W. Johnson, D. Segalman Jr. J. Non-Newtonian Fluid Mech., 1977; 2, 225. [12] N. Phan-Thien, R.I. Tanner. J. Non-Newtonian Fluid Mech., 1977; 2, 353. [13] N. Phan-Thien. J. Rheology, 1978; 22, 259–283. [14] J. Guillet. In Rheology of Polymer Melt Processing (Eds.: J.M. Piau, J.F. Agassant), Elsevier Science: 1996; pp.285–336. [15] T. Glomsaker, E.E. Hinrichsen, F. Irgens, P. Thorsteinsen. Rheol. Acto., 2000; 39, 80–96. [16] E. Slachmuylders, Polymer Flow Users Meeting, Brussels, Belgium, 1997. [17] H. Giesekus. J. Non-Newtonian Fluid Mech., 1982; 11, 69–109. [18] A.I. Leonov. Rheol. Acta, 1976; 15, 85–98. [19] A.N. Brooks, T.J.R. Hughes. Comp. Meth. Appl. Mech. Eng., 1982; 32, 199–259. [20] R.C. King, M.R. Apelian, R.C. Armstrong, R.A. Brown. J. Non-Newtonian Fluid Mech., 1988; 29, 147–216. [21] D. Rajagopalan, R.C. Armstrong, R.A. Brown. J. NonNewtonian Fluid Mech., 1990; 36, 159–192. 2009 Curtin University of Technology and John Wiley & Sons, Ltd. [22] J.L. White, A.B. Metzner. J. Appl. Polym. Sci., 1963; 7, 1867–1889. [23] H.A. Barness, G.P. Roberts. J. Non-Newtonian Fluid Mech., 1992; 44, 113–126. [24] L. Boltzmann. Ann. Phys. Chem., 1876; 7, 624. [25] A.S. Lodge. Elastic Liquids, Academic Press: New York, 1964. [26] A. Kaye. CoA Note 134, College of Aeronautics: Cranfield, 1962. [27] B. Bernstein, E.A. Kearsley, L.J. Zapas. Trans. Soc. Rheol., 1963; 7, 391–410. [28] R.S. Rivlin, K.N. Sawyers. Ann. Rev. Fluid Mech., 1971; 3, 117. [29] M.H. Wagner. Rheol. Acta, 1976; 15, 136. [30] M.H. Wagner. J. Non-Newtonian Fluid Mech., 1978; 4, 39. [31] M.H. Wagner, T. Raible, J. Meissner. Rheol. Acta, 1979; 18, 427–428. [32] A.C. Papanastasiou, L.E. Scriven, C.W. Macosko. J. Rheology, 1983; 27, 387–410. [33] X.L. Luo, R.I. Tanner. Int. J. Num. Meth. Eng., 1988; 25, 9–22. [34] A. Goublomme, M.J. Corchet. J. Non-Newtonian Fluid Mech., 1993; 47, 281. [35] P. Olley, R. Spares, P.D. Coates. J. Non-Newtonian Fluid Mech., 1999; 86, 337–357. [36] P. Olley. J. Non-Newtonian Fluid Mech., 2000; 95, 35–53. [37] R. Keunings. J. Non-Newtonian Fluid Mech., 1986; 20, 209–226. [38] R. Keunings. In Computer Modelling for Polymer processing (Ed.: C.L. Tucker), Oxford, 1989; pp.403–459. [39] G.G. Lipscomb, R.R. Keunings, M.M. Denn. J. NonNewtonian Fluid Mech., 1987; 24, 85–86. [40] M.R. Apelian, R.C. Armstrong, R.A. Brown, R. Keunings. J. Non-Newtonian Fluid Mech., 1988; 27, 299–321. [41] M.J. Crochet, R. Keunings. J. Non-Newtonian Fluid Mech., 1982; 10, 339–356. [42] J.M. Marchal, M.J. Crochet. J. Non-Newtonian Fluid Mech., 1987; 26, 77–114. [43] F. Debae, V. Legat, M.J. Crochet. J. Non-Newtonian Fluid Mech., 1994; 38, 421–444. [44] F.E. Swallow, Dow-Corning Technical Report TIS 1998, 1998; p.45778. Asia-Pac. J. Chem. Eng. 2010; 5: 732–742 DOI: 10.1002/apj

1/--страниц