Influence of thixotropic parameters on a non-Newtonian fluid flow

Thixotropy is an important characteristic of complex fluids. It is a time-dependent rheological property associated with the microstructural evolution of the materials during their flowing process. Indeed, two counter transition mechanisms are considered to take place simultaneously: the recovery and the breakdown; the former makes the materials more solid, while the latter makes them more liquid. Parameters that control the rate of these mechanisms are called thixotropic parameters. These fluids can be found in mud streams and lava flows in nature and are frequently encountered in engineering systems, such as fresh cement, ink, polymers, crude oil, pharmaceutical products, and food additives. In addition to thixotropy, these fluids can possess other rheological characteristics, such as viscoplasticity and viscoelasticity.

There have been few studies of viscoplastic fluid flows over an obstacle in the literature. Indeed, Mossaz et al.^1,2 simulated a nonthixotropic yield stress fluid flow over a circular cylinder at several (moderate) Reynolds and Oldroyd numbers. The simulation was validated using experimental data of Mossaz et al.³ Recently, Bertevas et al.⁴ adopted a two phase mixture model in the framework of Smoothed Particle Hydrodynamics (SPH) to simulate turbulent sediment transport. It was applied to a problem of sediment bed–plate interaction, in which only the yield stress characteristic was taken into account using the Bingham–Papanastasiou model, though experimental data of the sediment suspension exhibited complex behaviors with a yield stress, viscoelasticity, and thixotropy. Thixotropy was not considered, as the problem was limited to low clay concentrations. Likewise, Salmanidou et al.⁵ and Ren et al.⁶ investigated submarine landslide kinematics and tsunami; the Bingham model was employed in the former for the muddy debris flow, and the Herschel–Bulkley (HB) model in the latter as the material contained mostly small blocks and soil.

It is good to differentiate thixotropy from viscoelastic aging, which can be based on the relaxation of stress after the flow is ended, i.e., it is rapid for thixotropic liquids, and slow and continues to get slower for viscoelastic aging ones. A good description about this can be found in the work of Larson and Wei.⁷ Agarwal et al.⁸ investigated the possibility of aging in an aqueous dispersion of Carbopol 940 microgel. Results showed an increase in the viscosity (or stress) overshoot with an increase in the resting time, implying a possible aging under quiescent conditions. However, as the yield stress was found to remain constant, an open question was incurred that whether Carbopol is a simple yield stress fluid, or a thixotropic yield stress one, or a thixotropic constant yield stress one, i.e., a so-called third category.

Here, we focus on thixotropy only. To phenomenologically model thixotropic liquids, e.g., ceramic pastes, Moore⁹ introduced a transport kinetic equation for a structural parameter representing the microstructure state. Burgos et al.¹⁰ modified it by adding an exponential factor to the breakdown term for semisolid suspensions. More recently, Derksen¹¹ and Syrakos et al.¹² adopted Moore’s equation flow numerical simulations of thixotropic flows, and both showed that the flow patterns were affected greatly by the yield stress and thixotropy.

In this work, a brief review on thixotropic models is first carried out. Some of them are greatly different from one another in terms of how the breakdown and build-up processes are modeled. As an attempt to better understand to what extent the thixotropic parameters influence the flow behaviors of a thixotropic fluid past a cylinder, a parametric study is carried out wherein Burgo’s model is employed.

The rest of the paper is organized as follows: Governing equations, model validation for viscoplastic fluids, and a brief review on the thixotropic models are given in Sec. II; Sec. III presents results for the flow field in a wide range of thixotropic parameters; and remarks are made in Sec. IV.

A fluid flow is governed by a continuity equation and a momentum equation as follows:

where u and f are the velocity and body force vectors, respectively; ρ is the fluid density; σ is the total stress tensor with p being the pressure, I the unit tensor, τ the deviatoric stress tensor, and $γ̇$ the deformation rate tensor.

For a yield stress fluid with τ₀ being the yield stress, the Herschel–Bulkley (HB) model can be employed,

Here, K and n are the consistency and the power index, respectively; $γ̇$ is the magnitude of $γ̇$⁠. When n = 1, this model becomes a Bingham one, which is mainly investigated in this work.

Reynolds and Bingham numbers characterizing a Bingham fluid flow are defined, respectively, as

where u_∞ is the far-field velocity and D is the cylinder’s diameter.

To tackle the discontinuity of the HB or Bingham model, Papanastasiou’s regularization¹³ is utilized as

