2. Experimental viscosity measurements

Here, we consider the transient viscosity measurements of glycerol in a parallel-plate rheometer where the outer fluid interface is exposed to the atmosphere. Due to the hygroscopic nature of glycerol, it will absorb water vapour from the atmosphere at the outer boundary as shown in figure 1, which will subsequently lead to a local reduction in the viscosity of the fluid and a net reduction of the torque measured by the rheometer. Thus, the measured viscosity of glycerol in a parallel-plate rheometer is expected to decrease with time when exposed to atmosphere. Intuitively, the rate of decrease should depend on the water concentration/flux experienced at the outer boundary, which is influenced by the RH in the atmosphere. The rate of decrease should also depend on the secondary recirculating flows in the rheometer, since these redistribute the relatively less viscous fluid where the water concentration is higher, thereby altering the stress distribution and total torque experienced by the upper plate of the rheometer.

2.2. Experimental results

The experimental results for the transient viscosity measurements are presented in figure 2. First, the transient viscosity measurements of glycerol in a RH environment of 53 % are presented in figure 2(a). The experiments were performed with a rheometer with plate radius $R=2.5$ cm at an angular velocity of $\varOmega =0.4$ rad s$^{-1}$. Results are normalized by the initial viscosity $\mu _i$ to account for any small amount of water absorption by the glycerol while setting up the experiment. As can be seen, the measured viscosities decrease at varying rates depending on the gap thickness. These trends are all monotonic. However, the rate of viscosity decrease is seen to be a non-monotonic function of the gap thickness. That is, decreasing the gap height from $h_0=1.0$ to $h_0=0.5$ mm results in the viscosity dropping at a slower rate. However, subsequently decreasing the gap height results in the viscosity dropping at a faster and faster rate, until the viscosity drops sharply for the smallest gap height. Furthermore, experimental results for the final measured viscosity $\mu _f/\mu _i$ after $t=60$ minutes are shown in figure 2(b) as a function of gap height $h_0$ for two different RHs. Here, with ${\rm RH}=72\,\%$ the viscosities are seen to drop faster than with ${\rm RH}=45\,\%$ due to the increased mass flux of water vapour to the glycerol from the atmosphere. In addition, the final measured viscosity values also show the same non-monotonic behaviour as in figure 2(a). That is, decreasing $h_0$ first leads to an increase in $\mu _f/\mu _i$, followed by a sharp decrease for $h_0<0.5$ mm.

Figure 2.

Experimental viscosity measurements. (a) Transient decrease in the measured viscosity as a function of time, normalized by the initial viscosity $\mu _i$ with $\varOmega =0.4$ rad s$^{-1}$ and $R=2.5$ cm. Experiments were performed at a RH of $72\,\%$. (b) Measured final viscosities $\mu _f$ of glycerol at $t=3600$ s normalized by $\mu _i$ with $\varOmega =0.4$ rad s$^{-1}$ and $R=2.5$ cm. Experiments were performed at RHs of $45\,\%$ and $72\,\%$. (c) Normalized transient viscosity measurements over time for varying RHs at a gap height of 0.1 mm and a rotation rate of 0.4 rad s$^{-1}$. (d) Normalized transient viscosity measurements over time for varying rotation rates at a gap height of 0.5 mm and a RH of 65 %.

Next, the transient decreases in measured viscosity values are shown in figure 2(c) as functions of RH at a rotation speed of 0.4 rad s$^{-1}$ and a gap height of 0.1 mm. Here the results show a monotonic relationship in which increasing the RH leads to a faster decrease in viscosity measurements due to the increased flux of water into the glycerol. Finally, the transient measured viscosities are shown in figure 2(d) for varying rotation speeds at a gap height of 0.5 mm and a RH of 65 %. Here, a non-monotonic relationship with $\varOmega$ is observed. As the rotation speed is increased from 1.0 to 5.0 rad s$^{-1}$, the viscosity drops at a faster rate, which is consistent with the increasing inertial secondary flow in this regime. The key observations from the experimental measurements are as follows.

1. (i) Increased RH leads to a faster viscosity decrease for all measured cases.

2. (ii) The rate of viscosity decrease varies non-monotonically with the gap height. That is to say, for different gap heights, the measured viscosity decreases at a different rate over time.

3. (iii) There is a sharp decrease in the measured final viscosity for very small gap heights.

Observation (i) above is a natural and intuitive result, since the hygroscopic nature of glycerol leads it to absorb more water from the atmosphere at higher humidities. Since the viscosity of a glycerol–water mixture varies monotonically with the water mass fraction, the rate of viscosity decrease can be expected to vary monotonically with the RH. However, an intuitive explanation for observations (ii) and (iii) is not immediately obvious without a more careful consideration of the fluid dynamics and the coupled transport of dissolved water in the rheometer. In the following sections, we seek a physical explanation for these two experimentally observed behaviours using theory and numerical simulations.

3. Fluid dynamics of a glycerol–water mixture

Before proceeding to analyse the specific fluid dynamics and transient viscosity evolution in a parallel-plate rheometer, in this section we first introduce all of the relevant governing physics that applies to such systems. This includes the governing Navier–Stokes equations for variable-viscosity fluids, the advection–diffusion equation for the absorbed water concentration with variable diffusivity, the empirical relationships for the coefficients of viscosity and diffusivity of glycerol–water mixtures as well as the saturation concentration of absorbed water in glycerol as a function of the RH that will be used in the modelling.

3.1. Navier–Stokes equations for a variable-viscosity fluid

In a glycerol–water mixture, the viscosity is strongly dependent on the local mass fraction of water. In systems with a homogeneous concentration of absorbed water, the constant-viscosity form of the Navier–Stokes equations is appropriate. However, when gradients in water concentration exist in the system, due to circumstances such as mixing streams or the absorption of water at a gas–liquid interface, the viscosity of the fluid must be treated as a function of time and position. Note that for glycerol–water mixtures, the density is a weak function of the water concentration, ranging from approximately 1260 kg m$^{-3}$ for pure glycerol down to 1000 kg m$^{-3}$ for pure water. Here, we neglect this variation and use the density of pure glycerol, assuming that the water mass fraction does not get too large. The corresponding continuity equation for an incompressible flow is simply $\boldsymbol {\nabla }^*\boldsymbol {\cdot }\boldsymbol u^* = 0$. In such a system, extra stresses arise in the fluid that are related to the gradients of viscosity, and the appropriate form of the Navier–Stokes equations (for a Newtonian constitutive equation) is

(3.1a,b)\begin{gather} \boldsymbol{\nabla}^*\boldsymbol{\cdot }\boldsymbol u^* = 0 \quad \text{and}\quad \rho\frac{{\rm D}\boldsymbol{u}^*}{{\rm D} t^*}={-}\boldsymbol{\nabla}^* p^* + \mu^* {\nabla^*}^2 \boldsymbol u^* + \boldsymbol{\nabla}^*\mu^*\boldsymbol{\cdot }\boldsymbol{\nabla}^* \boldsymbol u^* + \boldsymbol{\nabla}^*\mu^*\boldsymbol{\cdot } \left(\boldsymbol{\nabla}^*\boldsymbol u^*\right)^\mathsf{T}, \end{gather}

where $\rho$ is the fluid density, $\mu$ is the fluid viscosity and superscript $\mathsf {T}$ denotes the transpose. Here, the asterisks denote dimensional variable quantities. For the case of the flow in a rheometer system, we represent the flow using cylindrical coordinates, and we non-dimensionalize the governing equations with

(3.2a–h)\begin{equation} \left.\begin{gathered} r=\frac{r^*}{R},\quad z=\frac{z^*}{h_0},\quad u_r=\frac{u_r^*}{\varOmega R},\quad u_\theta=\frac{u_\theta^*}{\varOmega R},\quad u_z=\frac{u_z^*}{\varOmega h_0},\\ p=\frac{p^*}{\mu_0 \varOmega R^2/h_0^2},\quad \mu=\frac{\mu^*}{\mu_0}\quad \text{and}\quad t=\varOmega t^*, \end{gathered}\right\} \end{equation}

where $\mu _0$ is a reference viscosity value. With these non-dimensionalizations, the component form of (3.1a,b) in cylindrical coordinates is given by

(3.3a)\begin{align} & {Re}\left(\frac{\partial u_r}{\partial t}+u_r\frac{\partial u_r}{\partial r}+ \frac{u_\theta}{r}\frac{\partial u_r}{\partial\theta}-\frac{{u_\theta}^2}{r}+u_z \frac{\partial u_r}{\partial z}\right)\nonumber\\ &\quad ={-}\frac{\partial p}{\partial r} + \mu\left(\frac{\epsilon^2}{r}\frac{\partial}{\partial r}\left(r\frac{\partial u_r}{\partial r}\right) +\frac{\epsilon^2}{{r}^2}\frac{\partial^2u_r}{\partial\theta^2} +\frac{\partial^2u_r}{\partial {z}^2}-\epsilon^2 \frac{u_r}{{r}^2}-\frac{2\epsilon^2}{{r}^2}\frac{\partial u_\theta}{\partial\theta}\right) \nonumber\\ &\qquad\ +2\epsilon^2\frac{\partial\mu}{\partial r}\frac{\partial u_r}{\partial r} + \frac{\epsilon^2}{r} \frac{\partial\mu}{\partial\theta}\left(\frac{1}{r}\frac{\partial u_r}{\partial\theta}+ \frac{\partial u_\theta}{\partial r}-\frac{u_\theta}{r}\right) + \frac{\partial \mu}{\partial z}\left(\frac{\partial u_r}{\partial z}+\epsilon^2 \frac{\partial u_z}{\partial r}\right), \end{align}

