Turbulent Pipe Flow of Thixotropic Fluids

In this study we consider the simplest representative scenario of fully developed turbulent flow of a single-phase incompressible and inelastic thixotropic fluid through a smooth long pipe. The flow domain consists of an axially-periodic straight pipe section of diameter $D$ and length $L_{z}=4\pi D$ with no-slip wall boundary conditions that is driven by a constant axial body force (Singh et al., 2018). Flow within the pipe is governed by the incompressible Navier-Stokes equations

where $\boldsymbol{v}$ is the fluid velocity, $\rho$ is density and $p$ static pressure, and the flow is driven by the constant axial body force

The shear stress tensor $\mathsfbi{\tau}=2\eta\mathsfbi{S}$ is the product of the rate of strain tensor $\mathsfbi{S}=1/2\left(\bnabla\boldsymbol{v}+\bnabla\boldsymbol{v}^{\top}\right)$ and the apparent viscosity

which is a function of local shear rate $\dot{\gamma}=\sqrt{2(\mathsfbi{S}:\mathsfbi{S})}$ and the structural parameter $\lambda$ that quantifies the impact of thixotropy upon the fluid rheology. A range of models for $\eta\left(\dot{\gamma},\lambda\right)$ have been proposed in the literature that span a wide range of materials (Burgos et al., 2001; Toorman, 1997; Farno et al., 2020; Mewis & Wagner, 2009; Barnes, 1999). As the fluid is inelastic, quantities such as the Deborah number and the thixoelastic number (Ewoldt & McKinley, 2017) are not relevant.

For the thixotropic evolution of $\lambda$ in a homogeneous flow, various kinetic models (Mujumdar et al., 2002; Mewis & Wagner, 2009; Larson, 2015; Barnes, 1997) have been proposed, most of which are of the form

where the kinetic functions $g(\dot{\gamma},\lambda)$ and $f(\dot{\gamma},\lambda)$ respectively represent the shear-induced breakdown and thermal and/or shear rebuild of the microstructure. Along with the shear viscosity $\eta$, determination of the kinetic functions form an integral part of the rheological characterization of thixotropic fluids. Under steady shear conditions ($\dot{\gamma}=\text{const.}$), these processes come into equilibrium and the structural parameter approaches its equilibrium value $\lambda\rightarrow\lambda_{\text{eq}}(\dot{\gamma})$ given by

where typically $\lambda_{\text{eq}}(\dot{\gamma})\rightarrow 1$ as $\dot{\gamma}\rightarrow 0$, $\lambda_{\text{eq}}(\dot{\gamma})\rightarrow 0$ as $\dot{\gamma}\rightarrow\infty$. In general, the viscosity $\eta(\dot{\gamma},\lambda)$ of thixotropic fluids varies strongly (typically decreasing) with both shear rate $\dot{\gamma}$ (shear thinning) and the structural parameter $\lambda$ (thixotropy) as the microstructure responds to imposed shear.

As this study aims to understand the fundamental mechanisms that govern thixotropic turbulence, we consider the simplest possible non-trivial thixotropic kinetic model, where the shear breakdown $g=k\,m\,\dot{\gamma}\,\lambda$ and thermal rebuild $f=m(1-\lambda)$ terms are simple linear functions of $\dot{\gamma}$ and $\lambda$, yet still capture the essential physics of more complex thixotropic models. Extension of (5) to heterogeneous flow such as turbulent pipe flow yields an advection-diffusion-reaction equation (ADRE) for the evolution of $\lambda$ as

where $\alpha$ characterizes the self-diffusivity of the microstructure which is typically very small, corresponding to Péclet number $\Pen\sim 10^{12}$ (Morris & Brady, 1996) for e.g. colloidal suspensions. However, a much larger value of $\alpha$ is typically used in simulations for numerical stability (Billingham & Ferguson, 1993). The thixotropic rate parameter $m$ governs the rate of convergence of $\lambda$ to the equilibrium value $\lambda_{\text{eq}}(\dot{\gamma})$. The parameter $k$ governs the relative rates of breakdown and rebuild, and impacts the equilibrium state as