with M being the regularization parameter. For the sake of computational efficiency, M = 40 000 is chosen here. Details about this can be found in the work of Bui and Ho.¹⁴ ANSYS Fluent is utilized for the computation; and a User-Defined Function (UDF) is developed to take Eq. (6) into account.

Figure 1 shows a circular computational domain employed here; its diameter is D_∞ with D_∞/D = 200. A structured mesh consisting of 92 000 elements is generated; mesh convergence has been verified and presented elsewhere.¹⁴ 

With the numerical model described above for nonthixotropic HB fluid, results for the streamline pattern and the size and shape of the unyielded zones are presented in Fig. 2 in comparison with those obtained numerically and experimentally by Mossaz et al.³ The flow was at Re = 45 and Od = 0.32. Good agreement is observed.

Thixotropy is a complex rheological property associated with the microstructure of a material, which evolves with flowing time. To model the microstructure evolution, a structural parameter λ ranging from 0 to 1 is introduced. The lower limit corresponds to a completely broken state, while the upper one to a fully recovered state of the microstructure. The evolution rate of λ can be expressed as follows:

Here, on the RHS, the first term is the recovery rate and the second term is the breakdown one; K₁ and K₂ are the rate parameters; and m and n are the rate indices. These indices have commonly been assumed to be 1, making Eq. (7) a first-order differential equation.^9,10,15–22 However, there exist other values found in the literature including m = 0,^23–25 m = 2,²⁶ and n = 2.²³ The recovery parameter K₁ can be a constant,^{9,10,15,16,23–25,27,28} a function of strain rate^{17,21,22,26,29} or a complex one of strain rate and time.^18–20 The breakdown parameter K₂ can be a linear function of strain rate,^{9,15–17,22,24,25,28} or a more complex one of strain rate,^10,21,26 e.g., exponential or power-law, strain rate, and time,^18–20 or strain rate and shear stress.^27,29 As a summary, a list of various thixotropic models is given in Table I.

It is evident that these thixotropic models are very much different for different materials. Therefore, special care needs to be taken when employing them. In an effort to study the effects of the thixotropic parameters on the flow behaviors, we choose to adopt the model by Burgos et al.¹⁰ It was proposed for semisolid metals; however, it could possibly be used for, e.g., fresh cement as both materials change their state from liquid to solid over time, though, the timescale can be different. Accordingly, the evolution of λ is written as

with α being the recovery parameter (i.e., K₁), and β₁ and β₂ being the breakdown parameter and exponent, respectively. Note that β₂ ≥ 0; when β₂ = 0, Eq. (8) becomes the first-order equation;^9,15,16,28 and when β₂ > 0, the breakdown rate is augmented.^10,30–32 These parameters are varied as α = 0.001–1, β₁ = 0.001–1, and β₂ = 0.0005–1. The range of α and β₁ typically spans from Newtonian to nonthixotropic Bingham fluid or vice versa. Furthermore, it is assumed that the yield stress is determined as τ₀ = λτ_y with τ_y being its maximum value, i.e., when the microstructure is fully developed or λ = 1.

An interesting characteristic of thixotropic fluids is stress overshoot when they are subjected to a sudden jump in the strain rate. For one-dimensional (1D) problem, from Eq. (3), the shear stress is calculated as $τ=Kγ̇+λτy$ for a Bingham fluid. When $γ̇1$ is suddenly increased to $γ̇2$ within a short time, λ would remain constant as the microstructure needs time to respond to this sudden change, the shear stress is therefore increased abruptly. As time goes on, the material is broken, which eases the stress and results in a smaller value of λ; a steady state is approached after a certain time. Likewise, a stress undershoot corresponds to a sudden decrease of the shear rate. The manifestation of the stress over/undershoot would depend on how fast and how much the strain rate is changed. For it to be better observed, or to have any important effects on the flow behaviors, with a small change of the strain rate, a short time interval may be required. The latter can be in the order of milliseconds or below. A gradual change of the strain rate would not incur the stress over/undershoot phenomenon. Another timescale associated with this phenomenon is the time for the stress after the strain rate has changed to approach its steady state value (or 95% of this). As reported by Skadsem et al.²⁰ for water- and oil-based drilling fluids, it was in the range of 1–7 s depending on the magnitude of the strain rate change.

