вход по аккаунту


Constitutive equations for modelling of polymeric materials and rubbers.

код для вставкиСкачать
Asia-Pac. J. Chem. Eng. 2010; 5: 732–742
Published online 13 October 2009 in Wiley Online Library
( 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
School of Mechanical and Aerospace Engineering, Queen’s University Belfast, Ashby Building, Stranmillis Road, Belfast, Northern Ireland, UK
Dow-Corning Corp., Midland, MI, USA
Dow-Corning Ltd, Barry, UK
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
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
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.
 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
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
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
There are various viscosity functions available to
describe the shear behaviour, the simplest of which is
the power-law model
η = 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
τ = 2ηD
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
τ = 2η(ID , IID , IIID )D
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 |
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
where η(IID ) is the viscosity function and Eqn (3) is the
GNF model. The model can be written in terms of the
shear rate γ̇
γ̇ = |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
Figure 1.
Asia-Pac. J. Chem. Eng. 2010; 5: 732–742
DOI: 10.1002/apj
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 =
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.
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]
τ +λ
= 2ηD
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
τ +λ
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
− (∇v)T · τ − τ · ∇v
where the substantial or material time derivative that
translates with the material particles is given by
+ v · ∇τ
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
+ τ · (∇v)T + ∇v · τ
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
and the lower-convected Maxwell (LCM) model
τ +λ
= 2ηD
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’
= 2η D + λτ
τ +λ
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
τ + λ = 2η D + λτ
Asia-Pac. J. Chem. Eng. 2010; 5: 732–742
DOI: 10.1002/apj
Asia-Pacific Journal of Chemical Engineering
and Schowalter[10] defined the Gordon and Schowalter
(GS) convected derivative, thus
and the Oldroyd-A fluid
τ + λ = 2η 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
where τp is the polymer stress and τs the solvent stress.
The polymer stress satisfies the UCM equation,
τp + λ
= 2ηp D
whilst the solvent stress satisfies the Newtonian relation
τs = 2ηs D
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
Apelian et al .[9] later developed the Modified UCM
model by making the relaxation time a function of trτ .
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:
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
τ +λ
The retardation time is related to the total shear viscosity
η = ηp + ηs and the relaxation time λr
τr =
ξ ξ∇
= 1−
2 τ
= 2ηD
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
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
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
τ2 = 2η2 D
If η2 ≥ 81 then the excessive shear-thinning is avoided
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
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 + λ
= 2η1 D
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
or the PTT(b) form
f (τ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
= 2η1 D
τ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
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
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
Table 1. Functions in Eqn. 35 for differential viscoelastic models.
Fluid model
f (τ1 )
No viscous
term, i.e.
τ2 = 0
UCM model if
η2 = 0
0<ξ ≤2
UCM model if
ξ =0
model if
0 < ξ ≤ 2 exp ελ
I + αλ
α = 0.5
λ(IID ), η1 (IID )
(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 + λ
= 2η1 D
τ2 = 2η2 D
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
τi + 2η2 D
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 )
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.
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.
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
Fluid memory can be found by integration of
Eqn (39) by parts and measure γ (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
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/λ
where λ is the relaxation time of the material. So
Eqn (39) gives
τ =
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) =
gi e −t/λi
i =1
the general linear viscoelastic model is obtained from
Eqn (43)
−∞ i =1
gi e −(t−t )/λi γ̇ (t ) dt τ=
−∞ 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
For an incompressible fluid I3 = 1 so
W = W (I1 , I2 )
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
The stress can now be expressed as
τ =
where ε is the Hencky strain
τ =2
e −(t−t )/λi
W =W
I1 = trC−1
1 −1 2 −2 trC
− trC
I2 =
= trC
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
τ1 =
m(t − t ) C−1 (t ) dt (47)
G trC−1
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
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
C −2
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
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
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 )
∂u(I1 , I2 , t − t )
C(t ) dt (56)
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 ) =
m(t − t )I1
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, t ) dt (58A)
τ12 =
N1 =
γ 2 (t, t ) dt (58B)
Therefore a kernel function is required
φ(γ , t − t ) ≡ 2
∂I1 ∂I2
which is found by performing step strain experiments
by imposing a shear strain γ at time zero to obtain the
shear stress given as
φ (γ , t − t ) dt = γ
φ (γ , s) ds
τ12 = γ
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
G(γ , t) = G(t)h(γ )
This implies that the kernel function also is factorable.
φ(γ , t − t ) = m(t − t )h(γ )
 2009 Curtin University of Technology and John Wiley & Sons, Ltd.
In simple shear flows, the separable KBKZ equation
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 )
The KBKZ equation takes on the general form
∂u(I1 , I2 ) −1 τ=
m(t − t ) 2
C (t )
∂u(I1 , I2 )
C(t ) dt (65)
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
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
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
B. CRAWFORD et al.
Asia-Pacific Journal of Chemical Engineering
from Eqn (51). The original Wagner equation used an
exponential function for h
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
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.
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 )
I = βI1 + (1 − β)I2
with the damping function
H (I ) = e −n(I −3)
So the advanced Wagner II equation became
gi −(t−t )/λi −n(I −3)1/2 −1 τ1 =
C (t ) dt −∞ i =1 λi
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
H (I ) =
α + (I − 3)
where α is a material parameter and I is the general
strain invariant of Eqn (71).
C−1 (t ) dt α
−∞ i =1
Again there are only two adjustable parameters, α
to fit the shear and β for the elongational behaviour.
τ1 =
This ratio had an important impact on swelling
behaviour. Thus the Luo-Tanner model is
t N
gi −(t−t )/λi
τ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.
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
[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
τ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
Asia-Pac. J. Chem. Eng. 2010; 5: 732–742
DOI: 10.1002/apj
Asia-Pacific Journal of Chemical Engineering
S 2 = I1 − 3
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.
The Weissenberg number describes the degree of elasticity in the flow and the importance of nonlinearities.
We = λc Vd /Rd = λc γw
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 .
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
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
B. CRAWFORD et al.
Asia-Pacific Journal of Chemical Engineering
equations and their use in computational dynamics for
two silicone systems.
[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,
[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,
[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;
[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,
[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,
[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,
[26] A. Kaye. CoA Note 134, College of Aeronautics: Cranfield,
[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,
[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,
[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,
[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,
[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
Без категории
Размер файла
167 Кб
constitution, equations, modellina, rubber, material, polymeric
Пожаловаться на содержимое документа