Modeling the viscoelastoplastic behavior of waxy crude
Received January
Modeling the viscoelastoplastic behavior of waxy crude
Teng Houxing 0
Zhang Jinjun 0
0 National Engineering Laboratory for Pipeline Safety/Beijing Key Laboratory of Oil and Gas Distribution Technology , China
satisfactory.

content of 21.5% and 26.3%, respectively. In the past
crystals interlock with each other to form a threedimensional
spongelike network structure, completing the transition
of crude oil from sol to colloidal gel (Visintin et al, 2005;
beginning of deformation (Chang et al, 1998). As the shear
strain increases, the three dimensional network composed of
loose clusters is disrupted into a small number of relatively
large flocs or aggregates, losing the connectivity of the
network structure (Chang et al, 1998; Visintin et al, 2005). In
this process, the viscoelasticity decays, and the mechanical
response gradually transfers from being dominated by the
elastic effect to being dominated by the viscous effect. After
shearing under some high shear rate for a period of time, the
internal network structure is nearly totally destroyed, thus the
viscoelastic behavior is very weak and even can be ignored,
et al, 2002; Zhang et al, 2010). The viscoplastic models
stress is greater than the yield stress. For this reason, they
and cannot account for the initial viscoelastic effect before
the yield point (Labanda and Llorens, 2006). In recent years,
some scholars attempted to describe both the viscoelastic and
(Barnes, 1997; Mujumdar et al, 2002; Dullaert and Mewis,
Jiang (1998) developed a constitutive equation based on the
principle of a mechanical analogy composed of a spring, a
However, the elastic stress of the model increases with the
shear strain unbounded in time. To overcome this problem, a
the elastic stress in the modified model increases from 0 to
Mujumdar et al (2002) proposed a nonlinear rheological
a rheological model that can fully describe the
viscoelastofor the timedependent elastic, viscous, and yielding
phenomena. The total shear stress is separated into elastic
and viscous components which originate from intrafloc
relative importance of the elastic stress to viscous stress is
determined by the structural parameter. Besides, to control
the transformation from viscoelastic behavior to yielding
behavior, an upper limit of elastic strain is adopted and its
relation. Following the line of the Mujumdar model, Huang
the timedependent nonlinear viscoelastic behavior of
LDPE melt. The elastic stresses at different shear strains are
depicted by a nonlinear damping function h( ) along with a
structural parameter , that is ( , ) 0 ( ) . However,
the model does not introduce the concept of elastic strain,
bringing about the drawback that the damping function will
decrease in the case of the internal structure building up. As
to the Mujumdar model and the Huang model, the viscous
stress under constant shear rate increases with a decrease
shear rate for a period of time, the internal structure is nearly
totally damaged and the elastic stress will be very small,
then the total shear stress is mainly composed of the viscous
stress. In this condition, the total shear stress predicted by the
Dullaert and Mewis (2006) developed a general structural
on inelastic suspending fluid, whose shear stress is divided
into a particle and a medium contribution, and the particle
contribution is further subdivided into an elastic and a viscous
hydrodynamic contribution. As for the elastic contribution, a
limit or critical elastic strain is also adopted and the evolution
of elastic strain is described by a differential kinetic equation.
Zhu and Smay (2011) proposed an engineering model to
as a simple algebraic equation. At the very beginning of
deformation, the shear stress predicted by the two models
(Dullaert and Mewis, 2006; Zhu and Smay, 2011) has elastic
and viscous components; while for the two models (Mujumdar
et al, 2002; Huang and Lu, 2005), the shear stress only has the
elastic component without the viscous component. This is one
distinction between them. Nevertheless, many papers showed
Burghelea, 2009), the initial mechanical response is elastic.
From this aspect, the initial mechanical response predicted
by the Mujumdar and Huang model seems more realistic and
plausible for these materials.
In this work, we proposed a viscoelastoplastic model based
elastic and viscous components. As to the evolution of the
equation, following the line of the model proposed by Zhu
and Smay (2011). The proposed model is validated by
stepwise shear rate tests and hysteresis loop tests of Daqing
applicability of the model is further checked by predicting the
parameters determined from the stepwise shear rate test.
2 Proposed model
Just like the common practice for the structural kinetics
models (Mewis and Wagner, 2009), the structuring level
of the microstructure of a material is represented by a
nonnegative scaled structural parameter . It varies between the
values of 0 for the completely brokendown structure and
1 for the fully developed structure. The proposed model is
composed of an equation of state and a kinetic equation for
the structural parameter.
2.1 Equation of state
As the model proposed by Mujumdar et al (2002), the
total shear stress is also composed of an elastic stress e
and a viscous stress v, that is = e+ v. The elastic stress e
is assumed to obey the Hookean elastic response, that is e =
G( ) e, where G( ) is the structuredependent shear modulus
and e is the elastic strain caused by the presence of the
deformable network structure. When the network structure
is complete ( =1), the mechanical response is linear elastic,
without a viscous response. With an increase in the total
shear strain , the network structure is disrupted into smaller
structures or aggregates, resulting in the decay of the shear
modulus G and the rise of the viscous behavior. Most of the
Dullaert and Mewis, 2006; Zhu and Smay, 2011) assume
that the shear modulus G is dependent only on the structural
parameter . Here, we also assume that the shear modulus G
is proportional to instantaneous value of , i.e. ( ) 0 ,
where G0 is the shear modulus of the completely structured
material (=1).
As to the elastic strain e, the model proposed by
Doraiswamy et al (1991) assumes that e will increase linearly
c equal to the critical
strain at the yield point, and it will remain constant as long
as the deformation process continues in the same direction.
For the two models (Mujumdar et al, 2002; Dullaert and
Mewis, 2006), prior to yielding c is the same as the elastic
strain in the Doraiswamy model; after yielding the Mujumdar
c is assumed to be a function of the
c co , while for the Dullaert
model, the evolution of e with shearing conditions and time
is described by a specially proposed differential equation.
The elastic strain e in the Doraiswamy and Mujumdar model
is continuous, but has an inflection point at the yield point.
Consequently, the elastic stress e is also continuous but has
e and e
should be continuous and smooth. As for the model proposed
the same as that of the Doraiswamy model, but to overcome
e
e
0
0 c 1 exp(
the yield point. G0 c
e
as e c 1 exp( ) , where p and q are characteristic
parameters. Nevertheless, if the shear rate changes stepwise,
the value of also changes stepwise, and so do the values
of e and elastic stress e. To overcome this problem, here
we assume that the evolution of e
c 1 exp(
) , where
is the total shear strain
and it has
( ) d . Then e
) , where c is the shear strain at
e 0 c
of the elastic yield stress. Finally, the elastic stress e has the
following form:
e
1 exp(
)
(1)
For the viscous stress v, it is assumed as in the model
k and a completely
unstructured consistency k. The dependence of v on shear rate
n1. The relative importance
of the contributions of e to v is determined by the structural
parameter (Mujumdar et al, 2002). Finally, the equation of
state takes the form as follows:
2.2 Kinetic equation
As for the kinetic equation for the structural parameter
, in analogy with chemical reaction kinetics, it is usually
assumed that its time evolution is controlled by the combined
result of structure buildup and breakdown rates. So far
many kinetic equations have been put forward, and some
of them introduce a common prefactor for the buildup and
breakdown terms (Mewis and Wagner, 2009; Jia and Zhang,
2012; Teng and Zhang, 2012). The prefactor does not change
the equilibrium value of at each shear rate. It only changes
the descending rate of coming to its equilibrium value. By
this way, it substantially improves the predictive capability
of the model. In this work, the prefactor 1 / (1 2 ) with
Zhang, 2012), is introduced. Then the structural parameter
is assumed to obey the following kinetic equation:
d
d
term of the structural kinetic equation is a function of the
where is the total shear strain; n2 is a positive dimensionless
constant; a is a kinetic constant for structure buildup. In this
equation, the first term on the righthand side represents
the structure buildup, while the second term stands for the
structure breakdown.
has the following
characteristics: (i) the breakdown rate is low at early times
when the microstructure is nearly undeformed; (ii) as time
elapses and the shear strain increases, the breakdown rate
increases; (iii) near the yield point the breakdown rate comes
term ( , ) is a function of shear rate only (Houska, 1981;
Toorman, 1997; Coussot et al, 2002; Mujumdar et al, 2002)
If the function ( , ) is taken to be dependent on shear rate
(e.g. ( , ) ),
2011), ( , ) = 0, and increases monotonically
with the increase of . For a completely structured material
initially at rest, if it is suddenly subjected to a constant shear
rate (artificially assumed here), the structure breakdown
it begins to decrease. However, if the function ( , ) is
dependent on shear rate only, ( , ) does not change with
time. As a consequence, the breakdown term ( , ) is a
=1
and ( , ) is constant), and decreases as decreases. This
is not consistent with the actual evolution of . Therefore, de
nonphysical to assume that the destruction of microstructure
is a function of shear rate. One thing to note is that the
nonphysical responses mentioned above lies in the viscoelastic
part before the yield point. Therefore, it does not show up in
considering the viscoelastic part before the yield point.
To overcome this problem, the two models proposed by
regime (away from the yield point) when the network
by an isolated floc is proportional to the shear rate raised
to a power (Barnes, 1997). Therefore, in this regime it
is reasonable to assume that the breakdown term of the
structural kinetic equation is dependent on shear rate. This
is also the reason why most of the kinetic equations in
the viscoplastic model assume that the breakdown term is
dependent on shear rate.
To take consideration of the initial elastic dominated
regime, here we assume that ( , ) is a function of a
combination of shear stress and shear rate, e.g. the rate of
energy dissipation , which is defined as in the
(3)
stress will increase with time until it reaches its peak value,
the energy dissipation rate will increase with time, and so
do the values of ( , ) and ( , ) , and the breakdown
term ( , )
shear stress point. This is in line with the actual evolution of .
Based on the above discussion, we propose the following
kinetic equation to depict the evolution of the structural
parameter with time:
d
d
1
1
2
where b is a kinetic constant for shearinduced breakdown
and m is a dimensionless constant.
In summary, the proposed model is composed of Eq. (2)
and Eq. (4), and contains ten adjustable parameters, i.e. y, k,
k, n1, n2, a, b, m, p, q. The values of the model parameters
curve, with the kinetic equation for the structural parameter
solved numerically by the fourth order RungeKutta
3 Materials and methods
The crude oils used in this study were Daqing and
oils in China. In consideration of the thermal and
shearthe oil specimens were pretreated to a temperature of 80 oC
of the oil for better repeatability and comparability of the
statically for 48 h before testing.
cylinder sensor system (Z41Ti) of a stresscontrolled
rheometer (HAAKE MARS III) and the temperature of
samples was controlled by a programmable water bath
(AC 200) with temperature control accuracy of 0.01 °C. A
pretreated oil specimen was heated to 50 °C and held at that
temperature for 20 min. It was then loaded into the measuring
cylinder preheated at 50 °C and held isothermally for 10 min,
and then statically cooled to the test temperature at a cooling
rate of 0.5 °C/min. At test temperatures, the oil specimen
was held isothermally for 45 min before testing to let the
measurements were performed isothermally. In this work,
temperatures were 34, 35, and 36 °C. For each test a fresh
specimen was used.
4 Model validation
In this section, the proposed model was validated
respectively by stepwise shear rate tests and hysteresis loop
tests. Furthermore, to check the validity and applicability of
the model, we used the model parameters obtained from the
hysteresis loop test.
4.1 Stepwise shear rate test
70
60
50
40
30
20
10
0
1.0
0.8
0.6
0.4
0.2
0
test, since the coupled effects of time and shear rate/stress
In this work, the stepwise shear rate test is adopted and
the chosen shear rates are 1, 2, 4, 8, 16, 32, and 64 s . The
(4)
at 34 °C are shown in Figs. 1 and 2, respectively. As we are
mostly interested in the viscoelastic regime at small shear
strains and the postyield region, the shear stress is plotted
against the shear strain in the logarithmic scale and against
the time in the linear scale. It can be seen in Figs. 1 and 2 that
within 2.0%.
20
15
a
P
,
s
s
trse 10
r
a
e
h
S
5
0
20
15
a
P
,
s
s
trse 10
r
a
e
h
S
5
0
100
Shear strain
From Figs. 1(b) and 2(b), it can be observed that at
the beginning of deformation, the mechanical response is
mainly viscoelastic, and the shear stress increases with an
increase in shear strain. The total shear stress predicted by
the model is mainly made up of the elastic stress e, with
Measured
Fitted
1000
2000
Time, s
the viscous stress v very small. As the shear strain further
increases, the structural parameter begins to decrease and
the mechanical response transfers from the viscoelastic
accompanying the yielding process. For the model, the
predicted total shear stress changes from being dominated by
e to being dominated by v. While at higher shear rates, is
very small, and the predicted shear stress is mainly composed
of v, with e very small. In summary, the model predicts the
of structural parameter with time is also shown in Fig. 3.
Figs. 1 and 3 indicate that at low shear strains, the breakdown
rate of is very small and it increases as the shear strain
increases. Near the yield point, the breakdown rate of
comes to its peak value, as can be observed in Fig. 1(b).
After that, the structure buildup rate starts to increase while
the structure breakdown rate decreases. As time goes on, the
structure buildup and breakdown gradually reach dynamic
equilibrium. Hence, the total shear stress gradually decays to
its equilibrium value. For the following shear rates, when the
shear rate is suddenly changed to another higher value, the
shear stress also suddenly jumps to a higher stress owing to
the jump of shear rate. The structure breakdown rate is also
stepwise elevated, and then gradually decays as gradually
comes to its equilibrium value at each shear rate. The total
shear stress at each shear rate also gradually decays to its
equilibrium, as shown in Fig. 1.
70
100
102
t
d
/
d
&
104
106
o
d /dt
104
102
100
Shear strain
102
104
4.2 Hysteresis loop test
The hysteresis loop measurement is a common method
(Barnes, 1997; Labanda and Llorens, 2006). The hysteresis
loop test in this study consists of two cycles of linear increase
of shear rate from 0 to 50 s1 in 200 s and then linear decrease
to 0 s1
total data range. The AADs are within 5.0% for Daqing and
Like the stepwise shear rate test, the initial shear stress
shows up a stress overshoot and undershoot, indicating the
the elastic stress e decreases sharply to a very small value,
and the total shear stress predicted by the model gradually
changes from being dominated by the elastic stress e to being
dominated by the viscous stress v
4.1.
To further estimate the proposed model, here we used the
model parameters obtained from stepwise shear rate test to
predict the transient flow curve of the hysteresis loop test.
at 33
behavior of the hysteresis loop test can be described by the
model on the whole. The AADs are within 18% for Daqing
Measured
Fitted
(a) Shear stress vs. time
200
600
800
104
102
100
102
104
Shear strain
(b) Shear stress vs. shear strain
and the measured data is relatively high. This is caused by the
fast decay of the structural parameter .
5 Conclusions
parameters is proposed, based on the Mujumbar model and
the Zhu model. The proposed model is validated by the data
from stepwise shear rate tests and the hysteresis loop tests of
the average absolute deviations of the two measurements are
within 2.0% and 5.0%, respectively. The model depicts the
whole rheological response before and after the yield point
well without a discontinuity, from the initial elastic dominated
region accompanied with the yielding process, and finally
to the viscous dominated region. Moreover, the model is
further validated by predicting the transient flow curve of
the hysteresis loop test with the model parameters obtained
from the stepwise shear rate test, and the results show that
the agreement between the predictions and measurements is
satisfactory.
Acknowledgements
National Natural Science Foundation of China (Grant No.
Petroleum (Beijing) (Grant No. LLYJ201155).
589
Mechanics . 1997 . 70 ( 12 ): 1  33 Industrial & Engineering Chemistry Research. 1998 . 37 ( 4 ): 1551  1559