Modeling the viscoelasto-plastic behavior of waxy crude

Petroleum Science, Aug 2013

Waxy crude oil exhibits complex rheological behavior below the pour point temperature, such as viscoelasticity, yield stress, and thixotropy, owing to the formation of a three-dimensional spongelike interlock network structure. This viscoelasto-thixotropic behavior is an important rheological behavior of waxy crude oils, determining the flow recovery and safe restart of crude oil pipelines. Up to now, the thixotropic models for waxy crude have been all viscoplastic models, without considering the viscoelastic part before the yield point. In this work, based on analyzing the variation of the elastic stress and viscous stress in the Mujumbar model, a new viscoelasto-plastic model is proposed, whose shear stress is separated into an elastic component and a viscous component. The elastic stress is the product of the shear modulus and elastic strain; the shear modulus is proportional to the structural parameter. For the elastic strain, we followed the line of Zhu and his coauthors and assumed that it may be expressed by an algebraic equation. The model is validated by stepwise shear rate tests and hysteresis loop tests on Daqing and Zhongyuan waxy crude. The results show that the model’s fitting and predictive capability is satisfactory.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs12182-013-0287-0.pdf

Modeling the viscoelasto-plastic behavior of waxy crude

Received January Modeling the viscoelasto-plastic 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 three-dimensional 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 time-dependent elastic, viscous, and yielding phenomena. The total shear stress is separated into elastic and viscous components which originate from intra-floc 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 time-dependent 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 viscoelasto-plastic 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 broken-down 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 structure-dependent 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 right-hand 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 un-deformed; (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 non-physical 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 shear-induced 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 Runge-Kutta 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 stress-controlled 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 post-yield 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 10-2 t d / d & 10-4 10-6 o d /dt 10-4 10-2 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 s-1 in 200 s and then linear decrease to 0 s-1 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 10-4 10-2 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. LLYJ-2011-55). 589 Mechanics . 1997 . 70 ( 1-2 ): 1 - 33 Industrial & Engineering Chemistry Research. 1998 . 37 ( 4 ): 1551 - 1559


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs12182-013-0287-0.pdf

Houxing Teng, Jinjun Zhang. Modeling the viscoelasto-plastic behavior of waxy crude, Petroleum Science, 2013, 395-401, DOI: 10.1007/s12182-013-0287-0