Thixotropic pumping in a cylindrical pipe

We consider the flow of a thixotropic fluid in a uniform cylindrical pipe, driven by an oscillating pressure gradient or a body force. For a variety of rheological models, solutions can be obtained by integrating ordinary rather than partial differential equations: we illustrate this approach for the thixo-viscoplastic Houska model and the thixo-viscous simplified Moore-Mewis-Wagner model. We present asymptotic results in the limits of small and large Deborah numbers, and numerical results for intermediate Deborah numbers. Under asymmetrical "sawtooth" forcing, thixotropy leads to the net transport of fluid along the pipe, even when there would be no net transport of the corresponding generalised-Newtonian fluid. We propose the name "thixotropic pumping" for this novel transport mechanism.


I. INTRODUCTION
We consider the flow of a thixotropic fluid in a uniform cylindrical pipe, driven by an oscillating pressure gradient or body force. (An everyday analog occurs when a bottle of tomato ketchup is gently shaken to break down the structure before it is poured.) If the oscillation is asymmetric, then the hysteretic response of the fluid can lead to net transport along the pipe, even though a nonthixotropic fluid would experience no net transport. This distinctive mechanism, which does not appear to have been investigated previously, may be referred to as thixotropic pumping.
Such paradigm problems can become particularly rich if processes with an externally imposed timescale interact with the internal thixotropic dynamics. A particularly relevant example is the thixotropic version of the Stokes problem in which rectilinear flow is driven by an oscillating wall. McArdle et al. [17] obtained asymptotic descriptions of this flow in the limits in which the thixotropic response was much faster or much slower than the oscillation, and presented numerical results in the intermediate regime in which the dynamics are most complicated. We will take a similar approach in the present work, although the dynamics are rather different: the crucial interaction in the Stokes problem is between thixotropy and inertia rather than between thixotropy and an imposed pressure gradient.
In this study we focus on the thixoviscoplastic Houška rheology. However, the approach can be generalized to other rheologies. We also briefly present illustrative results for the thixoviscous simplified Moore-Mewis-Wagner (sMMW) rheology, which allows us to explore the effect of shear-thickening (antithixotropic) as well as shear-thinning (thixotropic) behavior.

A. Governing equations
We consider axisymmetric rectilinear flow in a cylindrical pipe of radiusR, driven by a periodic pressure gradient with angular frequencyω and maximum valueĜ 0 . We present the governing equations in nondimensional form, where we have nondimensionalized using the length scaleR, the pressure-gradient scaleĜ 0 , the inverse timescaleω, a viscosity scale given by the minimum viscosityμ 0 of the fluid, and a velocity scaleĜ 0R 2 /μ 0 (see Ref. [18, chapter 7] for details). The Cauchy stress equation is where τ is the shear stress and G(t ) = −∂ p/∂z is the pressure gradient. We assume that G(t ) is periodic with period 2π and amplitude 1; unless otherwise stated, G(t ) > 0 for 0 < t < π (the "forward phase") and G(t ) < 0 for π < t < 2π (the "reverse phase"). We will consider two cases in particular: sinusoidal forcing, and "sawtooth" forcing, where and where the first peak of G(t ) occurs at t = T 0 . For 0 T 0 < π/2 the pressure gradient increases more rapidly than it decreases (we will refer to this as a short rising leg/long falling leg); for π/2 < T 0 π the pressure gradient decreases more rapidly than it increases (a long rising leg/short falling leg). We may integrate Eq. (1) subject to the symmetry condition τ = 0 on r = 0 to obtain We consider ideal thixotropic fluids [19], so the constitutive equation defines τ in terms of the shear rateγ (r, t ) = |∂w/∂r|, where w(r, t ) is the axial velocity, and a scalar "structure parameter" λ(r, t ) [20] which describes the local state of the fluid microstructure, The structure parameter evolves according to an equation of the form Here we have introduced the Deborah number, the ratio of the forcing frequencyω to the structure response ratef 0 . The function f (γ , λ) is a dimensionless net buildup/breakdown rate which depends on both the shear rateγ and the structure parameter λ. If the constitutive relation takes a suitable form, then we can use Eqs. (5) and (6) to write |∂w/∂r| in terms of r, t, and λ; Eq. (7) can then be integrated in t to obtain λ for any r. This is the approach taken in this study.
The velocity is zero at the pipe wall r = 1, No boundary condition is required for λ. Initial conditions will be specified when required, but our focus is on situations in which the solutions are periodic, We may expect this to be the long-term asymptotic behavior of the system under periodic forcing.

