Non-Newtonian fluids exhibiting complex rheological characteristics, such as yield stress and thixotropy, are frequently encountered in nature and industries. Thixotropy is a time-dependent shear thinning property, associated with the microstructural evolution of materials. During a flowing process, two microstructure 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. The microstructure is characterized by a dimensionless structural parameter, whose evolution is modeled by a rate equation consisting of two terms representing the rate of the two mechanisms. A brief review on thixotropic models for different materials is first carried out. It is then assumed that the recovery rate depends linearly on the structural parameter, and the breakdown one is a complex function of it and the shear rate. This work aims at investigating the influence of the parameters that control the recovery and breakdown rates on the flow of a thixotropic fluid past a circular cylinder. In addition, the Bingham and/or Herschel–Bulkley model with Papanastasiou’s regularization is utilized. Various flow characteristics, such as the microstructure evolution and the flow field including the yielded and unyielded zones, are analyzed and discussed in detail. The simulation results show that the size and shape of both static and moving unyielded zones are considerably affected by the thixotropic parameters.

## I. INTRODUCTION

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.*^{3} Recently, Bertevas *et al.*^{4} 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.*^{5} and Ren *et al.*^{6} 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.^{7} Agarwal *et al.*^{8} 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^{9} introduced a transport kinetic equation for a structural parameter representing the microstructure state. Burgos *et al.*^{10} modified it by adding an exponential factor to the breakdown term for semisolid suspensions. More recently, Derksen^{11} and Syrakos *et al.*^{12} 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.

## II. METHODOLOGY

### A. Governing equations

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 $\gamma \u0307$ the deformation rate tensor.**

*τ*For a yield stress fluid with *τ*_{0} being the yield stress, the Herschel–Bulkley (HB) model can be employed,

Here, *K* and *n* are the consistency and the power index, respectively; $\gamma \u0307$ is the magnitude of $\gamma \u0307$. 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.

### B. Numerical approach

To tackle the discontinuity of the HB or Bingham model, Papanastasiou’s regularization^{13} 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.^{14} 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.^{14}

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.*^{3} The flow was at Re = 45 and Od = 0.32. Good agreement is observed.

### C. Thixotropic models

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*_{1} and *K*_{2} 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,^{26} and *n* = 2.^{23} The recovery parameter *K*_{1} 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*_{2} 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.

References . | Materials . | m
. | n
. | K_{1}
. | K_{2}
. |
---|---|---|---|---|---|

Moore^{9} | Ceramic pastes | ||||

Toorman^{15} | Cohesive sediment suspension | 1 | 1 | α | $\beta \gamma \u0307$ |

Armstrong et al.^{16} and Horner et al.^{28} | Human blood | ||||

Pinder^{23} | Hydrate slurry | 0 | 2 | α | Β |

Worrall and Tulliani^{17} | Clay suspension | 1 | 1 | $\alpha \gamma \u0307$ | $\beta \gamma \u0307$ |

Burgos et al.^{10} | Semisolid metal | 1 | 1 | α | $\beta \gamma \u0307ea\gamma \u0307$ |

Coussot et al.^{24} | Bentonite, Laponite, hair gel | ||||

Roussel^{25} | Fresh concretes | ||||

Yziquel et al.^{27} | Colloidal suspensions | 1 | 1 | α | $\beta \tau \gamma \u0307$ |

Dullaert and Mewis^{18} | Suspension of silica and carbon black particles | ||||

Mortazavi-Manesh and Shaw^{19} | Crude oil | 1 | 1 | $t\u2212b\alpha 1+\alpha 2\gamma \u03070.5$ | $t\u2212b\beta \gamma \u0307$ |

Skadsem et al.^{20} | Drilling fluid | ||||

Galindo-Rosales and Rubio-Hernandez^{21} | Bentonite | 1 | 1 | $108\gamma \u0307\u22123.4$ | $0.0016\gamma \u03070.44$ |

Alessandrini et al.^{26} | Clay and Kaolin suspension | 2 | 1 | $\alpha 1+\alpha 2\gamma \u0307p$ | $\beta \gamma \u0307q$ |

Zhu et al.^{22} | Colloidal gels | 1 | 1 | $\alpha 1\gamma \u0307+\alpha \beta \alpha 2\gamma \u0307p+1$ | $\beta \gamma \u0307$ |

Teng and Zhang^{29} | Waxy crude | 1 | 1 | $\alpha 1+\gamma p$ | $\beta (\tau \gamma \u0307)q1+\gamma p$ |

References . | Materials . | m
. | n
. | K_{1}
. | K_{2}
. |
---|---|---|---|---|---|

Moore^{9} | Ceramic pastes | ||||

Toorman^{15} | Cohesive sediment suspension | 1 | 1 | α | $\beta \gamma \u0307$ |

Armstrong et al.^{16} and Horner et al.^{28} | Human blood | ||||

Pinder^{23} | Hydrate slurry | 0 | 2 | α | Β |

Worrall and Tulliani^{17} | Clay suspension | 1 | 1 | $\alpha \gamma \u0307$ | $\beta \gamma \u0307$ |

Burgos et al.^{10} | Semisolid metal | 1 | 1 | α | $\beta \gamma \u0307ea\gamma \u0307$ |

Coussot et al.^{24} | Bentonite, Laponite, hair gel | ||||

Roussel^{25} | Fresh concretes | ||||

Yziquel et al.^{27} | Colloidal suspensions | 1 | 1 | α | $\beta \tau \gamma \u0307$ |

Dullaert and Mewis^{18} | Suspension of silica and carbon black particles | ||||