which also corresponds to the limit of instantaneous thixotropic kinetics $m\rightarrow\infty$, where $\lambda\rightarrow\lambda(\dot{\gamma})$. Indeed, shear thinning behavior in non-Newtonian fluids may be conceptualized as thixotropic fluids with very rapid kinetics (Scott-Blair, 1951; Barnes, 1997).

Following the “eagle and rat” metaphor (Thompson, 2020) for homogeneous flows with time-varying shear rate, the parameter $k$ defines $\lambda_{eq}(\dot{\gamma})$ (position of the “rat”), while the parameter $m$ governs the rate at which $\lambda\rightarrow\lambda_{\text{eq}}$, analogous to the speed at which the “eagle” approaches the “rat”. However, for heterogeneous flows such as fully developed turbulence, this picture is more complicated as the shear rate varies spatially as well as temporally, in addition to the feedback mechanism that exists between the microstructure, rheology and turbulence.

Similar to the thixotropic model, for simplicity we also consider the simplest non-trivial rheological model given a purely viscous (non-elastic) Moore fluid (Moore, 1959)

where $\eta_{0}$ and $\eta_{\infty}$ represent the limiting viscosities for the structured $(\lambda=1)$ and unstructured $(\lambda=0)$ material. The coupling between the governing equations (1), (6), (8) directly quantify the feedback loop between the turbulence structure, microstructural evolution and the bulk rheology.

The governing equations (1), (6), (8) are non-dimensionalised with respect to the characteristic lengthscale $D$, advective timescale $\tau_{v}\equiv D/U_{b}$ and the wall viscosity

is chosen as a viscosity scale (Rudman et al., 2004), where subscript $w$ denotes values at the pipes walls. Henceforth all variables are non-dimensional and the same notation is used for simplicity. The resultant non-dimensional governing equations are then

In (10)-(12), the generalised Reynolds number

characterizes the ratio of inertial to viscous forces, the viscosity ratio $\eta_{r}\equiv\eta_{0}/\eta_{\infty}=2$ characterizes the ratio between the fully structured and fully unstructured viscosity, and the non-dimensional body force magnitude $F\equiv|d\langle{p}\rangle/dz|D/\rho U_{b}^{2}=1.8\times 10^{-2}$ characterizes the ratio of inertial to body forces. With respect to evolution of $\lambda$, the diffusion timescale is defined as $\tau_{\alpha}\equiv D^{2}/\alpha$ and the thixotropic reaction timescale is $\tau_{r}\equiv 1/m$, which characterize the Péclet and Damkhöler numbers as

which respectively characterise the timescales of advection and thixotropy to the diffusive timescale. A related dimensionless parameter known as the thixoviscous number (Ewoldt & McKinley, 2017)

which characterises the timescale of thixotropic kinetics relative to that of advection. Here $\Lambda\gg 1$ represents thixotropic fluids that evolve on timescales much faster than the flow, whereas $\Lambda\ll 1$ corresponds to thixotropic fluids that evolve over several advection times. The non-dimensional equilibrium parameter $K\equiv k\dot{\gamma}_{\text{nom}}$ in (11) characterises the equilibrium value of $\lambda$ under the nominal shear rate $\dot{\gamma}_{\text{nom}}$ as (7). This parameter although is defined as $K\equiv k(U_{b}/D)$, however, we have characterised it as $K\equiv k(8U_{b}/D)=0.4$ to ensure that the breakdown and the rebuild occur at roughly similar rates. Representative values of the key dimensionless variables used in the simulations are summarised in Table 1. The generalised Reynolds number $Re_{G}$ is in the moderately turbulent regime, and varies by roughly a factor of $\eta_{r}$ due to change in viscosity with fixed $F$. The thixoviscous number $\Lambda$ is varied from slow $(\Lambda\ll 1)$ to fast $(\Lambda\gg 1)$ kinetics. Although the Péclet number $Pe\gg 1$ indicates the transport of $\lambda$ is always advection dominated, the broad range of $\Lambda$ means that the Damkhóler number $Da$ is order unity for slow kinetics, meaning diffusion occurs on a similar timescale to thixotropy in this case.