B. Integrated measures of transport
It is helpful to define three quantities that describe the transport of fluid. The instantaneous flux of fluid along the pipe is defined by Integrating Q(t ) over time, we obtain the total volume of fluid transported during the forward phase, and the total volume of fluid transported during a full period, Under symmetric forcing, or in the absence of thixotropic effects, we expect V 2π to be zero. If V 2π = 0, then there is net transport in either the positive axial direction (V 2π > 0) or the negative axial direction (V 2π < 0).

C. Asymptotic limits
As in the thixotropic Stokes problem [17], we can make analytical progress in two limits of the Deborah number.

Fast adjustment (D = 0)
In the limit D = 0, the left-hand side of Eq. (7) is zero, so the structure parameter adjusts instantaneously to the equilibrium value λ = λ eq (r, t ) that satisfies f ∂w ∂r , λ eq = 0.
The fluid thus has a generalized-Newtonian rheology in which the shear stress τ is an instantaneous function of the shear rate |∂w/∂r|, defined implicitly by Eq. (14) together with Eq. (6).

Slow adjustment (D → ∞)
In the limit D → ∞, the structure parameter becomes independent of time. However, for periodicity to be maintained, we require the buildup and breakdown terms in Eq. (7) to balance when integrated over a period. This can be seen by expanding λ for large D as which leads to Equation (16) implies that at leading order the structure parameter depends only on position, λ 0 = λ av (r). Integrating Eq. (17) over a period and applying periodicity then yields and this condition implicitly defines λ av (r).

III. HOUŠKA RHEOLOGY
The Houška rheology [21,22] is a tractable and popular model [1,3,5,7,8] which combines thixotropy with a yield stress. The structure parameter evolves according to Eq. (7), with f given by a special case of Mewis and Wagner's [20] general buildup/breakdown function, Shear drives breakdown, so in an unsheared region the structure can only evolve by buildup; thus in a permanently unsheared region the only long-term equilibrium state of the fluid is fully structured, λ = 1. In general, λ lies within the range 0 λ 1. The constitutive relation is based on that of a Bingham fluid [23,Sec. 4.2]. In simple shear flow it becomes The parameter κ > 0 describes how rapidly buildup takes place relative to breakdown; the parameters τ y0 0, τ y1 0 and η H1 0 describe, respectively, the minimum yield stress and the rates of change of the yield stress and of the Bingham viscosity with the structure parameter. Note that the dimensionless minimum viscosity is 1. Under oscillatory forcing, it is necessary to keep track of up to three regions of the flow. Close to the center of the pipe is a region in which the shear stress given by Eq. (5) never exceeds the maximum value of the yield stress, τ y0 + τ y1 . In this core region, 0 r < r c , the fluid is fully structured and unsheared throughout the oscillation, λ = 1. Since max t |G(t )| = 1, the boundary of the core region is given via Eq. (5) by r c = 2(τ y0 + τ y1 ).
At any instant, there may also be a wider plug region, 0 r < r y (t ), in which the fluid is (instantaneously) unyielded, τ (r, t ) < τ y0 + λ(r, t )τ y1 . The outer boundary of the plug region, r y (t ), moves in and out as the driving force increases and decreases; when r y (t ) = 1 the entire flow becomes unyielded. Because the fluid in the annulus r c < r < r y (t ) is not in general fully structured, the location of r y (t ) depends on the solution for λ(r, t ) (except in the simplest case τ y1 = 0, in which the yield stress is independent of the structure parameter and so r y = r c throughout the oscillation).
Finally, at any instant there may also be a yielded region extending from the outer boundary of the plug as far as the wall, r y (t ) < r < 1. During phases in which the pressure gradient is small, this region ceases to exist because the fluid becomes unyielded everywhere across the pipe.

Fast adjustment (D = 0)
In the limit of fast adjustment, Eq. (19) yields Combining Eq. (21) with Eqs. (20) and (5), we find that the structure parameter is given by the positive solution of the quadratic equation when this lies between 0 and 1, and by λ eq = 1 otherwise. (The explicit solution for λ eq , which is presented by Ref. [3], is not informative for our purposes.)