(3.3b)\begin{align} & {Re}\left(\frac{\partial u_\theta}{\partial t}+u_r\frac{\partial u_\theta}{\partial r}+ \frac{u_\theta}{r}\frac{\partial u_\theta}{\partial\theta}+\frac{u_\theta u_r}{r}+u_z \frac{\partial u_\theta}{\partial z}\right) \nonumber\\ &\quad ={-}\frac{1}{r}\frac{\partial p}{\partial\theta} +\mu\left(\frac{\epsilon^2}{r}\frac{\partial}{\partial r}\left(r\frac{\partial u_\theta}{\partial r}\right)+ \frac{\epsilon^2}{{r}^2}\frac{\partial^2 u_\theta}{\partial\theta^2}\right. \nonumber\\ &\qquad \left.+\frac{\partial^2u_\theta}{\partial {z}^2}-\epsilon^2\frac{u_\theta}{{r}^2}+ \frac{2\epsilon^2}{{r}^2}\frac{\partial u_r}{\partial \theta}\right)+\epsilon^2 \frac{\partial\mu}{\partial r}\left(\frac{\partial u_\theta}{\partial r}+\frac{1}{r} \frac{\partial u_r}{\partial \theta}-\frac{u_\theta}{r}\right) \nonumber\\ &\qquad + \frac{2\epsilon^2}{r}\frac{\partial\mu}{\partial\theta}\left(\frac{1}{r} \frac{\partial u_\theta}{\partial\theta}+\frac{u_r}{r}\right) + \frac{\partial \mu}{\partial z}\left(\frac{\partial u_\theta}{\partial z}+ \frac{\epsilon^2}{r}\frac{\partial u_z}{\partial\theta}\right), \end{align}

(3.3c)\begin{align} & {Re}\,\epsilon^2\left(\frac{\partial u_z}{\partial t}+u_r\frac{\partial u_z}{\partial r}+ \frac{u_\theta}{r}\frac{\partial u_z}{\partial\theta}+u_z\frac{\partial u_z}{\partial z}\right)\nonumber\\ &\quad ={-}\frac{\partial p}{\partial z} +\mu\left(\frac{\epsilon^4}{r}\frac{\partial}{\partial r} \left(r\frac{\partial u_z}{\partial r}\right)+\frac{\epsilon^4}{{r}^2} \frac{\partial^2u_z}{\partial\theta^2}+\epsilon^2\frac{\partial^2u_z}{\partial {z}^2}\right) \nonumber\\ &\qquad \, + \frac{\partial\mu}{\partial r}\left(\epsilon^4\frac{\partial u_z}{\partial r}+\epsilon^2 \frac{\partial u_r}{\partial z}\right)+\frac{1}{r}\frac{\partial\mu}{\partial\theta}\left(\frac{\epsilon^4}{r} \frac{\partial u_z}{\partial\theta}+\epsilon^2\frac{\partial u_\theta}{\partial z}\right) + 2\epsilon^2\frac{\partial\mu}{\partial z}\frac{\partial u_z}{\partial z}, \end{align}

where the Reynolds number is defined as ${Re}=\rho \varOmega h_0^2/\mu _0$ and the gap aspect ratio is $\epsilon =h_0/R$, which is typically small in a parallel-plate rheometer. The corresponding continuity equation is given by

(3.4)\begin{equation} \frac{1}{r}\frac{\partial}{\partial r}\left(r u_r\right) + \frac{1}{r} \frac{\partial u_\theta}{\partial \theta} + \frac{\partial u_z}{\partial z} = 0. \end{equation}

Along with boundary conditions, these equations govern the motion of variable-viscosity fluids in a parallel-plate rheometer. The typical boundary conditions for such a system include no-slip conditions at the lower (stationary) and upper (rotating) plates, as well as a stress-free condition at $r=1$ at the gas–liquid interface.

3.2. Water absorption and transport

Along with the equations governing the fluid dynamics in the previous section, the system further requires a transport equation to model the absorption and transport of water in the glycerol. In particular, this transport can be modelled with an advection–diffusion equation that is given by

(3.5)\begin{equation} \frac{\partial c}{\partial t^*}=\boldsymbol{\nabla}^* \boldsymbol{\cdot } \left(D^*\boldsymbol{\nabla}^* c\right) - \boldsymbol{u}^*\boldsymbol{\cdot } {\boldsymbol{\nabla}^*} c, \end{equation}

where $c$ is the mass fraction of water in the glycerol, $D^*$ is the diffusivity of water in glycerol and $\boldsymbol {u}^*$ is the dimensional fluid velocity vector. Here, $D^*$ is a spatially/temporally varying function of the local water concentration $c$. Using the same non-dimensionalizations as above, (3.5) in cylindrical coordinates becomes

(3.6)\begin{align} Pe\,\epsilon^2\frac{\partial c}{\partial t} &= \frac{\epsilon^2}{r} \frac{\partial}{\partial r}\left(rD\frac{\partial c}{\partial r}\right)+ \frac{\epsilon^2}{r^2}\frac{\partial}{\partial\theta}\left(D\frac{\partial c}{\partial\theta}\right)+ \frac{\partial}{\partial z}\left(D\frac{\partial c}{\partial z}\right) \nonumber\\ &\quad -Pe\,\epsilon^2 \left(u_r\frac{\partial c}{\partial r}+\frac{u_\theta}{r} \frac{\partial c}{\partial \theta} + u_z\frac{\partial c}{\partial z}\right), \end{align}

where $D=D^*/D_0$ is the non-dimensional diffusivity with reference value $D_0$ and $Pe=\varOmega R^2/D_0$ is the Péclet number representing the ratio of the time scale for diffusion of water in the radial direction to the convective time scale. Examining (3.6), the diffusion in the $z$ direction is ${O}(\epsilon ^{-2})$ larger than the radial diffusion due to the separation in length scales: in the low-inertia, thin-gap limit, the water concentration will be approximately uniform in the $z$ direction.

Along with (3.6), the absorbed water concentration must satisfy certain boundary conditions in the system. In particular, the concentration satisfies a no-flux condition at both the upper and lower plates of the rheometer (i.e. ${\partial c}/{\partial z}=0$ in the parallel-plate case). In addition, a boundary condition for the water concentration is needed at the outer glycerol–air interface. In general, this condition could be represented by a water flux condition such as by $-D({\partial c}/{\partial r})|_{r=1}=j_w$, where the flux of water $j_w$ could be a function of the local water concentration in the glycerol at the interface as well as the water vapour concentration and distribution in the air near the interface. The solution of such a flux will typically require also solving the water vapour transport problem in the surrounding environment, since these transport processes are coupled at the interface. For example, the transport of water vapour in the surrounding atmosphere may be affected by the rotation speed of the rheometer, which can drive flow in the surrounding air that may alter the vapour transport at the interface. Here, we neglect these effects and assume that the water vapour transport in the atmosphere is fast relative to the water concentration transport in the glycerol. This assumption is valid due to the substantially higher diffusivity of water vapour in air than absorbed water in glycerol, provided the recirculation in the rheometer is not significantly fast. Thus, we assume that the water mass fraction at $r=1$ instantaneously reaches its saturation value based on the local RH in the surrounding air.

Considering both (3.3) and (3.6), we see how the dynamics of the problem are fully coupled. In the fluid problem, the Navier–Stokes equations are coupled to the water transport problem via the dependence of viscosity on water mass fraction. The water transport problem is further coupled to the Navier–Stokes equations through the dependence on the flow velocity. Finally, the transport is further complicated by the dependence of diffusivity on water mass fraction. Due to this two-way coupling and the empirical nature of both the $\mu (c)$ and $D(c)$ relationships, it is difficult to seek theoretical solutions to the coupled dynamics. Thus, here we primarily rely on numerical simulations to solve the coupled transport problem.

3.3. Physical properties of glycerol–water mixtures

Detailed empirical formulas for the viscosity of a glycerol–water mixture have been proposed by Cheng (Reference Cheng2008) which are valid for water mass concentrations in the range of 0 %–100 % and for temperatures ranging from 0 to 100 $^{\circ }$C. The viscosity of a glycerol–water mixture at 22 $^{\circ }$C varies from around $\mu _g=1.1$ Pa s for pure glycerol down to $\mu _w=0.96$ mPa s for pure water, spanning a range of approximately three orders of magnitude. Here, we choose a reference viscosity value corresponding to that of pure glycerol, $\mu _0=\mu _g$, such that the non-dimensional viscosity varies from an initial condition of 1 down to as low as ${\sim }8.73\times 10^{-4}$ if the water mass fraction were to approach 1. The empirical relationship determined between non-dimensional viscosity and water concentration used here is shown in figure 3(a).

Figure 3.

(a) Empirical non-dimensional viscosity and diffusivity relationships used in this study. (b) Saturation concentration of water in glycerol and corresponding specific gravity as functions of RH.

The diffusivity of water in glycerol is also a function of the water mass fraction, and so $D^*$ is expected to evolve as a function of both position and time as more water is absorbed and transported throughout the system. An empirical relationship for the diffusivity of water in glycerol has also been developed for mixtures at 25 $^{\circ }$C by D'Errico et al. (Reference D'Errico, Ortona, Capuano and Vitagliano2004), and is given by

(3.7)\begin{equation} D^*=\frac{1.024-0.91x_g}{1+7.5x_g}\times 10^{{-}9}\,{\rm m}^2\,{\rm s}^{{-}1}, \end{equation}

where $x_g$ is the mole fraction of glycerol that is related to $c$ via

(3.8)\begin{equation} x_g=\frac{M_w(1-c)}{M_w(1-c)+M_gc}, \end{equation}

where $M_w$ is the molar mass of water and $M_g$ is the molar mass of glycerol. For reference, the diffusivity of water in pure glycerol ($c=0$) is $1.341\times 10^{-11}$ m$^2$ s$^{-1}$, which increases to approximately $1.024\times 10^{-9}$ m$^2$ s$^{-1}$ as the water mass fraction approaches 1. Here we again choose a reference value equal to the diffusivity in pure glycerol $D_0=D_g$, such that $D$ varies from 1 at $c=0$ up to approximately 76.4. This empirical relationship from D'Errico et al. (Reference D'Errico, Ortona, Capuano and Vitagliano2004) is also shown in figure 3(a). Finally, the saturation concentration of water in glycerol as a function of RH is needed for the absorption boundary condition at $r=1$. These values were measured by the Glycerine Producers’ Association and are given in table 15 of Glycerine Producers’ Association et al. (1963). For reference, these values are plotted in figure 3(b) along with the specific gravity as functions of RH. With the full governing equations and empirical relationships for the physical parameters described above, we can now move on to consider the coupled fluid dynamics and water concentration transport in a rheometer. In the following section we first review the classical result for the axisymmetric parallel-plate rheometer with constant viscosity before moving on to cases with variable viscosity.

3.4. Note on the assumption of constant density