Mortazavi-Manesh and Shaw^{19} | Crude oil | 1 | 1 | $t\u2212b\alpha 1+\alpha 2\gamma \u03070.5$ | $t\u2212b\beta \gamma \u0307$ |

Skadsem et al.^{20} | Drilling fluid | ||||

Galindo-Rosales and Rubio-Hernandez^{21} | Bentonite | 1 | 1 | $108\gamma \u0307\u22123.4$ | $0.0016\gamma \u03070.44$ |

Alessandrini et al.^{26} | Clay and Kaolin suspension | 2 | 1 | $\alpha 1+\alpha 2\gamma \u0307p$ | $\beta \gamma \u0307q$ |

Zhu et al.^{22} | Colloidal gels | 1 | 1 | $\alpha 1\gamma \u0307+\alpha \beta \alpha 2\gamma \u0307p+1$ | $\beta \gamma \u0307$ |

Teng and Zhang^{29} | Waxy crude | 1 | 1 | $\alpha 1+\gamma p$ | $\beta (\tau \gamma \u0307)q1+\gamma p$ |

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

### D. Stress overshoot and undershoot

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 $\tau =K\gamma \u0307+\lambda \tau y$ for a Bingham fluid. When $\gamma \u03071$ is suddenly increased to $\gamma \u03072$ 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.*^{20} 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.

## III. RESULTS AND DISCUSSION

### A. Influence of the recovery parameter, $\alpha $

In this part, simulations for *α* = 0.001, 0.05, and 1 are carried out. The breakdown parameters are kept constant: *β*_{1} = 0.05 and *β*_{2} = 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*_{1} and *L*_{3} are distances from the cylinder to a point where *λ* = 0.99, considered to correspond to a fully structured state. *L*_{2} 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*_{1} and *L*_{3}. Moreover, although *L*_{2} 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.

α
. | L_{1}/D
. | L_{2}/D
. | L_{3}/D
. |
---|---|---|---|

0.001 | 1.799 | 1.633 | >99.5 |

0.05 | 1.686 | 3.736 | 80.209 |

1 | 1.032 | 2.399 | 5.976 |

α
. | L_{1}/D
. | L_{2}/D
. | L_{3}/D
. |
---|---|---|---|

0.001 | 1.799 | 1.633 | >99.5 |

0.05 | 1.686 | 3.736 | 80.209 |

1 | 1.032 | 2.399 | 5.976 |

Furthermore, unyielded zones defined by *τ* ≤ *τ*_{0} 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.

### B. Influence of the breakdown parameter, $\beta $_{1}

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

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

β_{1}
. | L_{1}/D
. | L_{2}/D
. | L_{3}/D
. |
---|---|---|---|

0.001 | 0.0202 | 2.399 | 10.14 |

0.05 | 1.686 | 3.736 | 80.209 |

1 | 7.667 | 3.092 | >99.5 |

β_{1}
. | L_{1}/D
. | L_{2}/D
. | L_{3}/D
. |
---|---|---|---|

0.001 | 0.0202 | 2.399 | 10.14 |

0.05 | 1.686 | 3.736 | 80.209 |

1 | 7.667 | 3.092 | >99.5 |

Moreover, the yield line is shown in Fig. 12. As can be observed, it is smaller (in the wake) for smaller *β*_{1}. The yield line for *β*_{1} = 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} = 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 *β*_{1}. When *β*_{1} 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 *β*_{1} = 0.001, like the yield line and the near-field rigid zones, the drag coefficient approaches that of the nonthixotropic Bingham flow.

### C. Influence of the breakdown exponent, $\beta $_{2}

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

It is noteworthy that the exponential factor $exp(\beta 2\gamma \u0307)\u21921$ when *β*_{2} → 0 and/or $\gamma \u0307\u21920$, and thus its contribution to the breakdown rate becomes negligible; the model of Eq. (8) becomes that of Moore^{9} and Toorman.^{15} Obviously, the difference in the flow field for *β*_{2} = 0.0005 [Figs. 15(a)] and *β*_{2} = 0.05 (two orders of magnitude greater) [Fig. 15(b)] is insignificant; and it becomes important only when *β*_{2} 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 *β*_{2} is insignificant when *β*_{2} is small. It becomes important when *β*_{2} is as great as 1.

β_{2}
. | L_{1}/D
. | L_{2}/D
. | L_{3}/D
. |
---|---|---|---|

0.0005 | 1.686 | 3.62 | 80.209 |

0.05 | 1.686 | 3.736 | 80.209 |

1 | 1.799 | 6.166 | 87.971 |

β_{2}
. | L_{1}/D
. | L_{2}/D
. | L_{3}/D
. |
---|---|---|---|

0.0005 | 1.686 | 3.62 | 80.209 |

0.05 | 1.686 | 3.736 | 80.209 |

1 | 1.799 | 6.166 | 87.971 |

Moreover, Fig. 17 presents the yield line for various values of *β*_{2}; 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 *β*_{2} is decreased. In addition, some of the moving zones disappear when *β*_{2} is as large as 1.

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

## IV. CONCLUDING REMARKS

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 (*β*_{1} and *β*_{2}) rates on a thixotropic flow past a cylinder has been investigated. Simulations with varying *α*, *β*_{1}, and *β*_{2} 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 *β*_{1} and exponent *β*_{2} were decreased. With the same order of magnitude, the influence of *β*_{1} was found to be more pronounced than that of *β*_{2}, which was significant only when *β*_{2} ≥ 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 AVAILABILITY

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

## ACKNOWLEDGMENTS

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.