energies Article A Nonlinear Viscoelastic Model for the Yielding of Gelled Waxy Crude Oil Mengran Sun 1,2,3, * , David Jou 4 and Zhihui Wang 5 1 National Engineering Laboratory for Pipeline Safety, Beijing 102249, China 2 MOE Key Laboratory of Petroleum Engineering, Beijing 102249, China 3 Beijing Key Laboratory of Urban Oil and Gas Distribution Technology, China University of Petroleum-Beijing, Beijing 102249, China 4 Department of Physics, Autonomous University of Barcelona, 08193 Barcelona, Spain; david.jou@uab.cat 5 Science and Technology Information Department, China Huanqiu Contracting and Engineering Co., Ltd., Beijing 100081, China; wangzhihui@hqcec.com * Correspondence: sunmengran@cup.edu.cn Abstract: We explore some rheological aspects of the yielding of gelled waxy crude oil on the basis of a fractal model for the structural description of the waxy gel and Marrucci’s model for the time evolution of the stress with mixed elastic and viscous effects. With some parameters of the model directly obtained from classic rheometry, and others by fitting the parameters to the experimental data of one shear-rate condition, the flow curves for another shear-rate condition are predicted. Both theoretical curves—the fitting and the predicted ones—share the basic features of the experimental ones. Comparison with results of Maxwell model shows that Marrucci’s model used here leads to much better results, as it incorporates nonlinear viscoelasticity of waxy crude gels in the stress evolution equation. The strain dependence of the elastic modulus also plays a relevant role on the prediction of the model, suggesting a double-network contribution for very small strain values. Due to the inertia of rheometric device, the actual shear rate is often found to depart from the setting  one, and modification of shear rate history can be necessary in model validation.  Citation: Sun, M.; Jou, D.; Wang, Z. A Keywords: nonlinear viscoelasticity; rheological model; waxy crude oil; yielding; fractal gels Nonlinear Viscoelastic Model for the Yielding of Gelled Waxy Crude Oil. Energies 2021, 14, 536. https://doi.org/ 10.3390/en14030536 1. Introduction The rheological behaviors of gelled waxy crude oils during yielding process is a Academic Editor: Mehrdad Massoudi very relevant problem for theoretical rheology, as it exhibits the consequences of drastic Received: 7 December 2020 structural modifications at a microscopic level, leading the system from a solid-like behavior Accepted: 6 January 2021 to a thixotropic fluid behavior. Furthermore, it is also very useful because it may help Published: 21 January 2021 engineers make smart decisions on the operation strategies for pipelines. For instance, before it enters the pipeline, waxy crude oil is often heated to a relatively high temperature Publisher’s Note: MDPI stays neu- to keep it from gelation. While during pipeline shutdown, heat supplement no longer tral with regard to jurisdictional clai- ms in published maps and institutio- enters the pipe and drastic temperature drop could happen. As a result, waxy components nal affiliations. once dissolved in waxy crude may precipitate from the oil when the temperature falls below its wax appearance temperature (WAT). If the shutdown lasts for a long time, irreversible gelation could happen in the pipe since a more stable structure has developed due to high fraction of wax precipitates. Quantitative description of the viscoelastic behavior Copyright: © 2021 by the authors. Li- of a waxy crude gel is helpful in estimating the minimum pressure and time needed to censee MDPI, Basel, Switzerland. restart the flow after a period of shutdown [1–4]. Such practical relevance and consequent This article is an open access article economic implications of yielding behavior add to its theoretical interest and experimental distributed under the terms and con- challenges, making it a topic of great significance. In this work, we aim to model this ditions of the Creative Commons At- nonlinear viscoelastic behavior of waxy crude oil as the yielding of a fractal gel. tribution (CC BY) license (https:// The yielding behavior of waxy crude oil has many complexities [5–11]. The flow curves creativecommons.org/licenses/by/ in Figure 1 (whose experimental details will be fully discussed in Section 3) are typical ones 4.0/). Energies 2021, 14, 536. https://doi.org/10.3390/en14030536 https://www.mdpi.com/journal/energies Energies 2021, 14, 536 2 of 25 of yielding. To quantitatively describe them, two forms of thixotropy models have been developed, both of which consist of an evolution equation for stress and for the structure, characterized by an abstract structure parameter λ, varying from 1 (pure elastic behavior) to 0 (pure viscous behavior). Depending on the rheology of the stress equation, the two categories are often named as the viscoplastic thixotropic models and the viscoelastic thixotropic models. (a) (b) Figure 1. Experimental results of the evolution of shear stress during yielding of an initially gelled waxy crude oil: (a) as a function of shear strain; (b) as a function of time t. Data provided by Zhang [12]. The viscoplastic thixotropic models widely adopted in the oil industry are developed on the basis of a pseudo-plastic model, such as the Bingham model or the Herschel-Buckley model [13–16]. These models consider the waxy gel as a rigid body before yielding and as a viscous fluid with decreasing viscosity with time immediately after a critical stress is loaded. The viscoelastic models, or the elasto-viscoplastic models, incorporate the description of viscoelasticity during pre-yielding period by using stress equation [17–20]. In most of these models, the yield stress term of a steady viscoplastic model is simply replaced with elastic stress. Since the structure-dependent yield stress well reflects the qualitative features of thixotropy during yielding, the predictive ability of these two kinds of models is usually found to be quite satisfactory [21]. The shortcoming of the viscoplastic thixotropic models is obvious. They fail to give a rheological description in the pre-yielding region. Moreover, the presence of the yield stress even after yielding in those models is theoretically problematic. Further problems of the existing viscoelastic-thixotropic models for crude oils are: (1) Their stress equations are not developed on mechanical basis, and the steady-state shear stress-shear rate relationship could not be recovered by setting the time derivatives in the evolution equation to zero; (2) Too many parameters with no specific physical meaning have to be determined through complex numerical fitting. Alternatively, the models based on a linear viscoelastic stress equation developed from classic mechanical models, such as Maxwell model or Jeffrey model, are able to describe the pre-yielding rheological behavior, and thus enhance their theoretical inter- est [22–27], among which the model proposed by de Souza Mendes and R. L. Thompson [27] has been modified and applied to describe the rheological behavior of waxy crude oil by Van Der Geest et al. [28]. However, quantitative description of the gel microstructure evolu- tion in the early stage of yielding does not contribute to any of these models, since most them are not developed on the basis of actual material, so that experimental aspects could not be taken advantage of in understanding the yielding mechanics. As a result, all of them use an ambiguous structure parameter λ to describe the structure, which seems very difficult to validate. Moreover, using a linear viscoelastic stress model (such as Maxwell Energies 2021, 14, 536 3 of 25 or Kelvin-Voigt model), the effects of the thixotropy caused by drastic structural changes during yielding process on the mechanical properties, i.e., the dependence of the elastic modulus G and the viscosity η on time t, is not accounted for. However, since thixotropy is caused by structural evolution of the gel, it should have been taken into consideration as a main character when modeling the yielding of structured gels such as waxy crude oil itself. Mendes et al. [29] have chosen a different way of modelling the rheological behavior of waxy crude oil by involving both flow and temperature history dependence, based on new findings about the structure destruction characteristics from the experimental results of waxy crude oil. Since waxy crude oils are complicated mixtures with special characteristics, their rheological modelling needs further studies. In this paper, for the evolution equation of shear stress, we employ a nonlinear vis- coelastic model proposed by Marrucci and used to model the viscoelasticity of structured polymers [30–34]. We propose an explicit identification for the structure parameter λ, which is often used in an abstract way in other thixotropic models [16–19,22–27], and in- corporate the structure change by relating λ to the several material coefficients in the stress equation. In order to do so, we consider the gel as a fractal structure with a characteristic fractal dimension D. Based on a model for the yielding of weakly flocculated suspen- sions [35,36], we are able to estimate the viscosity and the discussion of the yielding point from a microscopic point of view, in terms of a characteristic radius RF . Another particular aspect we consider here is the time variation of shear rate γ̇ (and therefore of shear strain γ). Instead of assuming γ̇ as an imposed constant along the whole process, we find the inertial transient behavior from 0 s−1 to the final imposed steady-state value of γ̇ with an overshoot in between. We use expressions closer to the experimental observations and analyze their influence on the results. This noteworthy influence occurs in the yielding experiments when the imposed shear rate is driven by a controlled-stress rheometer, which is not rare in References [7,8,20,29,37], considering the fact that controlled-stress rheometers are more common than controlled-strain ones. The article is organized as follows—in Section 2, we discuss the evolution of the structure during yielding from the microscopic perspective, based on which the theoretical justification of the model and its viscoelastic description for both elastic modulus G and viscosity η is proposed. In Section 3, experimental protocols and material properties used to obtain the rheological parameters of the model so as to validate the model are presented. In Sections 4 and 5, we summarize detailed computing method of the model, and then present the results under different hypotheses and comment on the merit and weakness of each one. In Section 6, the concluding remarks are given and the direction of future investigation is discussed. 2. Theory Our aim is to express the shear stress σ as a function of strain γ or of time t in the yielding region. First, we use a fractal model proposed by Snabre and Mills [35,36] to interpret the structural evolution of the waxy crude gel corresponding to its mechanical behavior in yielding, and to analyze the critical strain of yielding, which will be used in the following modelization. For the macroscopic description of the evolution of σ containing both elastic and viscous contributions, we use a thixotropy modelization similar to that in earlier works [23,25,32,38]. 2.1. Microscopic Description In this section, we use a blob model to characterize the flocculated clusters in the waxy crude gel. Under this assumption the critical transition conditions between two successive stages and microstructural changes during each stage will be analysed, among which the most important transition condition—the critical strain of fracture yielding—will be estimated theoretically. Energies 2021, 14, 536 4 of 25 2.1.1. Fractal Model of the Gel Pioneer studies in References [8,39,40] show that rheological properties, for instance, the viscosity and elastic modulus of waxy crude oil, may be well described by fractal models. Also, waxy crude gels and fractal systems have been found to have many morphological features as well as gelation mechanisms in common [9]. Therefore, from experimental point of view, waxy crude oil shares many basic features with fractal models. Here, we consider waxy crude gel as a fractal cluster network constituted by gel clusters of characteristic size RF accumulating primary spherical particles of smaller size a (precipitated wax crystals, in our case), and characterized by fractal dimension D, meaning that the number N of primary particles in each cluster depends on a, and RF as RF D N∼( ) . (1) a Furthermore, following Snabre and Mill’s work [35,36], we assume that the clusters in waxy crude oil may be conceptually considered as spherical blobs of radius RF0 , filling the volume of the system with a volume fraction of φA and a maximum packing volume fraction of φ0∗ (64% for random packing of rigid spheres). The sketch of the configuration of the blobs in different stage of yielding in Figure 2 is to help understand the yielding of the gel from the microscopic prospective. As long as it is in the gel state, the blob radius is assumed to stay RF0 , and the volume fraction occupied by the blobs φA will always be φ0∗ . After the blobs are broken and the system enters the sol state, φA will decrease with RF following the time evolution equation which may be obtained from the kinetic equation of the structure parameter λ (see Section 2.2.2). Because of the fractal assumption Equation (1), both φA and φ0∗ may be related to the volume fraction of the primary particles φ0 as [8] R F D −3 ∗ R φA ∼ φ( ) , φ0 ∼ φ( F0 ) D−3 . (2) a a Thus, the relationship between φA and RF , both before and after the breaking of the blobs, could be described by R F 3− D φA = φ0∗ ( ) . (3) a 2.1.2. Microstructure Evolution and Characteristic Transition Points As from the macroscopic point of view, the mechanisms during different stages of yielding—elastic deformation region, creeping stage, and fracture stage—are significantly different, we put our focus on interpreting what happens to the microstructure in each stage, and discuss the rheological behavior separately as well. From the flow curves in Figure 1a, four distinct regions can be identified: 1. The initial overlapping part of the two curves (corresponding to different nominal values of the shear rate γ̇ )—from the beginning up to point A—can be interpreted as an elastic deformation of the gel, with an infinite or extremely high viscosity; 2. The region between point A and peak-stress point B can be interpreted as to a viscoelas- tic creep, where elastic modulus decreases and viscosity becomes finite, but elastic behavior still dominates; 3. In the region beyond the peak point B up to the inflexion point C, viscous flow becomes increasingly dominant and the microstructure of the initial gel breaks in a steep way, leading to a reduction in both elastic modulus and viscosity; 4. The region between point C and the steady-state point D corresponds to thixotropic behavior of the the remaining structure, mostly a viscous behavior with the viscosity decrease slowing down until a steady-state is achieved. Upon the blob assumption of the microstructure of the gel proposed in Section 2.1.1, we may interpret the critical points on the flow curves in Figure 1a as follows: (1) point A Energies 2021, 14, 536 5 of 25 as the elastic-limit point where bonds between gel particles inside the blobs begin to break (when the spherical blobs are increasingly turned into ellipsoidal blobs); (2) point B as the structure fractural yielding point (corresponding to the strain value γy ) where the blobs become able to move relatively to each other because of the removing of the geometric constrictions to their motion so that the gel network collapses as a whole and the blobs begin to break; (3) point C (corresponding to the strain value γc ) as the inflexion point, where breaking process begins to slow down. From the beginning up to point A, the gel deforms elastically and one may imagine this as a relative deformation of the blobs with the density of bonds between them staying constant (Figure 2a). From point A to Point B, further deformation breaks internal bonds and stretches out the blobs from spheres into ellipsoids in the direction of shear. The break- ing of the bonds will lead to a reduction of the elastic effects, while the orientation of the blobs will make more space for them to move, thus generating viscous effect. But as long as γ < γy , the bond density still stays high enough to hold all those primary particles inside the blobs where they originally belong to (Figure 2b). And since the volume fraction of the blobs does not change, it will remain in the gel phase. For γy < γ < γc , that is, from point B to C, the bonds break with an increasingly reduced proportion and there are small particles leaving their original blob. Thus, the blobs become small enough to freely move with respect to each other, and the system enters the sol phase (Figure 2c). For γ > γc that is, from point C to D, the bonds between particles in individual blobs become weaker due to larger strain and more intense shear, and thus the radius of the blobs will decrease with time (Figure 2d). (a) (b) (c) (d) Figure 2. Blob model illustrating the microstructure evolution during yielding: (a) intact gelation structure; (b) microstructure during elastic dominant viscoelastic creeping; (c) microstructure with tangency geometry at critical point of structural yielding; (d) breaking of colloidal aggregates after the structural fracture of gel. 2.1.3. Estimation of the Yield Strain While the thixotropy models of waxy crude oil have always given a satisfactory description of the fourth stage (from C to D), our focus will, alternatively, be the first three Energies 2021, 14, 536 6 of 25 stages, that is, from the beginning to C, which are closely related to the transition from gel to sol state. To estimate the value of γy , let us consider four of the spherical blobs in contact in the compact gel, each one of radius RF0 , two in the upper line and the other two in the lower line (see Figure 3a). The separation l between√ the upper and lower horizontal planes in tangent contact with the spheres is l ∼ (2 + 3) RF0 . In the configuration of the Figure 3a, the upper sphere cannot move with respect to the lower ones, because of geometric restrictions. In Figure 3b, two deformed ellipsoidal blobs may stand vertically one over the other, and thus relative motion is enabled. In such critical configuration, the critical radius RFc of the inscribed circle of any√blob ellipse, whose semi-minor axis is equal to RFc , may be given by RFc = L/4 = ((2 + 3) RF )/4 = 0.93RF0 . A direct estimation of critical strain is given by the tangent of the angle θ 0 formed by the lines AC and BC in Figure 3a. In order that the upper sphere may slide over the lower one (with some vertical contraction of both of them), the contact point A of both√spheres must arrive at position B. The corresponding tangent is given by tan θ 0 = (2 + 3)−1 = 0.27, which provides an alternative estimation for the critical strain. In fact, the actual yield point is expected be somewhat smaller than 0.27 if one takes into account that the spheres will be compressed in the vertical direction as they are sheared one with respect to the other. Alternatively, we may estimate the yield strain in a more accurate way. The deforma- tion is a combined effect of shear and contraction. We assume a horizontal displacement ∆x along the upper surface due to simple shear, plus a vertical contraction. In order to start relative motion, the blob must be tilted with an angle θ, whose cosine may be ex- pressed as cos θ = RFc /RFt , with RFt half the distance between the two points (B’ and C’ in Figure 3b) on the elliptical blob in contact with neighboring blobs, and RFc equal to its critical vertical projection. The vertical projection 2RF = 2 × 0.93RF0 (as obtained in the former paragraph), and B0 C 0 may be taken 2RF0 as an approximation, in such a way that cos θ = 0.93, and γ = ∆x/2RF ≈ tan θ = 0.40. Actually, the distance between B’ and C’ is slightly lower than 2RF0 , so the actual yield strain must be lower than 0.40. We may estimate this strain more precisely with the aid of analytic geometry, and the results is shown to be 0.25 (see Appendix A). In our analysis we will take γc = 0.25 (see Table 3). The experimental results of yield strain 0.2 at 31 ◦ C to 1.0 at 35 ◦ C [41] are indeed in this range. (a) (b) Figure 3. Sketch for an estimation of the yielding shear strain: (a) before elongation; (b) at critical elongation. 2.2. The Model In this section, the theoretical concept of the development of the model will be given in detail. Based on a viscoelastic model for the description of mechanics, a kinetic equation for the structure and two material functions for the rheological parameters of the gel at a certain structured level are also proposed.In the viscoelastic model, two material parameters appear—the elastic modulus G and viscosity η, which depend on the structure of the material, and on γ and γ̇. To close the description, an evolution equation of the Energies 2021, 14, 536 7 of 25 structure and its relation to G and η are needed. Thus, two evolution equations, for the stress σ and for the structure parameter λ are needed, complemented with functions for G (λ, γ, γ̇), η (λ, γ, γ̇) appearing in the stress equation, and for the building-up and the destruction terms of structure K + and K − , appearing in the evolution equation for λ, respectively. For the sake of higher efficiency, we will slightly modify the general scheme of modeling the kinetics for the structure by expressing it in terms of γ. 2.2.1. Mechanical Description: The Nonlinear Viscoelastic Model In previous models, the evolution of shear stress σ is often described by linear vis- coelastic models like Maxwell-like or Jeffrey-like models, which do not account for the time-dependent behavior of G. In our model, instead, the shear stress σ will be described by the nonlinear viscoelastic model d(σ/G ) σ+η = η γ̇. (4) dt Equation (4) could be seen as the nonlinear version of Maxwell model. It has the same mechanical analogical model as Maxwell model, only with the time dependence, or the strain dependence of the material parameters considered. And it reduces to the Maxwell model upon constant η and G. This kind of model was first used by Marrucci et al. [30–34] for the description of structured suspensions viscoelasticity, but so far has not been applied to waxy crude oil, whose elastic modulus G dependence on time t or on shear strain γ should not be ignored as well due to the dramatic structural changes during its yield- ing. The shear modulus G (γ) and the viscosity η (γ̇) will be given in Section 2.2.3 as Equations (14) and (15) respectively. We will show that Equation (4) leads to better results than using the linear Maxwell equation in the discussion Section 5.2. 2.2.2. Kinetic Equation and Physical Interpretation of the Structure Parameter λ The second aspect to be considered is the structure evolution of the system, which in its turn influences G and η. Since the microstructure is very complicated, it is usually simplified as a scalar parameter λ which goes from λ = 1 for the completely structured gel to λ = 0 for the completely unstructured fluid [16–19,22–26]. Of course, this is a very radical simplification, but it may cause ambiguity in modelling and validation. So in this section, we will first make our definition for the structure parameter λ on the basis of the features of the microstructure of the gel. a. Physical interpretation of λ In searching for an explicit physical interpretation for the structure parameter λ, let us recall the microstructure of the gel, which has been discussed in Section 2.1.1. Since RF is the main feature that describes the gel structure, we define our structure parameter λ as RF λ= . (5) RF0 In the system of a fractal gel, while the microstructure is often characterized by size RF , the material parameters such as the elastic modulus and viscosity are mostly affected by the the concentration of gel clusters. Therefore, an expression for λ in terms of the volume fraction of the mentioned blobs φA , rather than in terms of their radii RF is needed. According to the relation between φA and RF described by Equation (3), the structure parameter λ may also be interpreted as 1 λ = (φA /φ0∗ ) 3−D . (6) Note that for the fully gel phase, RF = RF0 and λ = 1; and to the completely unstruc- tured fluid RF = Rss , and λ = λss , as usually required to recover the steady-state viscosity relation. Note, finally, that a yet more detailed analysis could assume a lamellar structure Energies 2021, 14, 536 8 of 25 of the gel instead of an isotropic fractal structure, as observed in microscopic analyses [42], but we will leave these details for future work. Furthermore, the fractal dimension D itself could change along the yielding process, but here we kept it constant. b. Kinetic equation for λ The evolution of the structure parameter λ consists of a build-up term K + describing structure formation and a destruction term K − describing structural degradation. Both terms depend on the shear conditions, and thus may be expressed in terms of stress [25–27], of the shear strain [37], or of the dissipation rate [18,23]. In the literature, this equation is often expressed as dλ 1 1 = (1 − λ)m1 − f (σ, γ̇)λm2 , (7) dt τ1 τ2 with τ1 , τ2 the times characterizing the growth and decay of λ, and m1 , m2 characteristic exponents. f (σ, γ̇) describes the structural evolution dependence on the instantaneous shear condition, which could be assumed, for instance, to be a univariate dependence on γ as References [17,19,37] did in their work. Departing from the assumption that λ is decreasing from the very beginning of loading (t = 0) made in other models [25–27], we assume that λ = 1 as long as the gel structure maintains intact, and only begins to decrease after the critical strain γy for yielding (which has been discussed in Section 2.1.2). Thus, f (σ, γ̇) in Equation (7) could be γ̇ n f (γ̇) = H (γ − γy )( ) , (8) γ̇0 in which H ( x ) is the Heaviside step function with a zero value for x < 0 and unit value for x >= 0, γ̇0 is the shear rate above which most of the blobs have become primary particles, and n is an exponent. If we take m1 = m2 = m for simplicity, we obtain the steady state value of λ 1 λss (σ, γ̇) = n 1 . (9) 1 + [( ττ12 ) γ̇γ̇ss0 ] m In the recent work of Mendes et al. [29], the microstructure evolution was found to depend more directly on shear strain, rather than on shear rate. Moreover, as is discussed in Section 2.1.2 , the yielding point may also be interpreted in terms of shear strain, therefore, instead of expressing the changing rate of structure in terms of t like Equation (7), here we choose to express the evolution of λ directly as a function of shear strain γ. They also proved that the extent of structure destruction depends on the most severe shear condition that the oil gel has experienced during its yielding process, which in our case would be the steady-state shear rate. Thus, we consider a same steady-state structure λss as given by Equation (9), and propose the evolution of structure parameter following relaxation equation of λ as dλ λ − λss (γ˙)ss =− , (10) dγ γ2 where γ2 is a value of the shear strain characterizing the change of λ from the initial value to the steady-state value of λ. Of course, the time-dependence rate of λ may be obtained as in Equation (10), by applying the chain rule dλ/dt = (dλ/dγ)(dγ/dt). This is especially easy if the dependence of shear rate on time γ(˙t) is explicitly given, for instance, if a constant, stepwise or oscillatory shear rate is applied on the system. However, in the particular situation under examination, where the shear rate depends on time in a more complicated way (see Section 4.1), it would take some serious effort to obtain the solution of Equation (7) as a function t, so we focus our attention on Equation (10) instead. Below, once γ̇(t) is given, we may obtain γ(t). This is a further original aspect of the present work, Energies 2021, 14, 536 9 of 25 which allows to consider the σ − γ curve in a more direct way. Anyway, after considering γ̇(t), in Section 6, we will also be able to express the time evolution of the shear stress σ (t). With the physical interpretation of λ in Equation (5), a shear-strain dependence of the evolution of RF , describing how the blob breaking in the region from point B to C in Figure 1a, may be obtained. Since the evolution equation of λ as a function of γ has been given by Equation (10), thus by solving it, one obtains the transition behavior of RF between RF,0 and RF,ss as γ − γ1 γ − γ1 RF (γ, γ̇) = RF0 exp(− ) + RFss (γ̇)[1 − exp( )], (11) γ2 γ2 in which RF0 is the initial value of the blob radius corresponding to the radius of the maximum packing situation (gel state), and RFss the steady-state radius. According to Equation (9), RFss is also to depend on the imposed γ̇ as RF0 RFss = , (12) 1 + ατ γ̇q in which ατ is phenomenological parameter related to the relaxation time of both blob break and particle reunion, γ0 is the critical strain above which most blobs have been breaking into primary particles, and q is an empirical exponent. This expression reflects the obvious fact that the characteristic blob radius will decrease with increasing γ̇, because primary particles will be stripped off more easily from the clusters. Thus, in the initial breaking region B ∼ C, RF decreases with γ, from its initial value RF to its steady-state value RFss . Practically, all the parameters in Equation (12) could be obtained by fitting the steady- state viscosity at several γ̇, which are related to RFss (γ̇) through relation Equation (15). Technically, these parameter could also be obtained through microscopic experiments. As for our analysis, it turns out γ0 = 100 s−1 , ατ = 0.125 and q = 0.2 fits the results well. 2.2.3. Material Functions and Rheological Properties The other two functions needed to describe the process are G and η appearing in Equation (4). The second one is dependent on RF , because viscosity of suspensions is known to strongly depend on the radius of particles. In contrast, no such strong dependence should be expected for G, because in the gel phase RF = R0F , and the only effect on G is that internal bonds between the blobs are increasingly broken with higher γ in the gel phase when the blobs deform. a. Elastic modulus G (γ) We assume the elastic modulus G is proportional to the number density of full chains of bonds times the elastic constant of each bond, and a certain number of such bonds break at a given shear strain. According to this assumption, the nonlinear elastic modulus G (γ) of the system is given by expression G0 G (γ) = , (13) 1 + (γ/γ0 ) p with G0 the equilibrium value of elastic modulus, γ0 the characteristic shear strain, and p the phenomenological exponent. Equation (13) of G exhibits a steep drop at γ0 after a plateau equal to G0 for γ < γ0 , and followed by a very small value for γ  γ0 , with p controlling the steepness of the reduction of G (γ) around γ0 . This form of G (γ) is analogous to that used in References [16,43], with the difference that they used elastic part of shear strain γe instead of the directly measured strain γ. G0 depends on the conditions of gelation, for instance, whether it was in a quiescent state or under shear, and how fast the cooling rate was [8]. It could be tested with rheometer by applying a small amplitude oscillatory shear to the waxy gel sample after it rests at the test temperature for a long time. Comparing the behavior of G (γ) to that of our situation, and since a sudden reduction in G reasonably corresponds to the collapse of the gel network, γ0 will be taken close to the yield strain γy . Energies 2021, 14, 536 10 of 25 As the imposed shear rate is as low as 0.1 or 1 s−1 , it is only reasonable to consider the initial behavior as purely elastic. Therefore, the elastic modulus G within the strain range γ  γ0 could be calculated by G = σ/γ. The calculated elastic modulus shown in Figure 4 is seen to start falling from a very high value at the very beginning of the deformation. This could arise from a pre-collapse of gel network at a rather small strain γ0,2 , which is hardly noticed because of its short-term existence. Thus, one may imagine, that besides the gel network which we have seen to break around the yield strain γy , with an equilibrium elastic modulus G0,2 which can be obtained from rheometric test, there is another more brittle and rigid network in the gel, with smaller characteristic strain of breaking γ0,1 and higher initial elastic modulus G0,1 . Such double network structure would have an elastic modulus equal to the summation of moduli of both networks, each of which is given in the form of Equation (13) G0,1 G0,2 G (γ) = + . (14) 1 + ( γγ ) p1 1 + ( γγ0,2 ) p2 0,1 In practical terms, the critical shear strain γ0,i would imply the strain for which the elastic modulus of the network i drops with a steepness characterized by pi (related to the fractal dimension of each network) due to the breaking of the corresponding type of bonds. Therefore, for small strains γ < γ0,1 , one would have a plateau of modulus G of approximately G0,1 + G0,2 , and between γ0,1 and γ0,2 , there might be another plateau of approximately G0,2 , whose width and position depend on the separation between γ0,1 and γ0,2 and on steepness of the two drop behaviors, that is, on the values of exponents pi . This behavior may be caused by asphaltene or resin macromolecules which coexist with wax crystals. Its influence will be further explored and justified on experimental grounds in Section 5.1. Figure 4. Shear modulus dependence on shear strain G (γ) in the early stage of yielding, obtained from the experimental data of σ and γ in Figure 1a by appling G (γ) = σ(γ)/γ . b. Viscosity η (γ̇) To estimate the viscous effect beyond point A, we use the viscosity model for suspen- sions with viscoelastic aggregates proposed by Snabre and Mills [35], which relates η to the volume fraction of aggregates ηA and deformation tensor D as 1 − φA η (φA ) = η0 φA , (15) [1 − φ0∗ (1+ βε xx ) ]2 where η0 is the solvent viscosity, φ0∗ the maximum packing volume fraction (64% for random close packing of rigid spheres), ε xx is the xx component of the deformation tensor Energies 2021, 14, 536 11 of 25 D, and β a parameter characterizing the viscous effect generated during creeping. For ε xx , we take a simple assumption that it is proportional to γ until the structural failure of the network, beyond which it becomes a saturated constant. Equation (15) is developed on the basis of a viscosity model for suspensions with rigidly flocculated blobs [35], by taking into account the effect of the elastic deformation of the blobs. Instead of being constant as in the model of the rigid-blob suspension, the effective maximum packing volume fraction φ∗ in Equation (15) is assumed to increase with deformation as φ0∗ (1 + βγ). This correction makes the volume fraction of the blobs φA smaller than φ0∗ so that a high but finite viscosity is expected before fracture yielding. It reflects the physics discussed in Section 2.1.1—since the blobs are elongated and become more orientated, bigger room is made for their movement. Thus, the internal flow through the network becomes less congested, leading to a higher fluidity, but still an inconspicuous viscous effect before the blobs are broken. 3. Materials and Experiments Rheometric tests are taken to obtain the data we need to verify the model and differ- ential scanning calorimetry (DSC) is used to detect the wax fraction precipitated from the oil varying with temperature. The material we use in the experiments is a typical waxy crude oil sampled from the Daqing oil field in the northeast of China, which is famous for its high wax content. In this section, a detailed description of the material as well as the procedure of the rheometry experiments will be given. 3.1. Materials The wax appearance temperature (WAT) of the specimen we use is 44.8 ◦ C. It has a total wax mass fraction of 21.7%, with only 1.1% mass fraction precipitated by the test temperature 34 ◦ C. In Figure 5, the concentration profile of alkanes with different number of carbon atoms obtained by gas chromatography shows more detailed composition of the waxy crude sample. In order to eliminate the memory effect of flow and thermal history, all the samples were pretreated before the measurements. The samples sealed in 100-mL-gas-tight bottles were heated in a 80 ◦ C water bath for 2 h, statically cooled to room temperature, and maintained for 48 h before use. Figure 5. Concentration profile of alkanes with different numbers of carbon atoms in the crude oil sample. 3.2. Differential Scanning Calorimetry The Differential Scanning Calorimetry (DSC) test of the crude oil sample was con- ducted with a modulated DSC instrument (TA 2000/MDSC 2910,New Castle, DE, USA), and according to the regulation SY/T 0545-2012, China’s standard method for differential scanning calorimetry of the oil and gas industry. With the data of heat flow transmitted at each tested temperature (Figure 6a), the mass fraction of the precipitated wax as a function Energies 2021, 14, 536 12 of 25 of temperature was obtained (Figure 6b). Figure 6a shows that the WAT (wax appearance temperature) of the waxy oil, corresponding to the point where heat flow slope turns non- zero, is 44.8 ◦ C. Figure 6b shows the sample has a total mass fraction of wax φtotal = 21.7%, but only 1.1% mass fraction has precipitated by the test temperature T = 34 ◦ C. (a) (b) Figure 6. Differential Scanning Calorimetry (DSC) curves of the waxy crude oil obtained in laboratory: (a) heat flow varying with temperature T. (b) mass fraction of precipitated wax varying with temperature T. 3.3. Small Amplitude Oscillatory Shear (SAOS) A controlled-stress rheometer (HAAKE MARS III, Karlsruhe, Germany) equipped with a sensor system Z38 Ti of geometry of concentric cylinder was used to test the viscosity and shear modulus. The pretreated specimen sealed in blue cap bottles were heated from ambient temperature to 50 ◦ C (higher than WAT) in the water bath and were kept for 30 min. Then it was loaded into the concentric cylinder gap preheated to 50 ◦ C. The oil samples were cooled at a rate of 0.5 ◦ C/min from 80 ◦ C to 30 ◦ C, loaded with a small amplitude oscillatory shear, of a frequency of 1 Hz and a dynamic shear strain of 0.0002, which is small enough to preserve the structure. Meanwhile, the storage and loss moduli as well as the complex viscosity were tested and recorded. The data are presented in Figure 7. According to Figure 7, the gelation point of the waxy crude oil is 36 ◦ C, and the order of dynamic moduli at this point are as small as 1 Pa, meaning that a fully developed gel structure has not been achieved and yielding behavior is not notable. The yielding experiments were carried out 2 ◦ C below, at which the tested storage modulus is 41 Pa, and the loss modulus is 14 Pa, and the equilibrium elastic modulus of the oil gel can be much higher than that as a result of ageing. One could also read from Figure 7 that the complex viscosity is 0.016∼0.028 Pa·s above WAT, which may be taken as a reference range for the solvent viscosity. 3.4. Viscosity-Temperature Relation Viscosities of the waxy crude oil sample were tested at every 5 ◦ C from 37 ◦ C down to 35 ◦ C with the rotational rheometer (Rheolab QS, Anton Paar, Graz, Austria). The pretreated oil sample was loaded in the rheometer which had been heated to 60 ◦ C for 15 min and then was cooled to each test temperature at a cooling rate of 0.5 ◦ C/min and stood still for 30 min. At each test temperature, viscometry was conducted under 3 different constant shear rate, 20 s−1 , 50 s−1 and 100 s−1 . The measured viscosity was recorded every 5 min until the difference between two successive tested viscosities was below 5%, that is, until the steady state was achieved. The steady-state viscosities at all the test temperatures are presented in Table 1. From Table 1, one could see that the viscosities of different shear rates at each temperature above coincide, while the viscosity of different shear rates at 37 ◦ C Energies 2021, 14, 536 13 of 25 have separated (not shown here), indicating the test temperatures between 40 ◦ C and 70 ◦ C are above the abnormal point, above which the flow behavior is completely Newtonian. Table 1. Newtonian viscosity at different temperatures, tested under different shear rates 20 s−1 , 50 s−1 , 100 s−1 . Data provided by Zhang [12]. Temperature ◦ C Viscosity Pa·s 70 0.00785 65 0.00810 60 0.00855 55 0.01005 50 0.01256 45 0.01580 40 0.02497 Figure 7. Storage and loss moduli of the waxy crude oil varying with temperature. 3.5. Rheological Flow Curves The flow curves presented in Figure 1, relating shear stress σ to shear strain γ and to time t, were obtained with the coaxial cylinder sensor system of a controlled-stress rheometer (HAAKE RS 150,Karlsruhe, Germany). The oil sample was statically cooled from 60 ◦ C to the test temperature 34 ◦ C and stood for 2 h until steady gel structure formed. Then the rheometer was set to a constant shear rate (the actual shear rate will increase from zero to its set value with delay due to system inertia, though) 0.1 s−1 and 1 s−1 respectively, and the waxy oil gel is exposed under shear. As a particular stimulus for our work, we consider the flow curves σ − γ and σ − t in the yielding region, corresponding to the shear strain range 0.001∼10. 4. Results This section deals with the computing methods of the shear stress σ as a function of shear strain γ, or of time t, for different shear rate histories γ̇(t) (or equivalently, for dif- ferent input functions for γ̇(γ). The connection between σ (γ) and σ(t) comes from the dependence of γ on t, in such a way that we may give σ(γ), or σ(γ(t)) = σ (t). Also, we will discuss an unusual finding related to the imposed shear rate history we encountered in analysing the experimental data. The evolution equations, material functions proposed in previous chapters are summarized in Section 4.2. Later, the parameters in each functions, both the directly and the indirectly obtained ones, are given respectively in Tables 3 and 4. Energies 2021, 14, 536 14 of 25 4.1. Imposed Shear Rate History The input function on the boundary, that is, γ̇ or γ as a function of time t, plays an influential role on the results of the model prediction. The different curves in Figure 4 are supposed to correspond to different constant values of γ̇, but it will always take some time to get the actual value of imposed shear rate γ̇ to its set value γ̇∞ , especially for the controlled- stress rheometers. And there might even be overshoots around the yielding point. In previous literature, this feature has only been briefly mentioned by Tracha et al. [11], and the influence of the particular history of γ̇(t) still needs detailed investigation. In experiments, the actual γ̇ starts from a small value close to 0, and for a short time stays very low, when the gel has an infinite (or very high) viscosity. When yielding happens, viscosity reduces by several orders of magnitude within rather narrow strain interval, and the stress applied to the system suddenly becomes higher than the viscous resistance, so the shear rate γ̇ will increase sharply around yield strain. Eventually, it will reach the time-independent set value γ̇∞ imposed by the controlled-stress rheometer. The dynamical reasons for the delay and the overshoot might be inertia of the experimental device, and it will be examined in the future. Here, for the shear rate history under examination, the experimental curves of γ̇(γ) will be mathematically approached as a empirical equation γ m1 γ̇(γ) = γ̇00 + (γ̇∞ − γ̇00 ) , (16) 1 − aγm2 + γm1 in which γ̇00 is a small initial, m1 , m2 phenomenological parameters controlling the increas- ing slope of the γ̇(γ) curve, and a relating to the overshoot value. If a ≤ 1, Equation (16) predicts a shear rate history γ̇γ starting from γ̇00 and increasing smoothly without over- shoot to its set value γ̇∞ , whereas if a > 1, it predicts a overshoot in shear rate history. To examine the role of shear rate plays in the model prediction, shear rate history with and without overshoot are both approximated with Equation (16). The values of the parameters are listed in Table 2 and the fitting results are presented in Figure 8. The dash line (color: navy blue) corresponds to shear rate history without overshoot, the solid line (color: light grey) to the shear rate history with overshoot, and the curve with discrete points (color: light blue) the observed one. The time evolution of shear strain γ(t) could be obtained by integrating Equation (16) with respect to time t. Table 2. Values of parameters of the function Equation (16) for shear rate history. m1 m2 a γ̇00 0.1s−1 no 1.65 0.1 1 0.0004 overshoot 1 s−1 no 1.5 0.15 0.7 0.003 overshoot 0.1 s−1 1.4 0.55 1.8 0.0004 overshoot 1 s−1 overshoot 1.4 0.9 1.55 0.003 4.2. Summary of the Model First, let us recall the evolution equation for stress σ in terms of the two material functions G (γ), η (φA ), and the evolution equation for the structure parameter λ = RF /RF0 , which is given by the particular solution Equation (9) of Equation (7), and related to φA Energies 2021, 14, 536 15 of 25 through the structural relation Equation (6). Thus, the model we have proposed in the above sections becomes d(σ/G ) σ+η dt = η γ̇ (17) G (γ) = ΣG0,i 1+(γ/γ 1 p (18) 0,i ) i 1−φ0∗ λ3− D η (λ, γ) = η0 3− D (19) (1− λ1+ βγ )2 dλ λ−λss (γ̇ss ) 1 dγ = γ2 , λss = 1+ατ (γ̇ss /γ̇0 )q (20) m1 γ̇(γ) = γ̇0 + (γ̇∞ − γ̇0 ) 1− aγγm2 +γm1 . (21) The input of the model is the shear rate γ̇(γ) described by Equation (21) with param- eters listed in Table 2 in Section 4.1, whereas the output is the shear stress σ (t) varying with time t, which is obtained in the experiment mentioned in Section 3.5. The time-, strain-dependent rheological variables necessary to obtain the output include the shear modulus G (γ) given by Equation (18) and the viscosity η (γ, λ), which may be obtained from Equation (19) along with Equation (20), in which λss is estimated to be 0.95 for γ̇ = 0.1 s−1 , and 0.97 for γ̇ = 1 s−1 by using Equation (12) and Equation (5). (a) (b) (c) (d) Figure 8. Fitting curves of shear rate history with and without overshoot in comparison with measured ones: (a,b) for γ̇ = 0.1 s−1 ; (c,d) for γ̇ = 1 s−1 . Equation (17) will be solved by the finite-difference method to obtain evolution of stress σt−dt + Gt ∆tγ̇t σt = , (22) 1 + Gt ∆t/ηt − ( γ̇∆t dG Gt )( dγ )t Energies 2021, 14, 536 16 of 25 with σt the stress at time t. 4.3. Summary of the Parameters The parameters of the presented model fall into two categories. Ones are rheological parameters with explicit physical meanings which can be directly measured or theoretically estimated under certain assumptions. Actually, their values are objective and are not supposed to be obtained through fitting process. The other category are the unknown parameters which are indirectly obtained through fitting of experimental curve, which in our case, is one of the curves for 0.1 s−1 in Figure 1. 4.3.1. Directly Obtained Parameters Five out of eleven parameters in the model may be obtained directly, including: (1) The equilibrium elastic modulus of the relatively resistant network G0,2 , the char- acteristic strains γ0,1 , γ0,2 in Equation (18); (2) Solvent viscosity η0 in Equation (19); (3) The characteristic strain γ1 and D in Equation (10). Among them, the characteristic strains γ0,2 and γ1 which are very much relevant to yielding are taken 0.25, according to the estimations of yield strain γy in Section 2.1.3. The equilibrium shear modulus of the resistent network is taken G0,2 = 177 Pa, the equi- librium elastic modulus of Daqing crude oil at 34◦ C as is tested in Teng and Zhang’s previous work [16]. We take the maximum packing volume fraction φ0∗ = 0.64 as Snabre and Mills [35,36] took. As for the fractal dimension of waxy crude gels, Yi and Zhng [44] obtained the box fractal dimension in experiments to be around 1.5, and thus it is taken here. In Section 5.4, we will take D = 2 to see the consequences of different values frac- tal dimension. For the solvent viscosity η0 we employ an Arrhenius-type equation η0 = Aexp( Ea /RT ) to extrapolate η0 at test temperature T = 34 ◦ C from the experimental Newtonian viscosities at several different temperatures above WAT in Table 3. The parameters are A = 5.6 × 10−8 and Ea = −33.4 kJ, obtained by fitting the experimental data. As a result, η0 in the condition of the experiment is theoretically extrapolated to be 0.026 Pa·s (Figure 9). Figure 9. Viscosity variation with temperature above wax appearance temperature (WAT) (exper- imental and ftting curves) and extrapolation of viscosity by Arrhenius equation, with parameters A = 5.6 × 10−8 and Ea = −33.4 kJ. Table 3. Values of directly obtained parameters of the model. G0,2 η0 γ0,2 γ1 D 177 Pa 0.026 Pa·s 0.25 0.25 1.5 Energies 2021, 14, 536 17 of 25 4.3.2. Indirectly Obtained Parameters The other category of indirectly obtained parameters includes 6 unknown parameters: (1) G0,1 , γ0,1 and p1 , p2 in Equation (18); (2) β in Equation (19); (3)γ2 in Equation (20). In this paper, we adjusted their values according to the flow curve σ (γ) in Figure 1a of the shear rate 0.1 s−1 . Since these parameters are restricted to a certain range, rather than confined by certain objective values (at least the values could not be validated through experiments as far as we are concerned), they do not have to be taken the exact same values in all the figures corresponding to different conditions. The parameter values used in each figures in Section 5 are slightly different under the different conditions. The values used for the parameters of the main results in Figure 10 are the listed in Table 4, and others will be given under the each figure being discussed in next section. Compared to a numerical way, these fitting parameters are relatively rough estimations to help us get the main features of the proposed model under a given shear condition (yet it is satisfactory enough as will be shown in the following section). Note that the curves corresponding to γ̇ = 1 s−1 have been predicted by using the same values obtained from the curve for γ̇ = 0.1 s−1 . Their reasonable agreement with the observed curves of agreement with the observed curves of γ̇ = 1 s−1 shows that our model is not only a fit of curves, but also has predictive capability. Table 4. Values of parameters estimated through rough fitting. G0,1 γ0,1 p1 p2 β γ2 1000 Pa 0.03 1.2 1.3 0.01 60 4.4. Fitting and Predicted Results: σ (γ) and σ (t) To obtain σ (t) we use the evolution equation of σ Equation (17), in combination with expressions Equations (18) and (19) for the viscoelastic material properties, and Equation (20) for structure parameter. To examining the actual shear rate in the experiment, and to remove influence of the fluctuations in experimental data, we use the shear rate history with overshoot predicted by the empirical equation Equation (21). In Figure 10 the fitting curves for γ̇ = 0.1 s−1 (red solid line) obtained by the model almost coincide with the experimental ones (red dotted lines). (b) (a) Figure 10. Fitting curves for shear rate γ̇ = 1 s−1 and γ̇ = 0.1 s−1 using the elastic modulus. model under the double- network assumption, and the shear rate history with inertial overshoot Equation (21) (with the fitting values in Table 2): (a) shear stress as a function of shear strain; (b) shear stress as a function of time. The curves for γ̇ = 1 s−1 are predicted by the model with the same parameters used to fit the 0.1 s−1 curves, and are presented in Figure 10 with navy blue solid lines. The overall Energies 2021, 14, 536 18 of 25 shapes of both σ (γ) and σ (t) curves are of good analog with the measured ones (navy blue dotted line), especially in the pre-yielding region. The small deviation after yielding could be raised by the difference between the shear history described by the input function Equation (21) and the actual history. Moreover, it could be seen that the predicted yield stress is 39 Pa, 1 Pa less than the measured one, while the predicted yielding time is 10 s, same as the measured ones, both of which are quite satisfactory. 5. Discussion From Sections 2–4, we present the theoretical justification from both microscopic and macroscopic prospective of the modelization of a nonlinear viscoelastic model used to describe the yielding of gelled waxy crude oil, and the experimental protocols used to obtain the data for model validation. The results of the model prediction shown in Figure 10 is found to be rather satisfactory. Several original aspects we have considered in modelling may have made a contribution: (1) the assumption of double-network gel; (2) the modification of the input shear history which deviates from the setting constant shear rate; (3) the employment of Marrucci’s viscoelastic model for the description of the stress evolution. In this section, we will plot the results of the model prediction in absence of the consideration of the above aspects and by comparing them with the results in Figure 10, the influence of each one of them will be evaluated respectively. 5.1. Double-Network vs. Single-Network Gel In obtaining the theoretical curves in Figure 11, same calculation is used as in Figure 10, except that G0,1 in Equation (18) is taken as 0 so that the gel is viewed as the conventional single network structure, and the value for the remaining network were changed to p2 = 1 to fit the curve as much as possible. By comparing with experimental curves, the theo- retical curves in Figure 10 obtained with the modulus model under the double network assumption (with an additional network whose elastic modulus is G0,1 = 1000 Pa, and char- acteristic shear strain is γ0,1 = 0.03) clearly exhibit a better description than the ones under the single network assumption. (a) (b) Figure 11. Fitting curves for shear rate γ̇ = 1 s−1 and γ̇ = 0.1 s−1 of shear stress as a function of shear strain and of time, taking the whole calculation same as in Figure 10, but using the elastic modulus under the assumption of a single network structure by G0,1 = 0 Pa in Equation (18) and changing value for the remaining network to p2 = 1: (a) shear stress as a function of shear strain; (b) shear stress as a function of time. Generally speaking, the main problem of the single-network elastic modulus model is that it results in a rather low stress before yielding. In Figure 11a, the σ (γ) curves for different imposed shear rates γ̇(γ) begin to separate at a strain of the order of 0.1, one Energies 2021, 14, 536 19 of 25 order of magnitude higher than experimental curves. The fitting curve for γ̇ = 0.1 s−1 reaches a maximum value equal to the experimental one, 30 Pa, while the predicted peak stress on the curve of γ̇ = 1 s−1 reaches a peak stress 3 Pa lower than the experimental one. Also, the shear strains corresponding to the peak stresses for both the fitting curve and the predicted curve are more than one order of magnitude higher than the experimental ones. In Figure 11b, the σ (t) curves show yielding times for both shear rates which are of the same order as experimental ones, with γ̇ = 1 s−1 curve around 10 s and the γ̇ = 0.1 s−1 curve around 100 s. But in the region before the stress peak, in both Figure 11a,b the predicted stresses are relatively small compared to the actual ones, so that the raise of the predicted curves is steeper than that of the experimental curves, while after the maximum value, the decreasing rates plotted with solid line are slower than the experimental ones. 5.2. Effect of the Non-Constant Input Shear Rate History In Figures 12 and 13, we explore the consequences of shear rate history, with all the parameters kept the same as in Figure 10. Rather than the actual shear rate history in experiments described by Equation (21), in Figure 13 we consider a constant shear rate and in Figure 12, we use a relaxational time history of shear rate without overshoot. This turns out to have relevant consequences. As we can see in Figure 13, the shapes of the curves change significantly upon using a constant shear rate input. Two most noticeable changes are, the small difference between the peak stress of different shear rate curves, and the very short yield times the model has predicted. This means that the fitness of previous models considered in the literature might be affected by the shear rate history, especially when the same type of rheometer (controlled-stress rheometer) is used. One problem of using the Maxwell-type stress equation is that the stress σ would have a non-zero initial value at time 0 s, as shown in Figure 13b, which however, does not agree with the experiments (discrete black points in Figure 8). With the inertial delay taken into account in Equation (21), this is no longer a limitation of the Maxwell-type model. The curves in Figure 12a and the region below the stress peak in Figure 12b are satisfactory, whereas the decreasing rates of shear stress σ for both the fitting and predicted curves (plotted with solid lines) in Figure 12b are slightly higher than the corresponding experimental ones (plotted with dash lines). Also, the decreasing rate of the consequent curves at both shear rates are not as closer to the experimental one as those in Figure 10b. This comparison shows influence of the overshoot of input on the predicted curves. (a) (b) Figure 12. Model predictions of shear stress as a function of shear strain and of time using a relaxational time history of shear rate Equation (21) (with the fitting values in Table 2: (a) shear stress as a function of shear strain; (b) shear stress as a function of time). Energies 2021, 14, 536 20 of 25 (a) (b) Figure 13. Model predictions of shear stress as a function of shear strain and of time using constant shear rate γ̇ without inertial delay: (a) shear stress as a function of shear strain; (b) shear stress as a function of time. 5.3. Marrucci Model vs. Maxwell Model In Figure 14, the curves are also obtained by using the same material functions Equations (14) and (10) as in Figure 10, and the same input shear rate history without over- shoot described by function Equation (16), but instead of Marrucci’s non-linear viscoelastic model we have used the Maxwell viscoelastic equation σ + (η/G )(dσ/dt) = η γ̇. We use a different parameterization from that used to obtain Figure 10, by taking the value of equilibrium modulus G0,1 = 700 Pa, p2 = 1.0, β = 0.17 and others kept the same. As a result, although the yield stress for the fitting curve of 0.1 s−1 is of the same order of the experiment alone, the height of the predicted curve of 1 s−1 is far from the experimental ones. Moreover, the peak heights of the fitting and predicted curves (plotted with solid lines) are not to the scale as the experimental curves (plotted with dash lines). In Figure 14, the decreasing rate after yielding obtained by Maxwell model are unrealistically slow. Comparison of Figures 10 and 14 shows the use of Marrucci’s model clearly improves the description of this situation. It is probably because the dynamics of shear modulus taken into consideration in Marrucci’s model has incorporated the thixotropic behavior of waxy gel and makes the description closer to its true mechanics. (a) (b) Figure 14. Predictions of shear stress as a function of shear strain and of time with Maxwell viscoelastic model for stress instead of Marrucci model Equation (17): (a) shear stress as a function of shear strain; (b) shear stress as a function of time. Energies 2021, 14, 536 21 of 25 5.4. Effect of Fractal Dimension Another factor to be considered related to structure is the value of the fractal dimension D. Generally speaking, fractal dimension is a parameter that may characterize the feature of the gel structure. The deviation of the value of fractal dimension from the dimension of Euclidean space D = 3 may be used as a measure for the complexity of the gel structure, which depends on the mechanism of aggregation. The actual value of fractal dimension found for waxy crude gel is often between 1.5 and 2.5 [40,44]. In the prediction results shown in Section 4.4, as well as in the anylasis in Sections 5.1–5.3, we have taken it as D = 1.5, which is noted in the literature as possible value of D [44]; here we compare it with the results obtained for D = 2. From Figure 15, we see that with lower fractal dimension assumed, the blob radius as well as the stress after yielding decrease faster with strain than curves for D = 2. The curves in Figure 15b share with Figures 10b and 11b the steep uprise for γ between 1 and 10. (a) (b) Figure 15. Model predictions of shear stress as a function of shear strain and of time with taking fractal dimensions D = 1.5 and D = 2.0: (a) shear stress as a function of shear strain; (b) shear stress as a function of time. 6. Conclusions We have presented a simplified theoretical model for the yielding of a waxy crude oil by assuming the gel as a fractal structure. The predicted results of the model show a very close form in comparison with experimental curves. It gives a quantitative description of the isothermal nonlinear viscoelastic behavior of waxy crude gels formed under a given constant cooling rate, and may be used to estimate the minimum pressure and time needed to restart the flow after a period of shutdown. It has sufficient simplicity and wide flexibility to explore the influence of different physical factors, such as fractal dimension, shear rate history or the initial elastic modulus. Moreover, since we have not used any specific details of a waxy crude oil in the modeling process, the main framework of our model is sufficiently general to be also applicable to the yielding of other fractal gels. Our model is based on evolution equations for stress σ and for structural parameter λ, and on a modelization for G (γ) and η (γ, λ, γ̇). It starts from a definite physical model for the structure parameter λ instead of from more abstract parameterizations. Compared to previous works, five novel factors especially considered in our work have consequent advantages, namely: (1) the use of Marrucci’s nonlinear viscoelastic model in descriptions for the stress evolution; (2) The use of an explicit microscopic interpretation for the structure parameter λ in terms of the radius RF of the gel clusters, which allows us to estimate a range of the yield strain with simple arguments; (3) the material function for shear modulus G inspired by the double-network model; (4) the application of well-established sophisticated viscosity model to the description of the viscosity η of a transient process; and (5) the Energies 2021, 14, 536 22 of 25 modification the actual imposed shear rate history γ̇(t) which deviates from the setting one due to device inertia. Though the main framework of the model is very general, we have focused our attention on a concrete application. Further studies should be the following ones: (1) An extension to more values of nominal shear rate should be examined. Furthermore, applications to other typical rheometric situations such as oscillating shear and small-step increase of strain should also be carried out. (2) Physical grounds for the inertia effect of device as well as the double-network structure of elastic modulus still need further exploration. The empirical function Equation (21) for the actual shear rate history and the material function Equation (18) for the elastic modulus should be improved accordingly. (3) Based on the presented model developed for the yielding region and applied to a fixed temperature history and composition situation, a comprehensive description of the rheological behavior may be extended to. First, the composition and temperature dependence of the material functions should be studied to compare the predictions of the model in different circumstances. Secondly, the structural evolution needs to be understood microscopically, leading to a more detailed evolution equation on the basis of a rate function incorporating the build-up and destruction of structure. Author Contributions: Conceptualization, M.S.; Formal analysis, M.S., D.J. and Z.W.; Funding acquisition, M.S.; Methodology, M.S. and D.J.; Software, M.S. and Z.W.; Validation, M.S. and Z.W.; Writing—original draft, M.S. and D.J.; Writing—review & editing, M.S., D.J. and Z.W. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by the Foundation of China University of Petroleum (Beijing) for Start-up Research Programs grant number 2462017YJRC021, and the Strategic Cooperation Technology Projects of CNPC and CUPB grant number ZLZX2020-05. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: The data presented in this study are available in the article. Acknowledgments: We would like to thank Langhui Zhang for the stress-strain/-time data in Figure 1 and the viscosity-temperature data in Table 1. Conflicts of Interest: The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. Abbreviations The following abbreviations are used in this manuscript: DSC Differncial Scanning Calormetry WAT Waxy Appearance Temperature SAOS Small amplitude oscillatory shear Energies 2021, 14, 536 23 of 25 Appendix A. Estimation of the Yield Strain with Analytical Geometry Figure A1. Geometric graph of elongated of blobs. We will use analytical geometry method to estimate the strain at the fracture (See Figure A1). The circle in the coordinate represents one of the blobs with its origi- nal shape and location before deformation, and the closed curve represents its status when flow begins. The origin of the coordinate is where the undeformed blob contacts its under- neath neighboring blob. When a shear rate equal to γ̇ is imposed on the blobs, the upper part will move faster than the lower part due to the gradient of velocity. Point O(0, 0) does not move during deformation, whereas point A (0, 2RF0 ), the top, moves ∆x = γRF . If the horizontal displacement of each point on circle is assumed to be proportional to its y-axis value, it will be translated with respect to the line y = x/γ. Assuming P(x, y) is any point on the circle, then we have ( x − γy)2 + (y − RF0 )2 = R2F0 (A1) Equation (A1), which can be identified as an elliptic equation, would be the equation of the blob after deformation. If we move the center of the ellipse to the origin of the coordinate, its equation will become ( x − γy)2 + y2 = R2F0 (A2) Note that the equation of an origin-centered ellipse, with semi-major axis l1 and semi-minor l2 , rotated by angle θ 0 has the following normalized form ( x cos θ 0 − y sin θ 0 )2 ( x sin θ 0 − y cos θ 0 )2 + =1 (A3) l12 l22 Equations (A2) and (A3) may be respectively rewritten as x2 y2 γ2 2 2 + 2 2 − 2xy( ) =1 (A4) RF0 RF / (1 + γ ) RF0 cos2 θ 0 sin2 θ 0 2 0 2 sin θ cos2 θ 0 1 1 x2 ( 2 + 2 ) + y ( 2 + 2 ) − 2xy sin θ 0 cos θ 0 ( 2 + 2 ) = 1 (A5) l1 l2 l1 l2 l1 l2 By comparing the correspondent coefficient in Equations (A4) and (A5), one gets the following equation set Energies 2021, 14, 536 24 of 25 cos2 θ 0 sin2 θ 0  1 + = 2     l 2 l22 RF0  sin21 θ 0   cos2 θ 0 1 + γ2  + = (A6)   l12 l22 R2F0 γ2 2   1 1  sin θ 0 cos θ 0 ( 2 + 2 ) = (  )   l1 l2 RF0 With the relation sin2 2θ 0 + cos2 2θ 0 = 1 and the restriction condition l2 ≤ 0.93Rt extrmF0, we eliminate l1 , l2 and RF0 in Equation (A6), and find the equation of γy as 3γy8 − 4.86γy6 + 2.76γy4 + 6.10γy2 − 0.38 = 0 (A7) Since the strain is rather small by the time of yielding, we neglect the high-order infinitesimals γy8 and γy6 , thus obtaining the approximate estimation of the critical shear strain γy = 0.25. References 1. Rønningsen, H.P. Rheological behaviour of gelled, waxy North Sea crude oils. J. Pet. Sci. Eng. 1992, 7, 177–213. [CrossRef] 2. Davidson, M.R.; Nguyen, Q.D.; Chang, C.; Rønningsen, H.P. A model for restart of a pipeline with compressible gelled waxy crude oil. J. Non-Newton. Fluid Mech. 2004, 123, 269–280. [CrossRef] 3. Vinay, G.; Wachs, A.; Frigaard, I. Start-up transients and efficient computation of isothermal waxy crude oil flows. J. Non-Newton. Fluid Mech. 2007, 143, 141–156. [CrossRef] 4. de Souza Mendes, P.R.; de Abreu Soaresa, F.S.; Zigliob, C.M.; Gonçalves, M. Startup flow of gelled crudes in pipelines. J. Non-Newtonian Fluid Mech. 2012, 179–180, 23–31. [CrossRef] 5. Cheng, D.C.H. Yield stress: A time-dependent property and how to measure it. Rheol. Acta 1986, 25, 542–554. [CrossRef] 6. Wardhaugh, L.T.; Boger, D.V. The measurement and description of the yielding behavior of waxy crude oil. J. Rheol. 1991, 35, 1121–1156. [CrossRef] 7. Chang, C.; Boger, D.V.; Nguyen, Q.D. The yielding of waxy crude oils. Ind. Eng. Chem. Res. 1998, 37, 1551–1559. [CrossRef] 8. Kane, M.; Djabourova, M.; Volle, J.L. Rheology and structure of waxy crude oils in quiescent and under shearing conditions. Fuel 2004, 83, 1591–1605. [CrossRef] 9. Visintin, R.F.G.; Lapasin, R.; Vignati, E.; D’Antona, P.; Lockhart, T.P. Rheological Behavior and Structural Interpretation of Waxy Crude Oil Gels. Langmuir 2005, 21, 6240–6249. [CrossRef] 10. Andrade, D.E.V.; da Cruz, A.C.B.; Franco, A.T.; Negrão, C.O.R. Influence of the initial cooling temperature on the gelation and yield stress of waxy crude oils. Rheol. Acta 2015, 54, 149–157. [CrossRef] 11. Tarcha, B.A.; Forte, B.P.P.; Soares, E.J.; Thompson, R.L. Critical quantities on the yielding process of waxy crude oils. Rheol. Acta 2015, 54, 479–499. [CrossRef] 12. Zhang, L. The Yield Properties of Gelled Waxy Crude Oil. Master’s Thesis, China University of Petroleum (Beijing), Bei- jing, China, 2011. 13. Houska, M. Inzenyrske Aspekty Reologie Tixotropnich Kapalin. Ph.D. Thesis, Techn. Univ. Prague, Prague, Czech Republic, 1980. 14. Doraiswamy, D.; Mujumdar, A.N.; Tsao, I.L.; Beris, A.N.; Danforth, S.C.; Metzner, A.B. The CoxMerz rule extended: Rheological model for concentrated suspensions. J. Rheol. 1991, 35, 647–685. [CrossRef] 15. Jia, B.; Zhang, J. Yield behavior of waxy crude gel: Effect of isothermal structure development before prior applied stress. Ind. Eng. Chem. Res. 2012, 2012, 10977–10982. [CrossRef] 16. Teng, H.; Zhang, J. A new thixotropic model for waxy crude. Rheol. Acta 2013, 52, 903–911. [CrossRef] 17. Dullaert, K.; Mewis, J. A structural kinetics model for thixotropy. J. Non-Newton. Fluid Mech. 2006, 139, 21–30. [CrossRef] 18. Teng, H.; Zhang, J. Rheological Behavior and Structural Interpretation of Waxy Crude Oil Gels. Ind. Eng. Chem. Res. 2013, 52, 8079–8089. [CrossRef] 19. Mujumdar, A.; Beris, A.N.; Metzner, A.B. Transient phenomena in thixotropic systems. Ind. Eng. Chem. Res 2002, 102, 157–178. [CrossRef] 20. Dimitrious, C.J.; McKinley, G.H. A comprehensive constitutive law for waxy crude oil: A thixotropic yield stress fluid. Soft Matter 2014, 35, 6591–6858. [CrossRef] 21. Mewis, J.; Wagner, N.J. Thixotropy. Adv. Colloid Interface Sci. 2009, 147–148, 214–227. [CrossRef] 22. Coussot, P.; Leonov, A.I.; Piau, J.M. Rheology of concentrated dispersed systems in a low molecular weight matrix. J. Rheol. 1993, 46, 179–217. doi:10.1016/0377-0257(93)85046-D. [CrossRef] 23. Yziquel, F.; Carreau, P.J.; Moan, M.; Tanguy, P.A. Rheological modeling of concentrated colloidal suspensions. J. Non-Newton. Fluid Mech. 1999, 86, 133–155. [CrossRef] Energies 2021, 14, 536 25 of 25 24. Quemada, D. Rheological modelling of complex fluids: IV: Thixotropic and thixoelastic behaviour. Start-up and stress relaxation, creep tests and hysteresis cycles. Eur. Phys. J. Appl. Phys. 1999, 86, 191–207. doi:10.1051/epjap:1999128. [CrossRef] 25. de Souza Mendes, P.R. Modeling the thixotropic behavior of structured fluids. J. Non-Newton. Fluid Mech. 2009, 164, 66–75. [CrossRef] 26. de Souza Mendes, P.R. Thixotropic elasto-viscoplastic model for structured fluids. Soft Matter 2011, 7, 2471–2483. doi:10.1039/c0sm01021a. [CrossRef] 27. de Souza Mendes, P.R.; Thompson, R.L. A unified approach to model elasto-viscoplastic thixotropic yield-stress materials and apparent yield-stress fluids. Rheol. Acta 2013, 52, 673–694. doi:10.1007/s00397-013-0699-1. [CrossRef] 28. Geest, V.C.B.G.C.V.D.; Merino-Garcia, D.; Bannwart, A.C. A modified elasto-viscoplastic thixotropic model for two commercial gelled waxy crude oils. Rheol. Acta 2015, 54, 545–561. doi:10.1007/s00397-015-0849-8. [CrossRef] 29. Mendes, R.; Vinay, G.; Ovarlez, G.; Coussot, P. Modeling the rheological behavior of waxy crude oils as function of flow and temperature history. J. Rheol. 2015, 59, 703–732. doi:10.1122/1.4916531. [CrossRef] 30. Marrucci, G. The Free Energy Constitutive Equation for Polymer Solutions from the Dumbbell Model. Trans. Soc. Rheol. 1972, 16, 321–300. [CrossRef] 31. Marrucci, G.; Titomanlio, G.; Sarti, G.C. Testing of a constitutive equation for entangled networks by elongational and shear data of polymer melts. Rheol. Acta 1973, 12, 269–275. [CrossRef] 32. Acierno, D.; Mantia, F.P.L.; Marrucci, G.; Titomanlio, G. A non-linear viscoelastic model with structure-dependent relaxation times. I. Basic formulation. J. Non-Newton. Fluid Mech. 1976, 1, 125–146. [CrossRef] 33. Acierno, D.; Mantia, F.P.L.; Marrucci, G.; Rizzo, G.; Titomanlio, G. A non-linear viscoelastic model with structure-dependent relaxation times. II. Comparison with l. D. polyethylene transient stress results. J. Non-Newton. Fluid Mech. 1976, 1, 147–157. [CrossRef] 34. Acierno, D.; Mantia, F.P.L.; Marrucci, G. A non-linear viscoelastic model with structure-dependent relaxation times. II. Comparison with l. D. polyethylene transient stress results. J. Non-Newton. Fluid Mech. 1977, 2, 271–280. [CrossRef] 35. Snabre, P.; Mills, P. I. Rheology of weakly flocculated Suspensions of rigid particles. J. Phys. III Fr. 1996, 6, 1811–1834.:1996215. [CrossRef] 36. Snabre, P.; Mills, P. II. Rheology of weakly flocculated Suspensions of viscoelastic particles. J. Phys. III Fr. 1996, 6, 1835– 1855.:1996216. [CrossRef] 37. Dimitrous, C.J.; McKinley, G.H.; Venkatesan, R. Rheo-piv analysis of the yielding and flow of model waxy crude oils. Energy Fuels 2011, 25, 3040–3052. doi:10.1021/ef2002348. [CrossRef] 38. Bautista, F.; de Santos, J.M.; Puig, J.E.; Manero, O. Understanding thixotropic and antithixotropic behavior of viscoelastic micellar solutions and liquid crystalline dispersions. I. The model. J. Non-Newton. Fluid Mech. 1999, 80, 93–113. [CrossRef] 39. Lorge, O.; Djabourov, M.; Brucy, F. Crystallisation and gelation of waxy crude oils under flowing conditions. Rev. L’institut Français Pétrole 1997, 52, 235–239. [CrossRef] 40. Gao, P.; Zhang, J.; Ma, G. Direct image-based fractal characterization of morphologies and structures of wax crystals in waxy crude oils. J. Phys. Condens. Matter 2006, 28, 11487–11506. [CrossRef] 41. Hou, L. Experimental study on yield behavior of Daqing crude oil. Rheol. Acta 2012, 51, 603–607. [CrossRef] 42. Venkatesan, R.; Ostlund, J.; Chawla, H.; Wattana, P.; Nyden, M.; Fogler, H.S. The effect of asphaltenes on the gelation of waxy oils. Energy Fuels 2003, 14, 1630–1640. [CrossRef] 43. Soskey, P.R.; Winter, H.H. Equibiaxial Extension of Two Polymer Melts: Polystyrene and Low Density Polyethylene. J. Rheol. 1985, 29, 493–517. [CrossRef] 44. Yi, S.; Zhang, J. Shear-induced change in morphology of wax crystals and flow properties of waxy crudes modified with the pour-point depressant. Energy Fuels 2011, 25, 5660–5671. [CrossRef]