Before moving on, we briefly comment on the assumption of constant density. As water is absorbed into the glycerol, the resulting density gradients introduce the possibility for buoyancy-driven flows. An estimate for the magnitude of such effects is given by considering that a vertical change in density ${\rm \Delta} \rho$ implies a radial pressure gradient of the order of ${\rm \Delta} \rho g h_0/R$. We can balance this radial pressure gradient with a radial viscous stress gradient $\mu u_b/h_0^2$, where $u_b$ is the characteristic magnitude of the buoyancy-driven flow. Thus, we have $u_b\sim {{\rm \Delta} \rho g h_0^3}/{\mu R}$. For buoyancy effects to be negligible, we need $u_b$ to be small relative to the magnitude of the inertial secondary velocity components which are ${O}({\rho \varOmega ^2h_0^2 R}/{\mu })$, as shown in (1.1). Thus, we need

(3.9)\begin{equation} \frac{{\rm \Delta}\rho g h_0^3}{\mu R} \ll \frac{\rho\varOmega^2h_0^2 R}{\mu}\longrightarrow \frac{{\rm \Delta}\rho}{\rho} \ll \frac{\varOmega^2 R^2}{g h_0}. \end{equation}

The left-hand side of this inequality has a maximum value of around 0.2 when the saturation water concentration approaches 100 %, and so it will typically take a smaller value depending on the RH. The right-hand side of the inequality depends on the system parameters, but a typical value with $\varOmega = 4.0$ rad s$^{-1}$, $R=2.5$ cm and $h_0=0.5$ mm is approximately 2.0. Thus, the buoyancy-driven flow can be expected to be less than 10 % of the inertial secondary flow. At larger gap heights and smaller rotation speeds, this inequality suggests that buoyancy effects become relatively more important, and could even dominate the dynamics. However, in the very slow rotation case, the dynamics are nearly one-dimensional (1-D), such that the density variation is almost entirely in the radial direction and the negligible vertical density gradient should not affect the radial pressure gradient. So, (3.9) should be considered only as an estimate of the importance of buoyancy effects.

4. Classical result: axisymmetric flow with constant viscosity

Before moving on to consider the full coupled dynamics of glycerol absorbing water in a parallel-plate rheometer, we first review the classical Newtonian, constant-viscosity flow solution in a parallel-plate rheometer. Understanding this flow is important because it illustrates the types of secondary flows we should expect in the rheometer, and also provides some initial insights into the non-monotonic relationship between measured viscosities and rotation speed and gap height in the full system. The axisymmetric parallel-plate rheometer has a well-known solution that has been previously described by multiple authors (see e.g. Middleman Reference Middleman1968; Bird et al. Reference Bird, Armstrong and Hassager1987). Here, we briefly reproduce this calculation in our notation for consistency with later sections. We consider a parallel-plate rheometer with radius $R$ and gap thickness $h_0$. The lower plate is stationary and the upper plate rotates at a rotation speed of $\varOmega$. We model the system using cylindrical coordinates with the origin located at the centre of the bottom plate. Thus, the governing equations are the axisymmetric and constant-viscosity forms of (3.3) and (3.4). The corresponding boundary conditions are no-slip at both the upper and lower plates, i.e. $(u_{r,{axi}},u_{\theta,{axi}},u_{z,{axi}})=(0,0,0)$ at $z=0$ and $(u_{r,{axi}},u_{\theta,{axi}},u_{z,{axi}})=(0,r,0)$ at $z=1$. For small gap heights $\epsilon \ll 1$ a solution for the velocity and pressure can be sought in the form of an expansion in powers of $\epsilon ^2$. Up to ${O}(\epsilon ^4)$, this axisymmetric solution is given by

(4.1a)\begin{gather} u_{r,{axi}}(r,z)={-}\frac{1}{12}r{Re}\,z(z-1)\left(-\frac{4}{5}+z+z^2\right)+{O}(\epsilon^4), \end{gather}

(4.1b)\begin{gather}u_{\theta,{axi}}(r,z)=rz - \frac{r{Re}^2z}{6300}(8+z^3(35-63z+20z^3)) + {O}(\epsilon^4), \end{gather}

(4.1c)\begin{gather}u_{z,{axi}}(z)=\frac{1}{30}{Re}\,z^2(z-1)^2(2+z)+{O}(\epsilon^4), \end{gather}

(4.1d)\begin{gather}p_{axi}(r)=\frac{3r^2}{20}{Re} + \frac{1}{30}{Re}\,\epsilon^2 z(4-9z+5z^3) + {O}(\epsilon^4). \end{gather}

Here, the subscript ‘$axi$’ denotes the axisymmetric case, and the expression for $u_{r,{axi}}$ is equivalent to the result first presented by Savins & Metzner (Reference Savins and Metzner1970), which was given above in dimensional form as (1.1). The primary flow is the $u_{\theta,{axi}}=rz$ component with an ${O}({Re}^2)$ correction, while the leading-order secondary flows in the $r$ and $z$ directions are both ${O}({Re})$. Since the flow of interest is axisymmetric, the primary velocity components of interest for redistributing absorbed species at the outer edge of the rheometer are the secondary velocity components, especially the radial component, since this will transport absorbed water from the outer edge inwards through the gap.

A visualization of the secondary velocity components is given in figure 4. Here, the $u_{z,{axi}}$ component shows that there is an upward drift that is independent of $r$ throughout the rheometer. The radial component shows that in the upper half of the gap the flow is directed radially outwards, and in the lower half of the gap the flow is directed radially inwards. Keep in mind that the theoretical solution presented in (4.1) must break down near the outer edge of the rheometer where $r\rightarrow 1$, since the lubrication approximation fails in that region. In the true system, near the outer edge the outward-travelling flow in the upper half of the rheometer must turn downwards for continuity and turn around to then travel radially inwards. The width of this turning region should be ${O}(\epsilon )$, and thus is progressively more confined at the outer edge as the gap height decreases. This effect is not captured by (4.1), although it may have an important role in the transport of absorbed water in the glycerol–water system.

Figure 4.

Secondary velocity components $u_{r,{axi}}$ and $u_{z,{axi}}$ in a parallel-plate rheometer in the small-gap limit ($\epsilon \ll 1$). The radial velocity proceeds outward along the upper half of the gap, reverses at the outer edge and proceeds radially inward along the lower half of the gap. The $z$ component shows an upward drift along the middle of the gap that is independent of $r$. Both secondary velocities are ${O}({Re})$. Results correspond to (4.1). Note that $u_{r,{axi}}$ is plotted normalized by $r$ and ${Re}$, and $\boldsymbol {u}_{axi}$ satisfies incompressibility.

This picture of the secondary flows appears to yield a satisfactory understanding of the relationship between the evolving viscosity measurements and the gap height and rotation speed, at least when the gap height is not too small. That is, for gap heights of 0.5 mm and above, increasing either the gap height or the rotation speed will contribute to greater radial secondary flows that pull relatively higher-water-concentration glycerol away from the outer edge, enhancing the flux of water into the glycerol at the outer edge. This contributes to a faster decrease in the measured viscosity as the high-water-concentration glycerol at the outer edge generates less shear stress on the upper plate. However, this picture of the secondary flows is inconsistent with our observed viscosity measurements at small gap heights. In the limit of ${Re}\ll 1$, the secondary flows in a parallel-plate rheometer are negligible. In this case, the transport equation simplifies to a purely 1-D radial diffusion problem, and the evolution of water concentration becomes independent of both the gap height and rotation speed. Furthermore, the viscosity distribution likewise is independent of $h_0$ and $\varOmega$, which is inconsistent with the sharp decrease in $\mu _f/\mu _i$ seen in figure 2(b) at very small gap heights. This will motivate us later in the paper to consider the potential role of misalignment. First, we examine in more detail the 1-D diffusive limit with variable viscosity. As a quick point of reference, with the definition of Reynolds number given by ${Re}=\rho _g\varOmega h_0^2/\mu _g$, with the characteristic density and viscosity based on values for pure glycerol, the experimental results presented above in figure 2 have Reynolds numbers ranging from ${\sim }1\times 10^{-6}$ up to ${\sim }0.05$. While these values seem small, recall that the dimensional viscosity can vary by over three orders of magnitude, such that a locally defined Reynolds number could be significantly larger.

5. One-dimensional diffusive limit with variable viscosity

First, we consider the evolution of the viscosity distribution and measured effective viscosity of glycerol absorbing water in the inertialess, 1-D diffusive limit. For small Reynolds numbers and gap heights, the axisymmetric form of (3.6) becomes

(5.1)\begin{equation} Pe\frac{\partial c}{\partial t}=\frac{1}{r}\frac{\partial}{\partial r} \left(rD\frac{\partial c}{\partial r}\right)+\frac{1}{\epsilon^2} \frac{\partial}{\partial z}\left(D\frac{\partial c}{\partial z}\right). \end{equation}

Considering that the water concentration boundary condition at $r=1$ is independent of $z$, along with the no-flux conditions at the upper and lower plates, when $\epsilon \ll 1$ it must be the case that $c$ is approximately independent of $z$, so that (5.1) further simplifies to

(5.2)\begin{equation} Pe\frac{\partial c}{\partial t}=\frac{1}{r}\frac{\partial}{\partial r} \left(rD\frac{\partial c}{\partial r}\right), \end{equation}

which is simply a 1-D radial diffusion equation with variable diffusivity $D(c)$. In this regime, a better choice for the characteristic time scale would be the characteristic radial diffusion time $R^2/D_0$, the use of which would yield the same equation without the $Pe$ factor on the left-hand side. For consistency, we continue to use the convective $1/\varOmega$ time scale as the characteristic time scale. Here, the only boundary conditions that are needed are symmetry at $r=0$ and the saturation water mass fraction at $r=1$, i.e. $c(r=1)=c_{sat}$.

As the water concentration evolves, the anticipated viscosity measurement from the rheometer can be predicted through the use of the viscosity distribution as follows. A parallel-plate rheometer cannot measure the viscosity distribution throughout the fluid layer, but rather simply infers an effective viscosity $\mu _{eff}^*$ by measuring the total torque exerted on the upper plate as it spins. In dimensional form, the azimuthal velocity at small gap heights is $u_\theta ^*=\varOmega r^* z^*/h_0$. This velocity profile is valid regardless of the viscosity distribution since $c$ is a function of $r$ and $t$ only. The total torque experienced by the upper plate is then given by