The DNS method utilized in this study is based on an generalised Newtonian (GN) extension nnewt (Rudman & Blackburn, 2006) of the spectral element code Semtex (Blackburn et al., 2019), which has been previously validated for the turbulent pipe flow simulations of GN fluids (Rudman & Blackburn, 2006; Singh et al., 2016; Yousuf et al., 2024). This code is further extended to simulate turbulent flow of thixotropic fluids, as is described below. The DNS method is well suited for simulation of thixotropic turbulence as it resolves the turbulent flow across all spatio-temporal scales, eliminating the need for sub-grid scale turbulent closures, and facilitating direct solution of the transport equation (11) for the evolution of $\lambda$.

The pipe cross-section is divided into discrete two-dimensional quadrilateral elements representing standard tensor-product Lagrange interpolants with Gauss-Lobatto-Legendre collocation points. The axial direction is effectively managed through the application of Fourier expansion, which naturally imposes periodic boundary conditions (Temperton, 1992). The temporal resolution is handled by the second-order time integration method (Karniadakis et al., 1991) to maintain numerical stability. The code strictly monitors the Courant-Friedrichs-Lewy (CFL) criterion as well as the average divergence of the numerical solution for diagnostics.

To simulate thixotropy, the advective non-linear terms in (10) and (11) are discretized explicitly, while the diffusive terms are handled implicitly to ensure numerical stability and convergence, following the methodology outlined by Rudman & Blackburn (2006). For the breakdown and rebuild terms in (11), a fully implicit novel numerical scheme has been integrated into the code to enhance numerical robustness, see Appendix A for more details.

In the present study, the thixoviscous number $\Lambda$ is systematically varied in DNS computations from fast ($\Lambda=10^{2}$) to slow kinetics ($\Lambda=10^{-2}$). Simulations at higher $\Lambda$ values were found to be infeasible, as sharp gradients in the $\lambda$ field led to numerical instabilities, even with the implicit solver detailed in Appendix A for the thixotropic kinetics. Similarly, computations for very small values of $\Lambda$ were also found infeasible, as the timescale to reach thixotropic stationarity becomes arbitrarily greater than the eddy turnover time, thus requiring extensive computational overhead. However, as shall be shown in §§4.2 - 4.3, these computational limits do not restrict understanding of thixotropic turbulence as the dynamics in the limits $\Lambda\rightarrow\infty$ and $\Lambda\rightarrow 0$ can be easily understood and the behaviour at endpoints of this finite range $\Lambda\in[10^{-2},10^{2}]$ is close to that given by these limits. To provide reference cases, we also compute two Newtonian turbulent flow simulations corresponding to fluids with fully structured ($\lambda=1$, $\eta=\eta_{0}$) and unstructured ($\lambda=0$, $\eta=\eta_{\infty}$) microstructures.

The mesh design for these DNS cases adhere to the established guidelines for Newtonian (Piomelli et al., 1997) and non-Newtonian turbulence (Rudman et al., 2004; Rudman & Blackburn, 2006). The pipe cross-sectional mesh has 300 spectral elements with ninth-order tensor-product shape functions, and 384 axial data planes, corresponding to 11.52 million nodal points. Note that the structural parameter diffusivity $\alpha$ was chosen to yield a moderate Péclet number ($\Pen=10^{3}$) to ensure numerical stability, but this value is much smaller than characteristic values ($\Pen\sim 10^{12}$) based on self-diffusivity of e.g. colloidal suspensions (Morris & Brady, 1996). A summary of the computational parameters is given in Table 2. The DNS code monitors the Courant-Friedrichs-Lewy (CFL) criterion which was kept at $~{}10^{-1}$ to maintain numerical stability and convergence. Each case was run on the Gadi resource on the National Computational Infrastructure (NCI), Australia, which is a HPC cluster comprised of 8 x 24-core Intel Xeon Platinum 8274 (Cascade Lake) with 3.2 GHz CPUs and 192 GB RAM per node.

The DNS computations for the thixotropic and the structured ($\lambda=1$) Newtonian cases were initialised by introducing isotropic Gaussian-distributed white noise (with variance 0.01) to the velocity field of a fully developed laminar pipe flow, and the $\lambda$ field is set to unity. For the unstructured ($\lambda=0$) Newtonian case, the $\lambda$ field is set to zero during initialisation. Under fully developed conditions the thixotropic turbulent pipe flow is symmetric and thus stationary in time, axial, and azimuthal directions but non-stationary in the radial coordinate. Each DNS computation was run until the velocity and $\lambda$ fields both achieved statistical stationarity (typically 30 wash-through times), after which Eulerian turbulent flow statistics were accumulated and analyzed for another 30 wash-through times. In addition, Lagrangian data (velocity, velocity gradient and $\lambda$) was also collected for each run for $10^{4}$ randomly placed tracer particles over a period of around 8 wash-through times.