In this work, we investigate a thixotropic flow past a circular cylinder. It is now convenient to look at the moving fluid from the Lagrangian frame of reference. As a parcel of fluid moves through different areas, it may experience different strain rates. Therefore, its microstructure is constantly broken and/or recovered, and a steady state would never be attained. It is attained only when the fluid parcel passes regions of a uniform strain rate for a long enough time, that is, beyond the yield line. Here, the flow is over a circular cylinder at a relatively small Reynolds number, i.e., Re = 45, and thus the moving speed is small and the variation of the strain rate is moderate as illustrated in Fig. 3. It is assumed that the fluid does not experience any abrupt change of the strain rate, and hence the over/undershoot phenomenon can be neglected.

It is good to distinguish the steady state for a moving fluid parcel with an equilibrium one for the whole domain defined instead in the Eulerian frame. The latter is attained when all the flow and fluid properties including the microstructure at any location in the flow field do not change with time. In what follows, results are presented and analyzed at this state.

In this part, simulations for α = 0.001, 0.05, and 1 are carried out. The breakdown parameters are kept constant: β₁ = 0.05 and β₂ = 0.05. Moreover, in this whole paper, Re = 45 and Bn = 0.5; it is noted that the latter is based on τ_y, the yield stress when λ = 1. Results for the streamline pattern and the distribution of λ are shown in Fig. 4.

It can be observed that when the recovery parameter α is small [Fig. 4(a)], the structural parameter λ is found to be small, approaching zero in the recirculation bubble behind the cylinder, i.e., the structure is almost completely broken. This is because the recovery process is weak, and dominated by the breakdown one. When α is great [Fig. 4(c)], λ is great especially in the cylinder’s wake, where λ → 1 and the structure is fully recovered. The recirculation bubble is found to be smaller for greater α. It is noted that when λ = 0, τ_y = 0 and the fluid becomes Newtonian; when λ is great, the yield stress and thus the effective viscosity, the term inside the brackets of Eq. (3) or Eq. (6), are great, which makes the flow more laminar. As a consequence, the recirculation bubble is smaller for greater α.

In addition, the distribution of λ along a horizontal line passing the cylinder’s center, i.e., the x-centerline, is shown in Fig. 5. As can be seen, λ drops down at some distance before and after the cylinder. In general, λ is smaller in the wake than in front of the cylinder except for the case α = 1 where the broken material in the front is quickly recovered when it moves downstream.

Three characteristic lengths are introduced as in Fig. 6. L₁ and L₃ are distances from the cylinder to a point where λ = 0.99, considered to correspond to a fully structured state. L₂ is the distance from the cylinder to the minimum point of λ in the wake, from where the recovery process predominantly takes place. Results for the characteristic lengths are shown in Table II. As can be noticed, the greater the recovery parameter α, the shorter L₁ and L₃. Moreover, although L₂ is relatively small, it seems to depend on a combination of α, the extent to which the material is broken before the cylinder and/or the flow structure behind it.

Furthermore, unyielded zones defined by τ ≤ τ₀ can be classified into three main types: the far-field one, the static one attached to the cylinder, and the moving one scattered in the flow field.^1–3 Figure 7 shows the yield line/boundary beyond which the fluid is unyielded, i.e., the far-field zone. This boundary is found to be slightly smaller (in the wake) for greater values of α. When α = 1, it approaches that of an equivalent nonthixotropic Bingham flow.

The static and moving unyielded zones are illustrated in Fig. 8. As can be observed, these zones are smaller for smaller α; some even disappear when α = 0.001. The static zones are found at the rear and two lateral positions on the cylinder. The latter tend to move downstream when α is increased, corresponding to a smaller recirculation bubble. When α = 1, like the far-field one, the near-field unyielded zones are almost identical to those of the Bingham flow.

Additionally, Fig. 9 shows the drag coefficient, C_d, on the cylinder. It is found to slightly increase with the recovery parameter, and approach that of the Bingham flow when α = 1.

Here, the breakdown parameter is varied, that is, β₁ = 0.001, 0.05, and 1; the other parameters are kept constant as Re = 45, Bn = 0.5, α = 0.05, and β₂ = 0.05. Results for the streamline pattern and the distribution of λ are presented in Fig. 10. It can be noticed that the effect of β₁ on these results is opposite to that of α. Indeed, when β₁ is small, which corresponds to a weak breakdown process, the material is well structured, as in Fig. 10(a) for β₁ = 0.001, λ ≈ 1 in almost the entire domain. On the contrary, when β₁ is great, λ is small, for instance, when β₁ = 1 [Fig. 10(c)], the material near the cylinder is almost completely broken. In addition, the recirculation bubble is larger for greater values of β₁.