(5.3)\begin{equation} T=\int_0^{2{\rm \pi}}\int_0^R\frac{\mu^*\varOmega}{h}{r^*}^3\, \mathrm{d}r^*\,\mathrm{d}\theta. \end{equation}

If the viscosity is constant and uniform, the total torque on the upper plate is then

(5.4)\begin{equation} T = \frac{{\rm \pi} \mu^* \varOmega R^4}{2h}\longrightarrow \mu^*=\frac{2hT}{{\rm \pi}\varOmega R^4}. \end{equation}

The rheometer assumes a constant-viscosity fluid and reports the ‘effective’ viscosity of the fluid that is calculated from (5.4) based on the measured torque. In the experimental system, the initial condition is assumed to be pure glycerol, such that $T_{init}={\rm \pi} \mu _g\varOmega R^4/(2h)$, and we have

(5.5) \begin{align} & \frac{\mu_{eff}^*}{\mu_g}=\frac{T(t)}{T_{init}}= \frac{2h}{{\rm \pi}\mu_g\varOmega R^4}\int_0^{2{\rm \pi}}\int_0^R \frac{\mu^*\varOmega}{h}{r^*}^3\,\mathrm{d}r^*\,\mathrm{d}\theta \nonumber\\ &\quad\longrightarrow \mu_{eff}=\frac{2}{\rm \pi}\int_0^{2{\rm \pi}}\int_0^1 \mu{r}^3\, \mathrm{d}r\,\mathrm{d}\theta, \end{align}

which simplifies to

(5.6)\begin{equation} \mu_{eff} = 4\int_0^1 \mu r^3\, \mathrm{d}r \end{equation}

for axisymmetric flow.

Here, we perform 1-D transient simulations of (5.1) using the finite-difference method with second-order accuracy in space and first-order accuracy in time. Convergence studies were performed in space and time to verify the results. Using these simulations, we compute the effective non-dimensional viscosity over time as water is absorbed at the outer edge and diffuses radially inwards. We perform these simulations over a range of $c_{sat}$ values which reproduces the effect of varying RHs. First, for comparison with experiments, the results are simulated for one hour to determine the degree of viscosity decrease that can be achieved via pure diffusion over the experimental time scale. These results are shown in figure 5(a). As can be seen, diffusion alone is sufficient to generate a significant decrease in measured viscosities over this time scale, although not to the degree seen in the experiments. For example, consider the experimental results in figure 2(c), which were performed at a Reynolds number of ${Re}=\rho \varOmega h_0^2/\mu _0=4.6\times 10^{-6}$ and aspect ratio of $\epsilon =h_0/R=4\times 10^{-3}$. Clearly, in such a regime the inertialess, 1-D model would be expected to apply. However, the experimental results show a much larger decrease in viscosity over this time scale. The $\textrm {RH}=72\,\%$ results (corresponding to approximately $c_{sat}=0.386$) drop to around $\mu =0.38$, and the $\textrm {RH}=45\,\%$ results (corresponding to $c_{sat}=0.185$) drop to around $\mu =0.6$. However, in the 1-D limit, the corresponding numerical predictions for $c_{sat}=0.4$ and $c_{sat}=0.2$ only decrease to around $\mu =0.78$ and $\mu =0.86$, respectively.

Figure 5.

Transient non-dimensional viscosity measurements in the inertialess, 1-D, axisymmetric regime. (a) Results simulated over 3600 s for comparison with experimental measurements. Results show substantial decreases in measured viscosities at large $c_{sat}$ values, but not as significant as those seen in the experiments. The inset shows the final non-dimensional viscosity $\mu _f$ versus $c_{sat}$. (b) Results extended to much longer times to show the final saturation of the glycerol, which corresponds to the curves levelling off. Clearly, higher values of $c_{sat}$ reach saturation more quickly.

Thus, the experiments show a much larger decrease in viscosity over this time scale than the 1-D model with pure diffusion. Furthermore, the results shown in figure 5(a) are clearly still evolving over this time scale, whereas in the long-time limit we expect all of the glycerol to homogenize at the saturation concentration based on the RH. Therefore, we extend these results to much longer times in figure 5(b), which shows that the measured effective viscosities level off as the water concentration saturates. Here we see the influence of the variable diffusivity on the time scale for the diffusive process. With the characteristic diffusivity $D_0$, the time scale for the process would be expected to be $t^*={O}(R^2/D_0)$. However, as can be seen, most of the cases have fully saturated well before this time scale, especially at larger $c_{sat}$ values. This is due to the enhanced diffusion at larger water concentrations. In fact, a much better prediction for the time scale of this 1-D diffusive process is to use the diffusivity based on $c_{sat}$, which we call $D_{sat}$. The rescaled results are shown in figure 6, which shows that for each case the water concentration in the glycerol has fully saturated over the time scale $t^*={O}(R^2/D_{sat})$. For each case, two regimes can be seen. In the early times, the water concentration in the glycerol is non-uniform, and so the diffusive transport in the domain proceeds with a spatially varying diffusivity coefficient. At late times, the water concentration throughout the system has nearly equilibrated at around the saturation concentration, such that the diffusion coefficient is nearly uniform and the results all decay exponentially with the same rate constant. This constant can be simply calculated by considering a 1-D radial diffusion problem with constant diffusivity (since this is nearly the case at long times), where the transport in dimensional form is governed by

(5.7)\begin{equation} \frac{\partial c}{\partial t^*}=D_{sat}\frac{1}{r^*}\frac{\partial}{\partial r^*} \left(r^*\frac{\partial c}{\partial r^*}\right)\quad \text{with} \left.\frac{\partial c}{\partial r^*}\right|_{r^*=0}=0\quad \text{and}\quad c(r^*=R)=c_{sat}. \end{equation}

The solution to this is given by

(5.8)\begin{equation} c(r^*,t^*)=c_{sat}+\sum_{n=1}^\infty a_n \exp({-D_{sat}t^*\lambda_n^2}) J_0(\lambda_n r), \end{equation}

where $a_n$ are coefficients that depend on the initial condition and $J_0$ is the zeroth-order Bessel function of the first kind. The $\lambda _n$ eigenvalues here are the roots of $J_0$ divided by the radius $R$. Thus, at late times we see that $c-c_{sat}\sim \exp (-\chi _1^2 t^* D_{sat} / R^2)$, where $\chi _1=2.40483$ is the first root of the $J_0$ function. Thus we see the $\chi _1^2=5.7832$ exponential decay seen in figure 6.

Figure 6.

Rescaled effective viscosities in the 1-D, inertialess, axisymmetric limit. For each case the water concentration fully saturates approximately over the time scale $t^*={O}(R^2/D_{sat})$, which is consistent with diffusion primarily occurring at the saturation concentration diffusivity. At late times, the rescaled viscosities all approach the saturation values exponentially with a rate constant of 5.78, consistent with the 1-D theory.

The previous results are strictly valid in the inertialess (${Re}\ll 1$), small-gap ($\epsilon \ll 1$) and axisymmetric limits. These calculations are significantly simplified compared with the solution for the inertial regime since (1) the water concentration profile can no longer be assumed to be independent of $z$ due to the secondary velocity components and (2) the fluid velocity profiles must be recalculated continuously as the concentration profile evolves while taking into account the spatial variations in viscosity. Nevertheless, it is clear that we must extend our results to the inertial regime, since the 1-D diffusion-dominated results cannot reproduce the same degree of viscosity decrease over the time scale of the experiments. These simulations are pursued in the following section.

6. Inertial regime with variable viscosity

Having explored the purely diffusion-dominated 1-D axisymmetric regime in the previous section, we now extend our results to the inertial, axisymmetric regime. Recall that in the inertial regime, the coupled dynamics is governed by four dimensionless parameters, which are

(6.1a–d)\begin{equation} Pe=\frac{\varOmega R^2}{D_0},\quad {Re}=\frac{\rho\varOmega h_0^2}{\mu_0},\quad \epsilon=\frac{h_0}{R}\quad \text{and}\quad c_{sat}, \end{equation}

whereas in the diffusion-dominated case the dynamics are governed only by $c_{sat}$. Thus, the system is governed by a relatively large parameter space. However, note that the ratio $\mu _0/(\rho D_0)=\mu _g/(\rho D_g)$ is fixed for glycerol, and the Péclet number can be written as

(6.2)\begin{equation} Pe=\frac{\varOmega R^2}{D_0}=\left(\frac{\rho\varOmega h_0^2}{\mu_0}\right) \left(\frac{R^2}{h_0^2}\right)\left(\frac{\mu_0}{\rho D_0}\right)=\frac{{Re}}{\epsilon^2} \left(\frac{\mu_0}{\rho D_0}\right), \end{equation}

so that the Péclet number is uniquely determined by the choice of ${Re}$ and $\epsilon$.

6.1. Numerical methods

Numerical simulations were performed using OpenFOAM (Weller et al. Reference Weller, Tabor, Jasak and Fureby1998) with an axisymmetric wedge-shaped mesh geometry with a wedge angle of 1$^{\circ }$. Local mesh refinement was used near $r=1$ to resolve the water concentration boundary layer. A sample mesh design is shown in figure 7. Simulations were performed using a custom in-house solver that iteratively updates the water concentration profile for 100 time steps using a time step of 0.01 s using second-order backward time-stepping and then recalculates the new steady-state velocity–pressure profiles using the SIMPLE algorithm (Jang, Jetli & Acharya Reference Jang, Jetli and Acharya1986; Barton Reference Barton1998). Thus, the solver assumes that the fluid velocity does not change much during one time step. Convergence tests were performed to confirm that recalculating the velocity every 100 time steps had a negligible impact on the calculated results compared with re-solving every time step. The solver assumes that the velocity–pressure profiles are quasi-steady and only evolve when the water concentration profile changes. The SIMPLE algorithm was used with relative pressure and velocity tolerances of $1\times 10^{-5}$. Convergence tests also confirmed the results were insensitive to these tolerances. Finally, grid resolution convergence tests were performed, and a final base grid of $375\times 30$ cells in the $r\times z$ directions was chosen. The cells within the region from $r=1-2\epsilon$ to 1 were all further refined by one level. Finally, the final layer of cells at $r=1$ was further refined by halving three times. Using this grid, convergence tests indicate that the errors due to spatial discretization should be less than 1 %. Torque measurements were calculated by integrating the wall shear stress over the upper plate.