Slow adjustment (D → ∞)
In the limit of slow adjustment, the structure parameter is determined by Eq. (18).
and thus λ av To evaluate I (r, λ) we require an expression for the shear during phases when the fluid is yielded. Combining Eqs. (5) and (20), we obtain Under sinusoidal forcing, Eq. (2), it is then helpful to define and to write Equation (24) can now readily be solved numerically to obtain λ av (r). Under sawtooth forcing, Eq. (3), the calculation is simpler, and we obtain 013303-5

Numerical integration for 0 < D < ∞
To obtain numerical results we substitute Eq. (25) into Eq. (19) to obtain a time-evolution equation for λ, Equation (30) can then be integrated forward in time using a standard method. In our calculations we took the initial condition λ(r, 0) = λ av (r) in all cases to ensure rapid convergence to the asymptotic solution for large D. All results presented here had been integrated for least four full periods, by which time periodicity had been reached.
Once λ(r, t ) has been obtained for a sufficient number of values of r, profiles of w(r, t ) and the flux Q(t ) can be obtained by quadrature. The results presented below were obtained using the inbuilt Runge-Kutta routine in Maple 2019 [24], and quadrature over r was carried out on a minimum of N = 100 values of r; results were robust to the choice of N. For sawtooth forcing with small values of T 0 , significantly higher t-resolution (up to 2000 points in contrast to 200 for larger values of T 0 ) was required in order to capture the fluxes accurately.

B. Sinusoidal forcing: Numerical results
We first consider the behavior under sinusoidal forcing, Eq. (2). For very small values of D, the structure parameter λ tracks the instantaneous equilibrium value λ eq closely [ Fig. 1(a)]. There are phases during which the fluid is fully structured right across the pipe, and small lags between λ and λ eq are visible just after yielding and unyielding.

Structure parameter and velocity
As D increases, the lag also increases and the fully structured phases disappear [Figs. 1(b)-1(d)]. The lag is somewhat longer for smaller values of r (larger values of λ), where the shear rates and thus the breakdown rates are lower; it is also most pronounced around the minima of the shear rate (t = nπ/2 for n ∈ N), when not only is the breakdown rate smallest but λ eq is varying most strongly.
Thixotropic effects are clearest for intermediate values of D [ Fig. 1(c)]. As D increases further [Figs. 1(d) and 1(e)] the lag increases further, but at the same time the variation of λ over a period diminishes so the lag becomes less apparent. Figure 1(e) demonstrates that by D = 1 the large-D asymptotic solution λ ∼ λ av (r) provides a reasonable approximation everywhere across the pipe. Figure 2 illustrates the corresponding profiles of λ and w across the pipe, during a quarterperiod when the pressure gradient is increasing, for the limit of rapid adjustment D The core region is 0 r < r c = 0.4 in each case. In the rapidly adjusting limit [ Fig. 2(a)], this represents the minimal extent of the fully structured region; the outer boundary of the fully structured region moves in and out through the oscillation, sometimes reaching the pipe wall r = 1. In this limit, this outer boundary corresponds to the outer boundary r = r y of the plug region, as can be seen by comparing Fig. 2(a) with the corresponding velocity profiles in Fig. 2(d). (Note the black squares that mark the boundary of the plug region.) In contrast, for the other cases plotted the fluid is never fully structured outwith the core phases during which the yield surface reaches the pipe wall and the velocity drops to zero across the pipe. The period-averaged structure parameter λ av at the wall is somewhat higher than the minimum value taken at the wall by the instantaneous equilibrium structure parameter λ eq (t ) [cf. Figs. 2(a) and 2(c)]. As a result the maximum velocity decreases somewhat with increasing D [Figs. 2(d)-2(f)]. The other notable feature is that for intermediate values of D the lag is different at different radial positions, leading to the rather complicated variation seen in Fig. 2(b). In the limit D = 0 the maximum flux is highest but the fluid is slowest to yield at the wall, giving the taller and narrower dashed peaks in Figs. 3(a)-3(c); conversely, in the limit D → ∞ the maximum flux is lowest but the fluid is quickest to yield at the wall, giving the shorter and wider dashed peaks in Figs. 3(a)-3(c). Between these limits, the maximum flux Q max decreases monotonically with increasing D [Fig. 3(d)].