Figure 11 shows the distribution of λ along the x-centerline. It is noticed that λ decreases significantly in front of the cylinder when β₁ is great, e.g., λ → 1 when β₁ = 1. Results for the characteristic lengths are provided in Table III. L₁ and L₃ are found to increase when β₁ is increased from 0.001 to 1. Typically, when β₁ = 1, full recovery takes place even outside of the computational domain (L₃/D > 99.5).

Moreover, the yield line is shown in Fig. 12. As can be observed, it is smaller (in the wake) for smaller β₁. The yield line for β₁ = 0.001 approaches that of the equivalent nonthixotropic Bingham flow. The same is observed for the near-field unyielded zones as shown in Fig. 13. In contrast, when β₁ = 1, these zones shrink substantially or even disappear.

Furthermore, the result for C_d is presented in Fig. 14. It is obvious that C_d decreases with the increasing β₁. When β₁ is as great as 1, however, C_d is larger than that of the Newtonian flow, which implies that the effect of the yield stress is still important. When β₁ = 0.001, like the yield line and the near-field rigid zones, the drag coefficient approaches that of the nonthixotropic Bingham flow.

In this part, the influence of β₂ is examined. To this end, simulations are realized for β₂ = 0.0005, 0.05, and 1. Other parameters are unchanged, i.e., Re = 45, Bn = 0.5, α = 0.05, and β₁ = 0.05. Results for the streamline pattern and the distribution of λ are shown in Fig. 15.

It is noteworthy that the exponential factor $exp(β2γ̇)→1$ when β₂ → 0 and/or $γ̇→0$⁠, and thus its contribution to the breakdown rate becomes negligible; the model of Eq. (8) becomes that of Moore⁹ and Toorman.¹⁵ Obviously, the difference in the flow field for β₂ = 0.0005 [Figs. 15(a)] and β₂ = 0.05 (two orders of magnitude greater) [Fig. 15(b)] is insignificant; and it becomes important only when β₂ is increased up to 1 [Fig. 15(c)]. Moreover, the difference is mainly in the lateral regions of the cylinder where the strain rate is supposed to be large, and hence the contribution of the exponential factor is significant.

The distribution of λ along the x-centerline is shown in Fig. 16, and the characteristic lengths in Table IV. As can be seen, the effect of β₂ is insignificant when β₂ is small. It becomes important when β₂ is as great as 1.

Moreover, Fig. 17 presents the yield line for various values of β₂; its relative difference is found to be negligible. The static and moving unyielded zones are shown in Fig. 18. The two static zones on the lateral sides are found to shift slightly downstream when β₂ is decreased. In addition, some of the moving zones disappear when β₂ is as large as 1.

Additionally, Fig. 19 shows the result for C_d on the cylinder as a function of β₂. It is found to be almost the same when β₂ ≤ 0.05, and it slightly decreases when β₂ is increased. The decrease is noticeable only when β₂ ≥ 0.1.

In this work, a brief review on the thixotropic models has been carried out. The influence of the parameters controlling the microstructural recovery α and breakdown (β₁ and β₂) rates on a thixotropic flow past a cylinder has been investigated. Simulations with varying α, β₁, and β₂ at Re = 45 and Bn = 0.5 were performed. Results for various flow characteristics including the yielded and unyielded zones, the structure state represented by the structural parameter λ, and the drag force coefficient on the cylinder were presented and discussed.

The results showed that when the recovery parameter α is increased, the fluid becomes more structured, and the wake is relatively narrower with smaller recirculation bubble; in addition, the near-field unyielded zones become larger, and the drag coefficient on the cylinder becomes greater, approaching that of the equivalent nonthixotropic Bingham flow. On the contrary, these trends were observed when the breakdown parameter β₁ and exponent β₂ were decreased. With the same order of magnitude, the influence of β₁ was found to be more pronounced than that of β₂, which was significant only when β₂ ≥ 0.1.

For future research, experimental characterization of the thixotropic parameters for, e.g., fresh cement and clay would be realized. Furthermore, the overshoot and undershoot phenomena would be taken into account for the estimation of the hydrodynamic forces.

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

This research was funded by Ho Chi Minh City’s Department of Science and Technology (HCMC-DOST) and Institute for Computational Science and Technology (ICST) under Grant No. 08/2018/HD-KHCNTT.