Figure 7.

Sample computational mesh design for the inertial, axisymmetric simulations. The grid has been coarsened by a factor of 3 in the $r$ and $z$ directions for visualization purposes. Local mesh refinement is used near $r=1$ to resolve the water concentration boundary layer. (a) Top-down view of the axisymmetric wedge mesh geometry. (b) Side view of the wedge mesh. (c) Zoomed-in image of the local refinement near $r=1$. Several extra layers of very thin cells exist on the right-hand side which are difficult to see in order to resolve sharp concentration gradients that can occur at the boundary when inertial effects come into play.

6.2. Results

Using the numerical methods described in the previous section, simulations were performed across a range of gap heights, ${Re}$ and RHs (through their proxy $c_{sat}$). Before introducing the final measured viscosity values for comparison with the experiments, we first present results illustrating the evolution and dynamics of the water concentration field in the glycerol over a range of gap heights and rotation speeds. A comparison of the evolving water concentration profiles at various rotation speeds is shown in figure 8, and we give the results in dimensional form to make more clear the relationship of the changes to the experimental results presented earlier. These simulations were performed at $c_{sat}=0.2$ and $\epsilon =(1.0\times 10^{-3}\,\textrm {m})/(2.5\times 10^{-2}\,\textrm {m})=0.04$ with rotation speeds of 0.4, 1.0, 4.0 and 10.0 rad s$^{-1}$. The corresponding non-dimensional parameters for these cases are summarized in table 1. Here, the Reynolds number ranges from ${Re}=4.58\times 10^{-4}$ up to $1.15\times 10^{-2}$. This seems counter-intuitive, since even the smallest rotation speed case shows some deviation from a purely 1-D, diffusive transport, as can be seen by the concentration variation in the $z$ direction, whereas the relatively small Reynolds numbers suggest inertial effects should be small for all of these cases. However, consider that $u_{r,{axi}}^*\sim \varOmega R{Re}$. Then the characteristic time for convection in the radial direction is $\tau _{conv,rad}=(\varOmega \,{Re})^{-1}$, while the characteristic time for diffusion in the radial direction is $\tau _{diff,rad}=R^2/D_0$. Thus, an appropriate radial Péclet number is $Pe_{rad}={\tau _{diff,rad}}/{\tau _{conv,rad}}={\varOmega R^2{Re}}/{D_0}$. These values are also tabulated in table 1. As can be seen, even for the smallest angular velocity case with $\varOmega =0.4$ rad s$^{-1}$, the radial Péclet number is still greater than ${O}(10^3)$, increasing up to ${O}(10^6)$ at 10 rad s$^{-1}$. Thus, even at relatively small Reynolds numbers, the radial transport will be dominated by convection due to the relatively low diffusivity coefficients. For reference, table 1 also tabulates the non-dimensional parameters based on the viscosity and diffusivity values associated with the saturation concentration $c_{sat}$ rather than reference values based on pure glycerol. Here, the Reynolds numbers are increased while both the convective and diffusive time scales are decreased due to the reduced viscosity and increased diffusivity at increased water concentrations.

Figure 8.

Numerical results for the evolving water concentration profile $c$ over time for different rotation speeds $\varOmega$ at $c_{sat}=0.2$ and $\epsilon =(1.0\times 10^{-3}\,\textrm {m})/(2.5\times 10^{-2}\,\textrm {m})=0.04$. Here, the angular speeds are (a) $\varOmega =0.4$ rad s$^{-1}$, (b) $\varOmega =1.0$ rad s$^{-1}$, (c) $\varOmega =4.0$ rad s$^{-1}$ and (d) $\varOmega =10.0$ rad s$^{-1}$. The corresponding non-dimensional parameters are summarized in table 1. Here, the Reynolds number (based on the saturation viscosity rather than $\mu _g$) ranges from $9.71\times 10^{-3}$ up to 0.243 as the role of secondary (inertial) flows clearly grows with $\varOmega$.

Table 1.

Summary of the simulation parameters used in figures 8 and 10. The time-scale values $\tau _{conv}$ and $\tau _{diff}$ have units of seconds. Here, the parameters with ‘sat’ subscripts are calculated based on the fluid properties at the appropriate saturation mass fraction of water, and parameters without this subscript are calculated based on the fluid properties of pure glycerol.

Examining the transport dynamics in figure 8, we see that the transport of water concentration is consistent with the axisymmetric, constant-viscosity flow picture described above. In particular, in the constant-viscosity case, flow proceeds radially outward along the upper plate, turns downward and then flows radially inward along the lower plate. This is shown in more detail in figure 9(a). Here, the arrows are colour-coded and scaled by the magnitude of the secondary velocity components $|\boldsymbol u^*_{sec}|$, and the background is colour-coded by the water concentration profile. Here, $\boldsymbol u^*_{sec}$ is the velocity field on the slice in the $r$ and $z$ directions. As can be seen, the water begins to diffuse inwards from the outer edge, but then the secondary velocity pulls the absorbed water down along the outer edge and then radially inwards along the bottom plate. This creates a very thin boundary layer region near the outer edge of the rheometer, which gets thinner as $Pe_{rad}$ increases. The full solute concentration profile corresponding to figure 9(a) is shown in 9(b), and the corresponding dimensionless diffusivity, viscosity, radial velocity and $z$ velocity are shown in figures 9(c)–9(f), respectively. This case corresponds to the parameters previously shown in figure 8(d) with the parameters shown in table 1, and all results are at the non-dimensional time $t=\varOmega t^*=(10\,\textrm {rad}\,\textrm {s}^{-1})(1000\,\textrm {s})=1\times 10^4$.

Figure 9.

Detailed look at the glycerol–water dynamics for the case corresponding to figure 8(d) taken at $t=1\times 10^4$ (see table 1 for all relevant parameters). (a) Secondary velocity vectors coloured and scaled by the magnitude of the secondary velocity superimposed on a colourmap of the water concentration profile. As can be seen, the water begins to diffuse inward from the outer boundary, where the secondary flow pulls the absorbed water downward and then radially inward along the bottom plate, leading to a steep concentration gradient at the outer edge as the rotation speed is increased. (b) Full water concentration profile over a full axisymmetric cross-section. (c–f) Non-dimensional diffusivity coefficient, viscosity, radial velocity component and $z$ velocity component.

As can be seen in figure 9, the regions of high water concentration correspond to the regions of increased diffusivity and decreased viscosity. The radial velocity component resembles the flow for the axisymmetric constant-viscosity case with outward radial flow along the upper half of the domain and inward radial flow along the lower half, except that the magnitude of both is increased throughout the extent of the low-viscosity region. Furthermore, at the front of the propagating front of water concentration, there is a steep gradient in viscosity that corresponds to an upward secondary velocity due to the viscosity gradients as seen in figure 9(d,f).

Finally, a last illustration of the water concentration dynamics at higher $c_{sat}$ is presented in figure 10 as a function of gap height. These results show the evolution of the water concentration profile over time at $c_{sat}=0.5$ and $\varOmega =4.0$ rad s$^{-1}$ for gap heights ranging from 0.05 to 2.0 mm. Again, all of the relevant non-dimensional parameters are given in table 1 based on both the pure glycerol reference values and the saturation values. Here, the Reynolds number based on the saturation parameters ranges from $2.30\times 10^{-3}$ up to 3.68, representing a transition from the inertialess regime into the moderate inertial regime. The radial Péclet numbers based on the saturation properties remain relatively high, increasing from $1.46\times 10^4$ up to $2.34\times 10^7$ as the gap height increases, suggesting that the radial transport of water is dominated by convection in these regimes. Nevertheless, the transport at the smallest gap heights is approximately 1-D, suggesting that the $Pe_{rad}$ threshold for this transition is at a relatively large magnitude in this particular system.

Figure 10.

Numerical results for the evolving water concentration profile $c$ over time for different gap heights $h_0$ at $c_{sat}=0.5$ and $\varOmega =4.0$ rad s$^{-1}$. Here, the gap heights are (a) 0.05 mm, (b) 0.1 mm, (c) 0.2 mm, (d) 0.5 mm, (e) 0.75 mm, (f) 1.0 mm, (g) 1.25 mm, (h) 1.5 mm and (i) 2.0 mm. The corresponding non-dimensional numbers are summarized in table 1. Over this parameter range, the Reynolds numbers (based on the saturation viscosity) range from $2.30\times 10^{-3}$ up to 3.68 and the gap aspect ratio ranges from 0.002 to 0.08. Thus, these cases capture the full transition from the 1-D, diffusive limit up to the inertial regime.

Note that the qualitative picture of the water concentration evolution is different in figures 8 and 10. In particular, figure 10(f) corresponds to the same gap height and rotation speed as figure 8(c), except with an increased $c_{sat}$ value of 0.5 versus 0.2, respectively. While this seems like a relatively minor change, the corresponding $\mu _{sat}$ value is an order of magnitude smaller at $c_{sat}=0.5$, which leads to an order-of-magnitude stronger secondary flow in the region of locally low viscosity. This generates an enhanced mixing that leads to a more homogeneously propagating front of water concentration. This can be visualized in figure 11 for the case corresponding to figure 10(f) at early times. As can be seen, the local secondary recirculation in the low-viscosity region completely dominates the expected global secondary recirculation for the axisymmetric, constant-viscosity case. In fact, the global velocity field (the axisymmetric constant-viscosity solution) is negligible in the figure. This enhanced local recirculation at larger $c_{sat}$ values explains why the water propagates as a more uniform front in that regime, as opposed to being pulled down and inward along the lower plate as illustrated in figure 9(a) for a smaller $c_{sat}$ value.

Figure 11.

Enhanced secondary recirculation in the low-viscosity region corresponding to larger $c_{sat}$ values. These results illustrate the enhanced mixing effect that is seen at early times with $c_{sat}=0.5$ for the parameters shown in figure 10(f). At large values of $c_{sat}$, the local viscosity drops in regions of large $c$ to such a degree that the local recirculation dominates the expected secondary motions for constant-viscosity, axisymmetric flow.