In this section we examine the characteristics of fully developed thixotropic turbulent pipe flow from an Eulerian perspective. Figures 1 and 2 respectively show representative cross-sectional plots of axial velocity $u_{z}$ and structural parameter $\lambda$ fields for for the two Newtonian and three thixotropic cases. At larger values of $\Lambda$, sharper and more fine-scale structures in the $\lambda$ field (and thus viscosity) are observed due to strong coupling between $\lambda$ and the instantaneous shear field, and from (7), in the limit $\Lambda\rightarrow\infty$ the structural parameter is solely dependent upon the local shear rate $\dot{\gamma}$, as demonstrated in Figure 1e. In this fast kinetics regime, the flow behavior effectively reduces to that of a shear-thinning rheology and the turbulence structure exhibits characteristics similar to that of shear-thinning turbulence (Rudman et al., 2004; Singh et al., 2017).

For smaller values of $\Lambda$, the structural parameter evolves more slowly along pathlines and the relevant shear rate history that governs $\lambda$ is longer. This leads to a reduction in fine-scale structures observed in the $\lambda$ field due to an effective averaging process backwards along pathlines. Although this results in more diffuse $\lambda$ distributions, this behavior persists in the limit of vanishing diffusivity ($\Pen\rightarrow\infty)$ due to the averaging process. For these computations, diffusion of $\lambda$ plays a minor role due to the moderate Péclet number ($Pe=10^{3}$) required for numerical stability. As the Damkhöler number scales with $\Lambda$ as $Da=\Lambda\Pen$, and $Da=10$ for $\Lambda=10^{-2}$, diffusion acts on a similar timescale as the thixotropic kinetics in this case, leading to the discrepancies observed in Figure 6c. Note that for complex fluids $Pe\gg 10^{3}$, and so for these flows smoothing of the $\lambda$ field is solely due to averaging of the Lagrangian shear rate. In the limit of $\Lambda\rightarrow 0$, the $\lambda$ field is almost uniform, as demonstrated in Figure 2, meaning that the viscosity is likewise and so the flow resembles that of Newtonian turbulence (Toonder & Nieuwstadt, 1997; Khoury et al., 2013).

Mean field and second order statistics also provide insights into the impact of thixotropy on turbulent pipe flow. The mean radial viscosity and $\lambda$ profiles for the three thixotropic cases are illustrated in Figure 3, and the Newtonian reference cases provide upper and lower bounds for these quantities.

For all thixotropic cases, shear degradation in the viscous sub-layer ($y^{+}\leq 5$) is more significant than in the pipe core due to increased shear near the pipe walls. This is more pronounced for the faster thixotropic kinetics, as the structural parameter $\lambda$ exhibits a stronger correlation with local shear rate $\dot{\gamma}$, leading to a monotone decreasing radial viscosity profile similar to that of turbulent pipe flow of shear thinning fluids (Rudman et al., 2004; Singh et al., 2017). Similar to the more diffusive profiles shown in Figure 6, this effect is less pronounced for the slower thixotropic kinetics, where the radial profiles of $\lambda$ and $\eta$ are significantly flatter, and the viscosity profile approaches a constant in the limit $\Lambda\rightarrow 0$.

The influence of viscosity on turbulence intensity is illustrated in Figure 4. For the Newtonian reference cases, the profiles of all the turbulence intensities (except $u^{\prime}_{zz}$) are higher for the unstructured case $\lambda=1$, with peaks for $u^{\prime}_{rr}$ and $u^{\prime}_{rz}$ moving further away from the walls, which is due to the increased $Re_{G}$ for the unstructured case. These observations are consistent with well-established trends of Newtonian turbulence  (Toonder & Nieuwstadt, 1997; Khoury et al., 2013). For the thixotropic cases, the azimuthal turbulence intensity $u^{\prime}_{tt}$ and the $u^{\prime}_{rz}$ Reynolds stress transition monotonically with increasing $\Lambda$ from the fully structured ($\lambda=1$) to the unstructured ($\lambda=0$) case, in terms of both turbulence intensity and peak shift.