Transport of fluid
In each limit the velocity, and thus the flux, is in phase with the pressure gradient and so the maximum fluxes in each direction occur at t = π/2 and t = 3π/2. Between these limits the flux lags behind the pressure gradient, because as the pressure gradient increases the structure is  The broadening of the peaks of Q(t ) with increasing D partly compensates for the reduction in Q max , with the consequence that the total volume of fluid transported during the forward phase V π decreases less strongly with D than does Q max [cf. Figs. 3(d) and 3(e)].

C. Sawtooth forcing: Numerical results
We now consider the behavior of the fluid under sawtooth forcing, Eq. (3) . Figures 4(a)-4(c) illustrate the flux Q(t ) for various values of D and T 0 , all in cases in which the pressure gradient has a short rising leg. Figure 4(d) shows the total volume of fluid transported during the forward phase, V π , while Fig. 4(e) shows the total volume of fluid transported over a full period, V 2π .
In both the small-D limit [ Fig. 4(a)] and the large-D limit [ Fig. 4(c)], the velocity at any position r is an instantaneous function of the pressure gradient G(t ). Consequently, the flux is in phase with G(t ) throughout the oscillation, and there is symmetry between the forward and reverse phases. Thus, in each limit there is no net transport of fluid over a full period [ Fig. 4(e)].
For intermediate values of D, thixotropy plays a role. During phases when |G(t )| is increasing, the structure is breaking down; consequently |Q(t )| increases at first more slowly than it would without thixotropy. Conversely, during phases when |G(t )| is decreasing, the structure is rebuilding; consequently |Q(t )| decreases at first more slowly than it would without thixotropy. The consequence is the "shark's-tooth" shape of Q(t ) apparent in Fig. 4(b), with a concave-outward phase as |Q(t )| increases followed by a convex-outward phase as |Q(t )| decreases. When this shark's-tooth response is combined with an asymmetrical sawtooth forcing, the result is that during the longer falling leg Q(t ) is convex-outward for most of the time when it is positive, and Q(t ) is concave-outward when it is negative. [This is clearest for T 0 = 0 in Fig. 4(b).] Consequently, the total volume of fluid transported during the forward phase, V π , is enhanced for intermediate values of D [ Fig. 4(d)], and the total volume of fluid transported during the reverse phase is reduced. The result is that over a period thixotropic pumping occurs: there is net transport in the positive direction, V 2π > 0 [ Fig. 4(e)], and the total volume of fluid transported over a period may be a substantial fraction of V π . (This pattern is, of course, reversed when the asymmetry of the sawtooth is reversed, π/2 < T 0 < π.) Figure 5 demonstrates that the basic mechanism of thixotropic pumping persists when either the shear stress or the viscosity depends on the structure parameter. Figures 5(a) and 5(b) show the total volume of fluid transported during the forward phase, V π [ Fig. 5(a)] and the total volume of fluid transported over a full period, V 2π [ Fig. 5(b)], when the yield stress is independent of λ, i.e., τ y1 = 0. The overall pattern is very similar to that shown in Figs. 4(d) and 4(e), except that because the apparent viscosity of the flux is lower the total volume of fluid transported during the forward phase is higher [cf. Fig. 5(a) with Fig. 4(d)], and because the shear stress depends less strongly on λ the effect of thixotropic pumping is weaker [cf. Fig. 5(b) with Fig. 4(e)].  Fig. 5(d)], when the viscosity is independent of λ, i.e. η H1 = 0. The overall pattern is again similar to that shown in Figs. 4(d) and 4(e), although the effect of thixotropy on the net flux during the forward phase for strongly sawtooth forcing is now more conspicuous and the lower yield stress leads to higher values of V π [cf. Fig. 5(c) with Fig. 4(d)]. The effect of thixotropic pumping is again weaker when V 2π [ Fig. 5(d)] is considered as a fraction of V π .

IV. SIMPLIFIED MOORE-MEWIS-WAGNER RHEOLOGY
The simplified Moore-Mewis-Wagner (sMMW) rheology provides a convenient description of purely viscous thixotropic behavior, and has been used in several previous studies [4,17,25]. The model is defined by in Eq. (7), together with the constitutive relation (in simple shear flow) The sMMW rheology has the feature that in the limit D = 0 the fluid becomes a power-law fluid. In contrast to the Houška rheology, λ is no longer bounded above, 0 λ < ∞. By combining Eqs. (5), (31), and (32), we obtain a time-evolution equation for λ at each point r, A. Asymptotic limits of the sMMW model