Finally, having characterized and visualized the coupled transport dynamics in the parallel-plate, axisymmetric, inertial regime, we now calculate the measured effective viscosities in these simulations to see if the proposed model can fully capture the trends seen in the experimental results. The final measured dimensionless viscosities $\mu _f$ at $t^*=3600$ s are presented in figure 12 as functions of gap aspect ratio $\epsilon$, saturation concentration $c_{sat}$ and angular rotation speed $\varOmega$. The corresponding angular rotation speeds in the figure are 0.4, 1.0, 2.0 and 4.0 rad s$^{-1}$. In the figure, the dashed lines correspond to the predictions of the 1-D, axisymmetric, diffusion-dominated results previously described in figures 5 and 6. Note that in the experiments, the reported viscosities were non-dimensionalized by $\mu ^*_i$, the initial measured viscosity at $t=0$, and in the simulations the reported viscosities have been non-dimensionalized by $\mu _g$. Here, several key relationships and trends emerge from the results. First, we see clearly that in every case, the results approach the 1-D diffusion-dominated limit as $\epsilon \rightarrow 0$ for constant $\varOmega$, and they appear to also approach this limit as $\varOmega \rightarrow 0$ for constant $\epsilon$. For a given $\varOmega$, deviations from this limit increase as $\epsilon$ increases, due to enhanced inertial effects, as well as for increased $c_{sat}$. This latter trend is also due to an increase in inertial effects, although indirectly through a decreasing in the local viscosity. Furthermore, increasing $\varOmega$ clearly leads to more significant deviations from the 1-D limit due to increasing secondary inertial flows.

Figure 12.

Compilation of all measured final dimensionless viscosities $\mu _f$ at $t^*=3600$ s from the axisymmetric, inertial, variable-viscosity simulations for comparison with the experimental results. Dashed lines indicate the 1-D, inertialess, diffusion-dominated results described in the previous section. Results are plotted separately by rotation speed with values of (a) 0.4 rad s$^{-1}$, (b) 1.0 rad s$^{-1}$, (c) 2.0 rad s$^{-1}$ and (d) 4.0 rad s$^{-1}$. Clearly, deviations from the diffusion-dominated limit increase with gap height and angular rotation speed due to the increase of inertial secondary flows, as well as with increasing $c_{sat}$ due to local reductions in viscosity (and consequent increases in inertial effects).

Comparing these results with the experimental results, the axisymmetric inertial simulations do seem to capture many features of the experimental results. In particular, we generally see decreased $\mu _f$ values at larger gap heights and larger $c_{sat}$ values (i.e. RH values), which are consistent with figure 2. Specifically, the numerical results in figure 12(a) correspond to the same rotation speed ($\varOmega =0.4$ rad s$^{-1}$) as figure 2(b). Similar trends are seen (except at small gap heights, which is discussed below), with slightly less significant decreases in viscosity in the simulations compared with the experiments. Increasing the rotation speed to 1.0 rad s$^{-1}$ in the simulations shows more significant viscosity decreases than the experimental results at 0.4 rad s$^{-1}$. So the experimental results at 0.4 rad s$^{-1}$ agree well quantitatively with numerical predictions slightly above 0.4 rad s$^{-1}$. One trend seen in the experiments that we do not see in the simulations is the non-monotonic relationship between viscosity decrease and angular rotation rate seen in figure 2(d). However, we do see evidence of a non-monotonic relationship with increasing inertial effects in the simulations. In particular, at $\varOmega =4.0$ rad s$^{-1}$ with $c_{sat}=0.1$ (orange curve in figure 12d), the final measured viscosity first decreases and then increases with increasing gap height, demonstrating that the axisymmetric case can demonstrate such trends.

The most significant experimental result that these simulations cannot explain is the large decrease in measured viscosity values at small gap heights. In fact, one of the most consistent results of the axisymmetric, inertial simulations is the approach to the 1-D, diffusion-dominated regime at small gap heights for any angular rotation speed. Thus, these axisymmetric simulations apparently fail to account for some effect that becomes dominant at small gap heights. We hypothesize that this is due to misalignment effects that only become significant at very small gap heights in practical parallel-plate rheometers. In the next section, we perform additional simulations based on a misaligned geometry in an attempt to validate this hypothesis.

7. Role of misalignment

In the previous section, we clearly saw that the axisymmetric, inertial, variable-viscosity model fails to account for the sharp decrease in measured viscosity at small gap heights. Thus, we must consider what possible sources of error could account for these effects. A variety of experimental challenges exist for performing accurate measurements with a rheometer, such as underfilling of the parallel-plate gap, instrument inertia and surface tension effects (Hellström et al. Reference Hellström, Samaha, Wang, Smits and Hultmark2014; Ewoldt et al. Reference Ewoldt, Johnston and Caretta2015). In addition to these, there are practical sources of error associated with the mechanical uncertainties in the rheometer itself. A key source of these errors comes from deviations in the geometry of the gap containing the fluid. These errors in the gap geometry could arise from non-parallelism, non-concentricity, non-flatness of the plates, non-zero slip lengths at the upper or lower plates, edge effects at the outer edge of the rheometer or errors in the gap-zeroing procedure (Connelly & Greener Reference Connelly and Greener1985; Kramer, Uhl & Prud'Homme Reference Kramer, Uhl and Prud'Homme1987; Kalika, Nuel & Denn Reference Kalika, Nuel and Denn1989; Henson & Mackay Reference Henson and Mackay1995; Davies & Stokes Reference Davies and Stokes2005; Andablo-Reyes, Hidalgo-Álvarez & de Vicente Reference Andablo-Reyes, Hidalgo-Álvarez and de Vicente2010; Andablo-Reyes, de Vicente & Hidalgo-Alvarez Reference Andablo-Reyes, de Vicente and Hidalgo-Alvarez2011). One reason the discussion of these sources of error arose was due to the experimental observation that as gaps decreased below several hundred micrometres, measured viscosities began to have systematic errors, typically decreasing with the gap height as also shown in figure 2(b) (Walters Reference Walters1975; Pipe & McKinley Reference Pipe and McKinley2009).

Based on these observations, a variety of studies suggest that a key factor in this discrepancy in our measurements and simulations could be the misalignment of the rotating plate. Although it is commonly assumed that the plates are perfectly aligned, a number of reports indicate that a small but finite misalignment is prevalent in parallel-plate and cone-and-plate rheometers (Andablo-Reyes et al. Reference Andablo-Reyes, de Vicente and Hidalgo-Alvarez2011). In fact, the gap height can vary over 50 $\mathrm {\mu }$m across a few centimetres in a parallel-plate rheometer due to the non-parallelism in the gap, causing a significant error in the viscosity measurements in narrow-gap, high-shear-rate experiments (Davies & Stokes Reference Davies and Stokes2008; Pipe, Majmudar & McKinley Reference Pipe, Majmudar and McKinley2008). This is due to the fact that the misalignment introduces additional lubrication forces in the fluid layer. A variety of semi-empirical techniques have been developed to account for these systematic errors at small gap heights. For example, a simple linear approximation has been proposed in which a simple gap error is defined to correct the measured values (Connelly & Greener Reference Connelly and Greener1985; Dudgeon & Wedgewood Reference Dudgeon and Wedgewood1994; Davies & Stokes Reference Davies and Stokes2008). Another technique involves using ultrasound time-of-flight measurements to detect the varying thickness of the fluid layer in the case of misalignment, which can be used to calculate the degree of misalignment (Rodríguez-López et al. Reference Rodríguez-López, Elvira, de Espinosa Freijo and de Vicente2013). A numerical solution of the flow in a misaligned parallel-plate rheometer was also presented by Andablo-Reyes et al. (Reference Andablo-Reyes, Hidalgo-Álvarez and de Vicente2010). Also, Clasen (Reference Clasen2013) introduced a system that can self-correct non-parallelism to a degree using hydrodynamic lubrication forces. Finally, a theoretical description of the velocity and stress profiles in a slightly misaligned cone-and-plate rheometer was achieved by Dudgeon & Wedgewood (Reference Dudgeon and Wedgewood1994) using a domain perturbation study in the limit of zero Reynolds number.

In this section, we consider the role of misalignment in the transport of absorbed water throughout the glycerol, and the effects of this misalignment on the measured viscosity values. In general, such an analysis would need fully three-dimensional simulations in a misaligned rheometer geometry. We considered performing such simulations, but found them to be intractable due to the extremely high computational cost of performing them. In particular, they require 2–3 orders of magnitude more grid cells than the axisymmetric simulations in order to resolve the water concentration boundary layer at the outer edge of the rheometer. Furthermore, all of the meshing techniques we tried that would maintain this resolution at the outer edge ultimately resulted in very high-aspect-ratio cells at some point in the domain that affect resolution and greatly increase the number of iterations needed to solve the velocity–pressure profile with the SIMPLE algorithm, which must be repeated continuously as the water concentration field evolves. For these reasons, we consider a simplified, depth-averaged case that is valid in the limit of small gap heights. This model and the corresponding simulations and results are described in the following sections.

7.1. Theory

Here, we consider a misaligned parallel-plate rheometer with radius $R$ and gap thickness $h(r,\theta,\phi )$, where the upper plate is slightly tilted by a small angle $\phi$. Once again, the lower plate is stationary, and the upper plate rotates at a rotation speed of $\varOmega$. The coordinate system and problem set-up are shown in figure 13. The boundary conditions for the system are $\boldsymbol {u}=\boldsymbol {0}$ at $z=0$ (on the lower plate), and

(7.1a)\begin{gather} u_r=r\cos\theta\sin\theta\sin\phi\tan\phi, \end{gather}

(7.1b)\begin{gather}u_\theta=r(\cos^2\theta\sec\phi+\cos\phi\sin^2\theta), \end{gather}

(7.1c)\begin{gather}u_z={-}r\epsilon^{{-}1}\sin\theta\sin\phi, \end{gather}

at $z= h(r,\theta,\phi )=1+\phi \epsilon ^{-1}r\cos \theta$ (the upper plate). Note that these simplify to $u_r=u_z=0$ and $u_\theta =r$ at $z=1$ as in the axisymmetric case when $\phi =0$. For small angles, the angle $\phi$ can range from 0 to a maximum of $\epsilon$. Thus, $\phi /\epsilon$ ranges from 0 to 1, and small values of $\phi /\epsilon$ correspond to small plate deflections. Here, $\phi /\epsilon =0$ corresponds to the case of no misalignment and $\phi /\epsilon =1$ corresponds to the case where the plates come in contact at one edge. Further, recall that as before we generally also need to apply boundary conditions at $r=1$. In a practical experiment, this boundary condition represents a fluid–air interface that is typically not flat and experiences surface tension effects. However, in the small-gap limit, we lose the ability to impose such a boundary condition, and we note that this contributes to the error in velocity–pressure profiles in the ${O}(\epsilon )$ region near $r=1$.