The profiles of the axial turbulence intensity $u^{\prime}_{zz}$ transition monotonically in terms of turbulence intensity, whereas, the profiles of the radial turbulence intensity $u^{\prime}_{rr}$ transition monotonically with in terms of peak shift.

The mean axial velocity profiles shown in Figure 5a demonstrate that all flow profiles roughly follow a linear relationship $U_{z}^{+}=y^{+}$ in the viscous sub-layer (Pope, 2001). Figure 5b provides a clearer view of these profiles, with the deviation monotonically increases with increasing $\Lambda$. The peak deviation occurs in the buffer layer ($30\leq y^{+}\leq 2000$) for all the DNS cases due to increased turbulent intensities in that region (see Figure 4). For the thixotropic cases, the area integral of the deviation over the pipe cross-section gives us the excess bulk flow rate (or drag reduction), thus drag reduction increases with $\Lambda$, which is consistent with previous studies (Escudier & Presti, 1996; Pereira & Pinho, 2001).

The results in §3 suggest that the shear history along pathlines plays a critical role in organising thixotropic turbulence. As the ADRE (11) for $\lambda$ is advection-dominated ($\Pen\gg 1$), the Lagrangian frame is well-suited to study the feedback mechanisms that govern thixotropic turbulence, as this naturally encodes the shear history experienced by fluid elements. In this section we develop a Lagrangian model of thixotropy that provides insights into the feedback mechanisms that govern thixotropic turbulence and use these insights to develop analogue time-independent rheological models.

We first consider evolution of the structural parameter $\lambda$ in the Lagrangian frame $\mathbf{X}$ as $\lambda(t;\mathbf{X},t_{0})=\lambda(\mathbf{x}(t;\mathbf{X},t_{0}),t)$, where $\mathbf{x}(t;\mathbf{X},t_{0})$ is the solution of the advection equation for a fluid particle initially at position $\mathbf{X}$ at time $t=t_{0}$:

In the limit of vanishing diffusivity ($Pe\rightarrow\infty$), the ADRE (11) can be transformed to the Lagrangian frame as

Solution of (17) shows that the structural parameter evolves with respect to Lagrangian shear history as

where $g(t,\mathbf{X},t_{0})\equiv 1+K\dot{\gamma}(t,\mathbf{X},t_{0})$. As we are interested in the long-time behavior of $\lambda(t)$, given by $t_{0}\rightarrow-\infty$ (or more accurately $t-t_{0}\gg 1/\Lambda$), the initial condition $\lambda_{0}$ has no impact on $\lambda$ and (18) may be recast as

Here $G$ in (19) represents a fading memory kernel $\hat{\mathcal{F}}$

that encodes the shear history from $s=0$ to $s\rightarrow-\infty$. Hence the most recent shear rates have the greatest impact, and the persistence of memory (analogous to relaxation time for viscoelastic fluids) is controlled by the thixoviscous number.

Equations (19), (20) also shows that advective transport is an important mechanism in non-stationary thixotropic turbulent flow. From a stochastic perspective, the Lagrangian shear rate history $\dot{\gamma}(t-s;\mathbf{X},t_{0})$ can be considered as a random process that is non-stationary in space and so is also dependent upon the history $s$ of Eulerian positions $\mathbf{x}(t-s;\mathbf{X},t_{0})$ along pathlines. For fully-developed turbulent pipe flow the shear rate $\dot{\gamma}$ is non-stationary in the radial direction $r$, hence the functional (20) simplifies to $G(r(t),t)=\hat{\mathcal{F}}[\dot{\gamma}(t-s;r(t-s))|_{s=0}^{\infty}]$, and the evolution equation (19) needs to be supplemented by an advective transport model for evolution of radial position $r(t)$. This stochastic approach will be considered further in §§5.2.