Fast adjustment (D = 0)
In the limit of fast adjustment, the fluid has a power-law rheology, in which the structure parameter and the stress are given by λ eq = κ 1/b ∂w ∂r n−1 Equation (5) may be integrated to give the velocity profile w(r, t ) = sign(G) |G| and hence the flux

Slow adjustment (D → ∞)
In the limit of slow adjustment, the structure parameter is determined by Eq. (18), and the shear rate is given by We thus obtain In the case of sinusoidal forcing, Eq. (2), we can evaluate the integrals in Eq. (39) to obtain where B(x, y) is the standard Beta function.
The forms of the velocity profiles in the quickly and slowly adjusting limits, given by Eqs. (35) and (42), are identical but their amplitudes and time-dependences are not. (It is only in these limits that the velocity profiles take an identical form.) For thixotropic (shear-thinning) cases in which c < a and so n < 1, the peaks in the quickly adjusting flux, Eq. (36), are thinner but higher than those in the slowly adjusting flux, Eq. (43); for antithixotropic (shear-thickening) cases in which c > a and so n > 1, the peaks in the quickly adjusting flux, Eq. (36), are broader but lower than those in the slowly adjusting flux, Eq. (43).

B. Numerical results
To obtain numerical results we integrate Eq. (33) forward in time using a standard method, and obtain the velocity w(r, t ) and the flux Q(t ) by quadrature, as described in Sec. III A 3. . This occurs because of a competition between the increasing magnitude of the pressure gradient during this interval and the increasing viscosity as the structure builds up; the maximum value of Q(t ) occurs when the pressure gradient is large but the structure has not yet built up enough to slow the flow.
We now consider sawtooth forcing, Eq. (3), with a short rising leg, 0 T 0 π/2.  Fig. 7(c)] rather than increasing it, and the consequence is that the total volume of fluid transported during a period is negative [ Fig. 7(d)] rather than positive. The effect of thixotropy is maximized for somewhat higher values of D than it was for the shear-thinning fluid.

V. DISCUSSION AND CONCLUSIONS
We have investigated the pumping of a thixotropic fluid by an oscillating pressure gradient in a cylindrical pipe. This paradigm problem allows us to investigate the interaction between externally imposed forcing and the thixotropic dynamics, which may operate on very different timescales. At the same time, it is sufficiently simple that for a variety of rheological models we can obtain solutions by integrating ODEs rather than PDEs; it therefore offers high-precision benchmark solutions as well as dynamical insight.
As in the thixotropic Stokes problem [17], we find a rapidly adjusting limit in which the fluid has a generalized-Newtonian rheology, and a slowly adjusting limit in which the structure depends on position but not on time and is controlled by a balance between the period-averaged buildup and breakdown rates. As the Deborah number increases, the structure of the fluid first tracks its instantaneous equilibrium value; then lags it, with the lag largest when the flow conditions are changing most rapidly; and finally approaches the slowly adjusting limit. For shear-thinning fluids the flux also lags behind the pressure gradient, most strongly for intermediate values of the Deborah number.
For a thixoviscoplastic Houška fluid, there are up to three flow regions: a fully structured core region; a partly structured but unyielded plug region; and an annular yielded region. The boundary between the plug and yielded regions moves in and out over the course of an oscillation; this behavior recalls that predicted in Stokes' third problem for a fluid with distinct static and dynamic yield stresses [26], which can be thought of as a limiting case of thixotropy.
Under asymmetrical sawtooth forcing with a longer period of acceleration in one direction than the other, thixotropic pumping leads to a net flux of fluid over the course of one period. For intermediate Deborah numbers and for strongly asymmetric waveforms, the total volume of fluid transported over a period can be a substantial fraction of the total volume of fluid transported in each direction over the corresponding half-period.
A purely thixoviscous fluid, described using the simplified Moore-Mewis-Wagner model, behaves analogously to a Houška fluid when it is shear-thinning. However, a shear-thickening fluid behaves in the opposite manner: the flux leads the pressure gradient for intermediate values of the Deborah number, and the net transport of a shear-thickening fluid under sawtooth forcing is in the opposite direction to that of a shear-thinning fluid.
It is plausible that the mechanisms elucidated here are generic rather than being artefacts of particular rheological models. Interesting directions for future work might include exploring how they are modified in more complex rheological models, such as those that exhibit viscosity bifurcations [27], and investigating the differences or similarities between the responses to asymmetrical forcing of thixotropic and elastoviscoplastic fluids [15].