Figure 13.

Misaligned parallel-plate rheometer geometry and coordinate system. The upper plate is misaligned by a small deflection angle $\phi$ and rotates at angular speed $\varOmega$. With $z$ non-dimensionalized by $h_0$ and $r$ non-dimensionalized by $R$, the $z$ coordinate defining the upper plate is $h(r,\theta,\phi )=1+\phi \epsilon ^{-1}r\cos \theta$. Note that for small angles, the angle $\phi$ can range from 0 to a maximum of $\epsilon$.

The governing equations are again the Navier–Stokes equations with variable viscosity and the continuity equation, which are given by (3.3) and (3.4), respectively, as well as the water concentration advection–diffusion equation given by (3.6). As mentioned above, the numerical simulation of the full system of coupled equations in a well-resolved three-dimensional geometry is computationally expensive. In the limit of narrow gap heights $\epsilon \ll 1$ and negligible inertia ${Re}\ll 1$, the Navier–Stokes equations simplify to

(7.2a)\begin{gather} 0={-}\frac{\partial p}{\partial r} +\mu\frac{\partial^2u_r}{\partial {z}^2} + \frac{\partial \mu}{\partial z}\frac{\partial u_r}{\partial z}, \end{gather}

(7.2b)\begin{gather}0={-}\frac{1}{r}\frac{\partial p}{\partial\theta}+\mu\frac{\partial^2u_\theta}{\partial {z}^2} + \frac{\partial \mu}{\partial z}\frac{\partial u_\theta}{\partial z}, \end{gather}

(7.2c)\begin{gather}0=\frac{\partial p}{\partial z}. \end{gather}

Furthermore, in the thin-gap limit, the water concentration can be assumed to be approximately uniform in the depth direction, which gives ${\partial c}/{\partial z}\approx 0$ and ${\partial \mu }/{\partial z}\approx 0$. This gives

(7.3a–d)\begin{equation} \left.\begin{gathered} 0={-}\frac{\partial p}{\partial r} +\mu\frac{\partial^2u_r}{\partial {z}^2},\quad 0={-}\frac{1}{r}\frac{\partial p}{\partial\theta}+\mu\frac{\partial^2u_\theta}{\partial {z}^2},\\ 0=\frac{\partial p}{\partial z}\quad\text{and}\quad\frac{1}{r}\frac{\partial}{\partial r}\left(r u_r\right) + \frac{1}{r}\frac{\partial u_\theta}{\partial \theta} + \frac{\partial u_z}{\partial z} = 0. \end{gathered}\right\} \end{equation}

Note here that $c(r,\theta,t)$, $\mu (r,\theta,t)$ and $p(r,\theta,t)$ in this limit. With the fact that $\partial p/\partial z=0$ in this limit, the next leading-order form of the $z$ component of the Navier–Stokes equations becomes

(7.4)\begin{equation} 0=\mu\frac{\partial^2 u_z}{\partial z^2}+\frac{\partial\mu}{\partial r} \frac{\partial u_r}{\partial z} + \frac{1}{r}\frac{\partial \mu}{\partial\theta} \frac{\partial u_\theta}{\partial z}. \end{equation}

By examining the form of the gap height distribution $h(r,\theta,\phi )=1+\phi \epsilon ^{-1}r\cos \theta$ we see that the magnitude of the perturbation is ${O}(\phi /\epsilon )$. This suggests the use of a solution given by

(7.5a)\begin{gather} u_r(r,\theta,z,t)=\left(\frac{\phi}{\epsilon}\right)u_{r,1}(r,\theta,z,t)+ \left(\frac{\phi}{\epsilon}\right)^2u_{r,2}(r,\theta,z,t)+\cdots, \end{gather}

(7.5b)\begin{gather}u_\theta(r,\theta,z,t)=rz+\left(\frac{\phi}{\epsilon}\right)u_{\theta,1}(r,\theta,z,t)+ \left(\frac{\phi}{\epsilon}\right)^2u_{\theta,2}(r,\theta,z,t)+\cdots, \end{gather}

(7.5c)\begin{gather}u_z(r,\theta,z,t)=\left(\frac{\phi}{\epsilon}\right)u_{z,1}(r,\theta,z,t)+ \left(\frac{\phi}{\epsilon}\right)^2u_{z,2}(r,\theta,z,t)+\cdots, \end{gather}

(7.5d)\begin{gather}p(r,\theta,t)=\left(\frac{\phi}{\epsilon}\right)p_1(r,\theta,t)+ \left(\frac{\phi}{\epsilon}\right)^2p_2(r,\theta,t)+\cdots, \end{gather}

which is valid in the limit $\phi /\epsilon \ll 1$. Substituting this expansion into (7.3a–d) and (7.4) and applying the boundary conditions gives

(7.6a)\begin{gather} u_{r,1}(r,\theta,z,t)=\frac{z(z-1)}{2\mu}\frac{\partial p_1}{\partial r}, \end{gather}

(7.6b)\begin{gather}u_{\theta,1}(r,\theta,z,t)={-}r^2 z \cos\theta + \frac{z(z-1)}{2 r \mu}\frac{\partial p_1}{\partial\theta}, \end{gather}

(7.6c)\begin{align} u_{z,1}(r,\theta,z,t)&={-}\frac{z^2}{12r^2\mu}\left\{6r^3\mu^2\sin\theta - (2z-3) \left[\frac{\partial\mu}{\partial\theta}\frac{\partial p_1}{\partial\theta} + r^2 \frac{\partial\mu}{\partial r}\frac{\partial p_1}{\partial r}\right]\right. \nonumber\\ &\quad \left. +\, (2z-3)\mu\left[\frac{\partial^2 p_1}{\partial \theta^2}+r\left(\frac{\partial p_1}{\partial r}+r\frac{\partial^2 p_1}{\partial r^2}\right)\right]\right\}. \end{align}

This procedure also yields a partial differential equation governing the pressure distribution that is given by

(7.7)\begin{equation} \mu\left\{\frac{\partial^2 p_1}{\partial\theta^2} + r\left[\frac{\partial p_1}{\partial r} + r\left(6r\mu\sin\theta + \frac{\partial^2 p_1}{\partial r^2}\right)\right]\right\} = \frac{\partial\mu}{\partial\theta}\frac{\partial p_1}{\partial\theta} + r^2 \frac{\partial\mu}{\partial r}\frac{\partial p_1}{\partial r}. \end{equation}

With some known distribution of viscosity in the system, a numerical solution of (7.7) yields the pressure distribution in the gap. This in turn can be used to calculate the velocity profiles from (7.5) and (7.6). The velocity profiles can then be used to update the water concentration distribution via the advection–diffusion equation. With the assumption that $c$ is independent of $z$ (valid in the small-gap limit), the solute transport equation becomes

(7.8)\begin{equation} Pe\frac{\partial c}{\partial t}=\frac{1}{r}\frac{\partial}{\partial r} \left(rD\frac{\partial c}{\partial r}\right)+\frac{1}{r^2}\frac{\partial}{\partial\theta} \left(D\frac{\partial c}{\partial\theta}\right)-Pe\left(u_r\frac{\partial c}{\partial r}+ \frac{u_\theta}{r}\frac{\partial c}{\partial \theta}\right). \end{equation}

Since $c$ is independent of $z$, we consider solving the depth-averaged version of this equation instead, which is simply

(7.9)\begin{equation} Pe\frac{\partial c}{\partial t}=\frac{1}{r}\frac{\partial}{\partial r} \left(rD\frac{\partial c}{\partial r}\right)+\frac{1}{r^2}\frac{\partial}{\partial\theta} \left(D\frac{\partial c}{\partial\theta}\right)-Pe\left(\bar u_r\frac{\partial c}{\partial r}+ \frac{\bar u_\theta}{r}\frac{\partial c}{\partial \theta}\right), \end{equation}

where bars denote depth-averaged quantities, $\bar c = c$, $\bar D = D$ and

(7.10a)\begin{gather} \bar u_r = \frac{1}{1+\dfrac{\phi}{\epsilon}r\cos\theta} \int_{z=0}^{1+({\phi}/{\epsilon})r\cos\theta} u_r(r,\theta,z)\,\mathrm{d}z ={-}\frac{1}{12\mu}\frac{\phi}{\epsilon}\frac{\partial p_1}{\partial r} + {O}\left(\frac{\phi}{\epsilon}\right)^2, \end{gather}

(7.10b)\begin{gather}\bar u_\theta = \frac{1}{1+\dfrac{\phi}{\epsilon}r\cos\theta} \int_{z=0}^{1+({\phi}/{\epsilon})r\cos\theta} u_\theta(r,\theta,z)\,\mathrm{d}z = \frac{r}{2} - \frac{1}{12 r \mu}\frac{\phi}{\epsilon}\frac{\partial p_1}{\partial\theta} + {O}\left(\frac{\phi}{\epsilon}\right)^2. \end{gather}

Thus, having solved the pressure distribution due to the misalignment $p_1$ from (7.7), the depth-averaged velocity components can be calculated from (7.10), which can in turn be used to advect the solute concentration.

7.2. Numerical methods

We perform numerical simulations of the coupled transport equations described in the previous section. Recall that we seek solutions of the water transport and associated viscosity measurements in a misaligned parallel-plate rheometer that are valid at small gap heights and in the limit of negligible inertia. The numerical approach for solving these systems is as follows:

1. (i) Solve (7.7) for the pressure perturbation due to misalignment subject to the boundary condition $p\rightarrow 0$ at $r=1$.

2. (ii) Calculate the depth-averaged radial and azimuthal velocity components from (7.10).

3. (iii) Advance the water concentration profile in time by numerically integrating (7.9) for one or more time steps.

4. (iv) Calculate the new viscosity and diffusivity fields and iterate back to step (i).