Conversely, for statistically stationary turbulent flows (20) simplifies to $G(t)=\hat{\mathcal{F}}[\dot{\gamma}(t-s)|_{s=0}^{\infty}$. If the Lagrangian shear rate follows a simple autocorrelation process with decorrelation time $\tau_{\dot{\gamma}}$

then $\dot{\gamma}(t-s)$ may be modeled as a Markovian random process, leading to a fairly simple random walk model for the evolution of $\lambda$ via (19).

To test the Lagrangian assumption, the evolution of $\lambda$ along different pathlines is compared in Figure 6 between DNS results and computation of $\lambda$ via in (17). For the fast thixotropic kinetics ($\Lambda=10^{2}$) with a high Damköhler number ($Da=10^{5}$), the DNS results closely match the Lagrangian solution (17), with a relative $L_{2}$ error of 0.18% over $10^{4}$ trajectories. For intermediate kinetics ($\Lambda=1$), diffusion is significant although thixotropic kinetics still dominate ($Da=10^{3}$), and the Lagrangian solution (17) has relative $L_{2}$ error of 1.79%. For the slow kinetics ($\Lambda=10^{-2}$), the timescales of diffusion approach those of the thixotropic kinetics ($Da=10$), leading to noticeable discrepancies in pathlines and relative $L_{2}$ error of 4.03%. The discrepancies in Figure 6c are due to to the moderate Péclet number ($\Pen=10^{3}$) used for numerical stability, corresponding to $Da\sim 10$ for $\Lambda=10^{-2}$. As previously discussed, both $\Pen$ and $Da$ are much higher for complex fluids, hence the Lagrangian framework is broadly applicable to thixotropic flows. Although the Lagrangian assumption does not strictly hold in this study for $\Lambda\ll 1$, we shall show in the following that this does not significantly alter the overall dynamics of the flow, and and accurate analytic closures for the limiting cases of $\Lambda$ can be developed.

In the fast kinetic regime, $\Lambda\rightarrow\infty$, the thixotropic timescale $\tau_{r}$ is much shorter than the eddy turnover time $\tau_{v}$ that governs the Lagrangian shear rate decorrelation time, hence only the very recent shear history dictates $\lambda$. Following (19), we derive an expression for the structural parameter $\lambda$ in the limit $\Lambda\rightarrow\infty$ by first performing a series expansion of $\dot{\gamma}(t-s;\mathbf{X},t_{0})$ about $t$ as

where from the definition of $\tau_{\dot{\gamma}}$, the coefficients $\dot{\gamma}_{n}(t;\mathbf{X},t_{0})\equiv\dot{\gamma}(t;\mathbf{X},t_{0})^{(n)}\tau_{\dot{\gamma}}^{n}/n!\sim 1$, and $x^{(n)}(t)$ denotes the $n$-th derivative of $x$ with respect to $t$. Introducing the rescaled thixotropic time $t_{r}=t\Lambda$, the structural parameter evolves according to

where the history function $G(t;\mathbf{X},t_{0})$ in the limit $\Lambda\rightarrow\infty$ is then

Application of L’Hopital’s rule to (23) yields

Hence the structural parameter in Eulerian space for $\Lambda\rightarrow\infty$ is given by the equilibrium value

For $\Lambda\rightarrow\infty$ the structural parameter instantly equilibrates with respect to the local shear rate, hence the effective viscosity of the thixotropic psuedo-Newtonian fluid is a time-independent shear-thinning generalised Newtonian fluid (Scott-Blair, 1951; Barnes, 1997) that corresponds to a Cross fluid with unit index

Figure 7 verifies (25) via comparison against DNS computations for $\Lambda=10^{2}$. Figure 7a shows that the fast kinetics model (25) for evolution of $\lambda$ along sample trajectories matches DNS results very well, with a relative $L_{2}$ error of 0.1% across $10^{4}$ Lagrangian trajectories. Figure 7b also shows that the conditional p.d.f. $p_{\lambda|r}(\lambda|r)$ of $\lambda$ at various radial positions from (25) also matches the DNS results very well, verifying both the Lagrangian assumption and the closure model (25).

To test the closure (25) and the corresponding effective viscosity model (27), DNS computations using (27) are compared with those using the thixotropic kinetics with $\Lambda=10^{2}$ in Figures 3 - 5. Computations based on the effective viscosity model (27) show excellent agreement with the full thixotropic model, with relative errors of under 0.1% for the mean structural parameter $\lambda$, viscosity $\eta$ and axial velocity $U_{z}$. The errors in turbulent statistics ($u^{\prime}_{rr},u^{\prime}_{tt},u^{\prime}_{zz},u^{\prime}_{rz}$) are also all less than 0.5%. These results confirm that, as expected, for $\Lambda\gg 1$ thixotropic fluids exhibits turbulent fluid dynamics akin to that of shear-thinning fluids. Note that this result extends to a broad range of rheologies and thixotropic models described by (3), (5) respectively.

To analyse thixotropic rheology in the limit of vanishingly slow thixotropic kinetics, $\Lambda\rightarrow 0$, we decompose the memory kernel into an ensemble average $\langle\dot{\gamma}\rangle$ that, due to ergodicity, corresponds to the Lagrangian average

and a fluctuating component $\dot{\gamma}^{\prime}$ as

where the scaling $t/\tau_{\dot{\gamma}}$ ensures $\partial\dot{\gamma}^{\prime}/\partial t\sim 1$. Using this decomposition, the memory kernel (19) is then

and integration by parts yields

where $t_{2}=t\tau_{\dot{\gamma}}$. In the limits $t_{0}\rightarrow-\infty$, $\Lambda\rightarrow 0$, the structural parameter given by (19) then simplifies to

Hence for $\Lambda\rightarrow 0$, the thixotropic kinetics are so slow compared to the fluctuation rate of the shear rate that the structural parameter effectively samples the ensemble of shear rates during convergence to equilibrium, and so the structural parameter is steady and is governed by the ensemble averaged shear rate $\langle\dot{\gamma}\rangle$. Hence the structural parameter in Eulerian space for $\Lambda\rightarrow 0$ is given by the equilibrium value at the mean shear rate as:

and so the thixotropic viscosity for $\Lambda\rightarrow 0$ simplifies to the Newtonian viscosity

Figure 8 verifies the slow kinetics analytic closure (32) by comparison with DNS computations for $\Lambda=10^{-2}$. Despite that the Lagrangian framework does not strictly apply for this case, Figure 8a shows that time-series of $\lambda$ values along sample trajectories computed via DNS fluctuate close to the equilibrium value given by (32), with a relative $L_{2}$ error of 1.08% across $10^{4}$ Lagrangian trajectories. Similarly, Figure 8b shows that the conditional p.d.f. $p_{\lambda|r}(\lambda|r)$ from the DNS computations at different radial positions $r$ all have small variance around the equilibrium value given by (32). These results show that despite breakdown of the Lagrangian approximation (17) for $\Lambda\ll 1$ due to significant diffusion of $\lambda$ (i.e. $Da=10$), the slow kinetics closure (32) is still accurate as the sampling of shear rates that governs convergence of $\dot{\gamma}(t;\mathbf{X},t_{0})\rightarrow\langle\dot{\gamma}\rangle$ in (32) is not affected by the diffusivity of $\lambda$.

To test the closure (32) and the corresponding effective viscosity model (34), DNS computations using (34) are compared with those using the thixotropic kinetics with $\Lambda=10^{-2}$ in Figures 3 - 5. Computations based on the effective viscosity model (34) show good agreement with the full thixotropic model, with relative errors of 2.24% for the mean structural parameter $\lambda$, 0.95% for the mean viscosity $\eta$ and 0.29% for the axial velocity $U_{z}$. The slightly higher errors than for the fast kinetics closure may be partially due to breakdown of the Lagrangian approximation in the slow kinetics case. The errors in turbulent statistics ($u^{\prime}_{rr},u^{\prime}_{tt},u^{\prime}_{zz},u^{\prime}_{rz}$) are also all less than 0.7%. These results confirm that, as expected, for $\Lambda\ll 1$ thixotropic fluids exhibits turbulent fluid dynamics akin to that of a generalised Newtonian fluids. Note that this result extends to a broad range of rheologies and thixotropic models described by (3), (5) respectively.

As discussed in §2, fully-developed turbulent pipe flow is radially non-stationary, hence the random process governing the Lagrangian shear rate history $\dot{\gamma}(t-s;\mathbf{X},t_{0})$ is also radially non-stationary. Specifically, the microstructural evolution equation (19) may be expressed as

where the Lagrangian shear rate $\dot{\gamma}(t,r(t))$ is considered as a random variable that is non-stationary with respect to $r$ with conditional p.d.f. $p_{\dot{\gamma}|r}(\dot{\gamma}|r)$. Hence we require coupled stochastic models for both the radial position $r(t)$ and for the Lagrangian shear rate $\dot{\gamma}(t,r(t))$ which is conditional on $r(t)$. From (35), the dependence of the structural parameter $\lambda$ at time $t$ and position $r(t)$ on the Lagrangian shear history can be encoded as

and so the evolution of $\lambda$ follows a non-Markov process, while the evolution of $\dot{\gamma}(t,r(t))$ and $r(t)$ may be Markovian in an appropriate frame. Following the GN models developed in §4, we propose a simple viscosity model for arbitrary $\Lambda$ where, based on the radial non-stationarity of the flow, the local viscosity at any radial position $r$ is based on the conditionally averaged structural parameter $\langle\lambda|r\rangle$ as

where

and $p_{\lambda|r}(\lambda|r)$ is the conditional probability of $\lambda$ occurring at radial position $r$. This viscosity model (37) effectively assumes that the variation in $\lambda$ about its conditionally average has minimal impact on the flow, which shall be tested in due course. Discretizing the $s$-domain as $s_{n}=n\Delta s$, with $\Delta s\ll\tau_{\dot{\gamma}}$, the functional (36) can be expressed as

where $\dot{\boldsymbol{\gamma}}$ is the vector $\dot{\boldsymbol{\gamma}}\equiv(\dot{\gamma}_{0},\dot{\gamma}_{1},\dots,\dot{\gamma}_{q})$ of Lagrangian shear rates $\dot{\gamma}_{n}\equiv\dot{\gamma}(t-n\Delta s;\mathbf{X},t_{0})$. As $G(t-s)\sim\exp(-\Lambda s)$, if $q\gg 1/(\Delta s\Lambda)$ then the shear history for $s>q\Delta s$ has no impact on $\lambda(t,r(t))$. Similarly, the discrete Lagrangian radial history can be encoded as $\mathbf{r}=\{r_{0},r_{1},\dots,r_{q}\}$, where $r_{n}\equiv r(t-n\Delta s)$ for $n=0:q$.

From (36), the conditional average (38) can be conceptualised as a path integral (Kleinert, 2006; Wio, 2013) of $\mathcal{F}$ over all the shear $\dot{\boldsymbol{\gamma}}$ and radial $\mathbf{r}$ histories backwards in time from a particle at position $r(t)$ at current time $t$ to $r(t-q\Delta s)$ at time $t-q\Delta s$ as

where $\mathcal{P}[\dot{\boldsymbol{\gamma}},\mathbf{r}|r_{0}=r]$ is the conditional probability of the shear history $\dot{\boldsymbol{\gamma}}$ and radial history $\mathbf{r}$ given current radial position $r$. This general formalism describes evolution of the structural parameter over a broad range of thixotropic and rheological models.

Via Bayes’ theorem, the conditional probability $\mathcal{P}[\dot{\boldsymbol{\gamma}},\mathbf{r}|r_{0}=r]$ is related to the joint probability $\mathcal{P}[\dot{\boldsymbol{\gamma}},\mathbf{r}]$ of shear rate history $\dot{\boldsymbol{\gamma}}$ and radial history as

where $\mathcal{P}[\mathbf{r}|r_{0}=r]$ is the conditional probability of $\mathbf{r}$ given $r_{0}=r$. Under the assumption that the shear rate and radial position evolve in a Markovian manner in time, the joint probability $\mathcal{P}[\dot{\boldsymbol{\gamma}},\mathbf{r}]$ can then be expanded as

where the backwards propagator $P_{\Delta s}(\dot{\gamma}_{n+1},r_{n+1}|\dot{\gamma}_{n},r_{n})$ is the conditional probability of a material element being at position $r_{n+1}$ with shear rate $\dot{\gamma}_{n+1}$ at time $t-(n+1)\Delta s$ given it is at position $r_{n}$ with shear rate $\dot{\gamma}_{n}$ at time $t-n\Delta s$. Using the Chapman-Kolmogorov equation for Markov processes