We attempt a numerical implementation of this process using a finite-difference implementation in MATLAB. However, several complications emerge due to the extremely large Péclet numbers in the system. In particular, the diffusion of the water concentration field is so slow that the field effectively propagates with a very sharp front. In order to resolve this and avoid spurious oscillations in the concentration field, we use slope-limited finite differencing based on the minmod limiter function to switch to first-order spatial differencing at the steep gradient (Roe Reference Roe1986). This avoids the oscillations that result in a pure second-order differencing scheme and still allows nearly second-order accuracy in space globally. The more serious difficulty that we encountered with this numerical approach is solving (7.7) for the perturbation pressure field. These solutions do not behave nicely due to the sharp viscosity gradients on the right-hand side of the equation. We were not able to resolve this issue using slope limiters.

As an alternative approach and to illustrate the qualitative dynamics that can be expected with a misaligned upper plate, we instead assume a constant-viscosity model for the purposes of calculating the velocity profile, since this solution is well behaved. Note that (7.7) has an analytical solution when $\mu =1$ which is given by

(7.11)\begin{equation} p_1(r,\theta) ={-}\tfrac{3}{4}r\left({-}1+r^2\right)\sin\theta. \end{equation}

We use this theoretical result at constant viscosity to calculate the depth-averaged velocities in the misaligned rheometer and use these to update the concentration profile. When calculating the effective measured viscosity and torque, we always use the viscosity distribution that corresponds to the concentration profile. This limitation to a velocity profile based on constant viscosity is a clear limitation of our results, but nevertheless they capture qualitative features of the experiments that the axisymmetric model could not predict, and we leave a full solution with evolving velocity profiles based on spatial variation of viscosity to future work.

7.3. Results

Here, we introduce the numerical results achieved for the misaligned rheometer system based on the methodology described in the previous section. First, we highlight the role of the misalignment in the water concentration field in figure 14. Here, the depth-averaged water concentration profiles are shown at $t^*=3600$ s as a function of misalignment for $c_{sat}=0.5$ with misalignment ranging from $\phi /\epsilon =0$ to 0.95, where 0 is the perfectly aligned parallel-plate case, and 1 is the limit where the plates come into contact at one edge. As can be seen, the misaligned cases all show a non-axisymmetric concentration profile. This is due to the non-axisymmetric secondary velocity components due to the misalignment. In particular, the radial component $\bar u_r$ transports water towards or away from the outer edge, and is $\bar u_r \sim ({\phi }/{\epsilon })({\partial p_1}/{\partial r})$. With $p_1\sim \sin \theta$, this represents a radially inward flow on one half of the rheometer and a radially outward flow on the other half and explains why in figure 14 the concentration profile appears to be pulled in from the right edge and pushed towards the left edge. Furthermore, this secondary velocity component is proportional to $\phi /\epsilon$, so doubling the degree of misalignment doubles the radial advective fluxes in both directions.

Figure 14.

Depth-averaged water concentration profile at $t^*=3600$ s as a function of misalignment for $c_{sat}=0.5$ and $\varOmega =0.4$ rad s$^{-1}$. Results correspond to $\phi /\epsilon$ values of (a) 0 (perfectly aligned), (b) 0.1, (c) 0.2, (d) 0.5, (e) 0.75 and (f) 0.95 (plates are nearly contacting). As can be seen, as the misalignment increases, the concentration profile becomes no longer axisymmetric, and there is a significant increase in total water transport into the fluid layer from the outer edge due to the misalignment-driven secondary flows.

In order to quantify the effect of misalignment on the measured viscosity values, simulations were performed across a range of $\phi /\epsilon$ and $c_{sat}$ values. The final measured viscosity values from these simulations are presented in figure 15. Here, the results correspond to an angular speed of 0.4 rad s$^{-1}$. The dashed lines in the figure correspond to the 1-D diffusion-dominated regime, and the results asymptotically approach these limits as $\phi /\epsilon \rightarrow 0$. Furthermore, the results show a steep drop-off in measured final viscosities at large misalignments, which possibly explains the sharp decrease in measured viscosity in the experiments at small gap heights (e.g. figure 2b) that was not captured in the axisymmetric model (e.g. figure 12).

Figure 15.

Final measured viscosity values $\mu _f/\mu _i$ at $t^*=3600$ s as functions of the misalignment $\phi /\epsilon$ and $c_{sat}$ at an angular rotation speed of 0.4 rad s$^{-1}$. Dashed lines correspond to the 1-D pure diffusion limit. Results asymptotically approach the 1-D diffusion-dominated limit as $\phi /\epsilon \rightarrow 0$. Furthermore, results show a steep drop-off in measured viscosity values as the misalignment increases.

Finally, a comparison between all of the three different proposed models (i.e. the 1-D pure diffusion limit, the axisymmetric inertial limit and the misaligned inertialess small-gap limit) is shown in figure 16. Here, the dashed lines indicate the 1-D diffusion-dominated limit, dot-dashed lines correspond to the inertial, axisymmetric regime and the solid lines show the results for the misaligned, inertialess, small-gap limit. Results were calculated for angular rotation speeds of 0.4 and 1.0 rad s$^{-1}$. The results show that at small gap heights, the misalignment effects dominate the inertial effects. This becomes more clear when considering that the radial secondary velocity due to inertial effects in an axisymmetric case is ${O}({Re})$, whereas the secondary velocity components due to misalignment are ${O}(\phi /\epsilon )$. For fixed angular rotation speed, the secondary velocities must dominate the inertial secondary velocities as the gap height decreases. Furthermore, the results also show that the opposite is true at large gap heights. For a fixed misalignment angle and rotation speed, $\phi /\epsilon$ decreases as the gap height increases, whereas the inertial effects increase, such that the large-gap regime is dominated by inertial effects. This cross-over explains the non-monotonic relationship between $\mu _f$ and gap height reported in figure 2(b). Thus, there is a critical $\epsilon$ value at which the measured viscosity values switch from being misalignment-dominated to inertia-dominated. Numerical simulations using the previously described models and numerical methods can be used to estimate this transition, as shown in figure 16.

Figure 16.

Comparison of the final measured viscosity values at $t^*=3600$ s for the models corresponding to each of the three regimes: (1) the 1-D axisymmetric, pure diffusion limit (dashed lines), (2) the axisymmetric, inertial regime (dot-dashed lines) and (3) the misaligned, inertialess, small-gap limit (solid lines). Results correspond to angular rotation speeds of (a) 0.4 rad s$^{-1}$ and (b) 1.0 rad s$^{-1}$. Misaligned cases are calculated with a misalignment angle of 0.0005 rad. All cases were performed with a rheometer of radius $R=2.5$ cm.

8. Comparison between experiments and simulations

Here, we provide a quantitative comparison between the experimental and numerical results. Figure 17 shows a detailed comparison between both the transient viscosity measurements and the final measured viscosities at $t^*=3600$ s. First, figure 17(a) shows the transient evolution of the measured and predicted viscosities normalized by the initial viscosity at $t=0$. Solid lines correspond to the experimental measurements and are measured at RHs of 5 %, 45 % and 72 %. Dashed lines correspond to the numerical results based on the misalignment-dominated model and are calculated at values of $c_{sat}=0.05$, 0.2 and 0.4, respectively. These results are calculated at a gap height of 0.1 mm and rotation speed of 0.4 rad s$^{-1}$. The misalignment parameter is $\phi /\epsilon =0.3$, which corresponds to a misalignment tilt angle of approximately 0.0012 rad, which is slightly higher than typical estimates. As can be seen, the misaligned modelling is able to roughly capture the qualitative features of the experimental measurements, as well as a rough order of magnitude of the effects on the transient viscosities.

Figure 17.

Comparison between experimental results and numerical simulations. (a) Transient viscosity measurements normalized by the initial viscosity at RHs of 5 % (black), 45 % (grey) and 72 % (blue). Solid lines correspond to experiments and dashed lines correspond to simulations. Numerical simulations are calculated via the misalignment-dominated modelling and use $c_{sat}=0.05$, 0.2 and 0.4, respectively, for the black, grey and blue curves. Here, the gap height is 0.1 mm and the rotation rate is 0.4 rad s$^{-1}$. (b) Measured final viscosities $\mu _f$ of glycerol at $t^*=3600$ s normalized by $\mu _i$ with $\varOmega =0.4$ rad s$^{-1}$ and $R=2.5$ cm. Symbols correspond to experimental results at $\textrm {RH}=72\,\%$, and the solid and dash-dotted blue lines correspond to predictions from the misalignment-dominated and the inertia-dominated modelling, respectively, with $c_{sat}=0.5$.

There are several likely sources of error in this modelling that could be improved to achieve a better fit. First, the numerical results show a rapid decrease in the measured viscosity at early times. This is likely due to the fact that the lubrication-style depth-averaged modelling is not able to satisfy the boundary conditions at the outer edge of the rheometer. Instead, the modelling predicts a relatively strong radially inward flow at the outer boundary that rapidly pulls in absorbed water from the saturated boundary condition. This effect would likely be lessened with a more comprehensive model of the outer interface. Next, while the numerical predictions all begin with exactly pure glycerol, the initial measured viscosity in the experiments actually decreases some with increasing RH. This is due to the time it takes to prepare the experiment, in which the glycerol is able to absorb some water from the air. The consequence of this is that the viscosity has already decreased to a small degree in the experiments, which may contribute to relatively stronger fluid advection and enhanced water transport. Finally, our modelling assumes the dynamics are in either the inertia-dominated regime or misalignment-dominated regime. In fact, there is a transition region where both inertia and misalignment contribute to enhanced water transport relative to the pure diffusion limit. In this transition region, both models will underpredict the rate of water absorption, and thus predict slower-than-expected rates of decrease for the viscosity.

Next, figure 17(b) shows a comparison between the experimentally measured final viscosities (symbols) at $t^*=3600$ s and the predictions for both the misalignment-dominated (solid blue line) and inertia-dominated (dot-dashed line) modelling. Here, $\varOmega =0.4$ rad s$^{-1}$, $R=2.5$ cm, $c_{sat}=0.5$ in the modelling, and $\textrm {RH}=72\,\%$ in the experiments. The misalignment angle is assumed to be 0.0005 rad. As can be seen, the experimental results demonstrate a smooth transition between the two regimes, and together both models capture the qualitative shape of the data, as well as the rough order of magnitude of the effects.