Yield stress fluids solidifying in capillary imbibition

Abstract

When subjected to external stress exceeding their yield stress ($\sigma _{Y}$), yield stress fluids (YSFs) undergo transition from solid to liquid. While the flow dynamics of capillary flows – including wetting, coating, spreading, and wicking – can be predicted accurately for fluids without yield stress, only a few asymptotic solutions have been validated against experiments for YSFs. In this study, we explore experimentally and theoretically the entire dynamics and solidifying process during capillary imbibition of YSFs, representing the fundamental aspect of capillarity-driven flow. We observe that during the gradual deceleration of capillary rise, YSFs exhibit apparent flow below $\sigma _{Y}$. This unexpected behaviour is linked to wall slip, manifested as plug flow in particle image velocimetry analyses. We introduce a model with numerical solutions, grounded in the rheological properties of YSFs and slip mechanisms, to accurately describe this liquid–solid transitional behaviour.

1. Introduction

A yield stress fluid (YSF) keeps its shape like a solid until a sufficiently large stress ($>\sigma _{Y}$) initiates a liquid-like flow (Coussot et al. 2002; Bonn et al. 2017). Such YSFs are widely found in biofluids (cell medium and mucus) (Hu et al. 2015), as well as in various industries, such as paintings, cosmetics and foods. For YSFs, the solid-to-liquid change is the key flow behaviour commonly denoted as yielding (Møller, Fall & Bonn 2009; Balmforth, Frigaard & Ovarlez 2014; Malkin, Kulichikhin & Ilyin 2017; Coussot 2018; McConnell et al. 2022; Dey et al. 2023) or fluidization (Bocquet, Colin & Ajdari 2009; Benzi et al. 2019, 2021; Garg et al. 2021; Jadhav, Rossi & Karimfazli 2021). The yielding depends on the strain rate and stress, and shows a gradual transition rather than a sudden event (Oldroyd 1947; Bonn & Denn 2009; Coussot & Rogers 2021). Since plastic deformation takes place after elastic response, stress overshoot appears as a consequence of microstructural rearrangement, and transient/permanent shear banding follows due to its inhomogeneity (Benzi et al. 2021).

Contrary to yielding, which is a transition from solid to liquid under stress exceeding $\sigma _{Y}$, solidification in YSFs represents a liquid-to-solid transition as stress decreases below $\sigma _{Y}$. The disparity between dynamic and static yield stresses suggests that the solidification process differs fundamentally, requiring a distinct analytical approach (Divoux, Barentin & Manneville 2011; Chaudhuri, Berthier & Bocquet 2012; Balmforth et al. 2014). Solidification processes have been studied across diverse contexts, ranging from time-dependent simple shear scenarios to the cessation of motion in particles, bubbles, and droplets (Huilgol, Mena & Piau 2002; Matson & Hogg 2007; Tokpavi, Magnin & Jay 2008; Tsamopoulos et al. 2008; Sikorski, Tabuteau & de Bruyn 2009; Dimakopoulos, Pavlidis & Tsamopoulos 2013; Damianou et al. 2014; Jalaal et al. 2018; Huilgol, Alexandrou & Georgiou 2019). It has been shown that the criteria for cessation are influenced by factors including yield stress, stress during solidification, surface tension, and relaxation time. Given its prevalence in capillarity-driven flows like wetting, coating, spreading and wicking (Huilgol & Kefayati 2018), solidification warrants specific attention within the scope of capillary imbibition. While a few models address the flow dynamics of YSFs during solidification (Bertola 2009; Géraud et al. 2014; Huilgol & Kefayati 2018; Balmforth 2019) and compare them with experimental results (Balmforth, Ghadge & Myers 2007; Jalaal, Stoeber & Balmforth 2021; van der Kolk, Tieman & Jalaal 2023), mostly they focus on the asymptotic solution for the final state ($t \rightarrow \infty$). Notably, none has been compared directly with the experimental dynamics throughout the entire wetting process.

In this work, our primary focus is the solidification of YSFs in capillary rise (figure 1a). From the conservation of momentums acting on the liquid (density $\rho$), the wall shear stress on a capillary inner surface with a sufficiently narrow radius $R$ is a function of rise height ($h(t)$) as expressed by

(1.1)\begin{equation} \sigma_{w}(t)=\frac{{\varGamma}\cos \theta_{E}}{h(t)}-\frac{\rho gR}{2}, \end{equation}

with the constitutive equation for Newtonian and shear thinning fluids,

(1.2)\begin{equation} \dot{h}(t)=\frac{R}{s+2}\left(\frac{\sigma_{w}(t)}{K}\right)^{s-1}, \end{equation}

where ${\varGamma }$ is the surface tension of the liquid–air interface, $\theta _{E}$ is the equilibrium contact angle between the liquid and the capillary, $g$ is the gravitational constant, $K$ is the consistency parameter, $n$ is the thinning power, and $s=1/n+1$ (Lucas 1918; Washburn 1921; Rideal 1922). Equation (1.2) can be generalized for the Newtonian fluid by putting $n=1$ and $K$ as the viscosity. Because $h(t)$ increases in time, $\sigma _{w}(t)$ finally converges to zero for typical constitutive equations (Newtonian or shear thinning fluid). This corresponds to a monotonic rise curve of aqueous glycerol in figure 1(b). The YSFs, however, undergo solidification at $h_{p}={\varGamma }\cos \theta _{E}/(\sigma _{Y}+\rho gR/2)$, the height where $\sigma _{w}=\sigma _{Y}$ (Balmforth 2022). We found that the solidification slows down YSFs at $h_{p}$, but is followed by another onset of capillary rise, as in figure 1(c). This contradicts a previous study that says a YSF stops at $h_p$ (Géraud et al. 2014). Instead, slip attributes the flow after the solidification. Our proposed model predicts the complete capillary rise of YSFs quite accurately, and explains what determines the flow dynamics after solidification based on the relevant length and time scales.

Figure 1. (a) Schematic for the capillary rise experiment. Capillary rise $h(t)$ of (b) 85 $\%$ glycerol (Newtonian fluid, $\sigma _{Y}=0$) and (c) 0.12 wt. $\%$ Carbopol (YSF, $\sigma _{Y}=10$ Pa). Dashed lines in (b,c) are (1.1) combined with (1.2) for Newtonian and simple power law fluids ($\sigma _{Y}=0$), respectively. The solid line is (1.1) with (2.1). The shading is an error bar. Experimental parameters follow table 1.

2. Results and discussion

Figure 1(b) shows that Newtonian fluids show a monotonic increase until reaching equilibrium height ($h_{E}=2{\varGamma } \cos \theta _{E}/\rho gR$), where Laplace pressure is balanced by gravity (see Appendix A for experimental details). The model curve aligns with Newtonian fluid data when Hagen–Poiseuille flow is assumed. However, in figure 1(c), the model fails to match for YSFs with a simple power-law fluid, possibly due to the presence of yield stress and inability to account for yielding or solidification impacts.

Accordingly, a well-known Herschel–Bulkley (HB) model is implemented to replace (1.2). For $\sigma >\sigma _{Y}, \sigma =\sigma _{Y}+K\dot {\gamma }^{n}$, where $\dot {\gamma }$ is the shear rate. When $\sigma \le \sigma _{Y}$, $\dot {\gamma }=0$. The apparent rise velocity $\dot {h}(t)$ is obtained by averaging the shear rate in radial direction: $\dot {h}(t)=\int _{0}^{R}2{\rm \pi} r\,\dot {\gamma }(r,t)\,{\rm d} r$ (Sochi & Blunt 2008; Pérez-González et al. 2012). By replacing the shear rate with the HB model, it rearranges and explicitly links $\dot {h}(t)$ to the wall shear stress ($\sigma _{w}(t)$) as

(2.1)\begin{equation} \dot{h}_{HB}(t)=R\left(\frac{\sigma_{w}(t)}{K} \right)^{s-1}f(\epsilon(t),s), \end{equation}

where

(2.2)\begin{equation} f=(1-\epsilon(t))^s\left\{\frac{(1-\epsilon(t))^2}{s+2} +\frac{2\epsilon(t)\,(1-\epsilon(t))}{s+1}+\frac{\epsilon(t)^2}{s}\right\}, \quad \epsilon(t)=\frac{\sigma_{Y}}{\sigma_{w} (t)}<1. \end{equation}

For regime $\epsilon (t)\ge 1$, $\dot {h}=0$. Upon applying $\epsilon =0$, (2.1) reduces to (1.2), where $f$ becomes ${1}/({2+s})$.

With the implementation of (2.1), the numerical solution with measured $\sigma _{Y}$ is shown as a solid line in figure 1(c) (see Appendix B for the rheology measurement result). The plateau height, denoted as $h_{p}$, was determined precisely by the model. This indicates that the plateau is a result of the solidification process of the YSF. For all the theoretical curves, dynamic contact angle (DCA) theory is employed to account for the wetting friction of the meniscus (Kim et al. 2020). The DCA description and its fitting are available in Appendix C and figures 5(a,b) below. The YSF surface tension is estimated from the Laplace pressure at $h_{E} = h(t\rightarrow \infty )$ due to the lack of reproducibility of the previously proposed method, which heavily relies on the selection of geometry (Jørgensen et al. 2015; Jalaal et al. 2021). The measured surface tension is given in table 1, where its value varies by polymer concentration. These values can be considered as apparent surface tensions acting on a capillary imbibition system, since the value can be hindered dynamically by the yield stress (Martouzet et al. 2021).

Table 1. Rheological properties and slip layer thicknesses: $\sigma _{HB}$ is fitted from a rheology measure, while $\sigma _{Y}$ is fitted from capillary rise $h(t)$. Theoretical prediction follows $\sigma _{Y}$.

Interestingly, stress overshoot has negligible impact on the macroscale dynamics of capillary rise. Recent experimental and theoretical studies have demonstrated that the transition from solid to liquid is accompanied by a stress overshoot, as evidence of plastic deformation of YSF (Donley et al. 2020; Kamani, Donley & Rogers 2021). In a case of the capillary rise, the good fit of (2.1) to the initial rising speed indicates that stress overshoot has little effect on the macroscale flow dynamics; with a sufficiently large Laplace stress ($\sigma \sim 100$ Pa) at the beginning, the characteristic time scale for yielding rate (${\sim }10^{-2}$ s) is an order of magnitude faster than the time scale of our measurements (${\ge }10^{-1}$ s) (Benzi et al. 2021).

As shown in figure 1(c), the experimental curve of YSFs shows a significant deviation from the saturation predicted by (2.1) (solid line) and continues to rise until $h_{E}$. Derived from (1.1), which is relevant to all instances of capillary rise irrespective of the fluid's rheological characteristics, it is observed that the wall shear stress diminishes with an increase in penetration length. Hence the rise beyond $h_{p}$ is indicative of a decline in wall shear stress. Ultimately, this anomalous flow behaviour is associated with the culmination of the liquid-to-solid transition, during which shear stress continues to decrease even after surpassing $\sigma _{Y}$. This makes solidification distinct from the re-entrance effect under long-term creep tests (Landrum, Russel & Zia 2016) and thixotropic aging (Choi, Armstrong & Rogers 2021). One possibility is that YSFs may not undergo complete solidification, as seen in the analogy of delayed yielding (Landrum et al. 2016). This approach emphasizes the rheological inhomogeneity of YSF and incorporates local mobility as a function of shear stress and shear rate (Benzi et al. 2021). Another possibility is that a YSF flows after $h_{P}$ as a ‘gliding solid block’. The experiments can be explained by the presence of a shallow lubrication layer near the wall, where high shear stress is concentrated. The slip of YSFs is common, regardless of their yielding mechanism, chemical composition, geometry or wall surface roughness (Piau 2007; Damianou et al. 2014; Zhang et al. 2017; Jalaal et al. 2018).

In order to verify the aforementioned hypotheses, the flow velocity is visualized with the particle image velocimetry (PIV) method. We expect a plug flow in the case of an entire solid, otherwise a curved velocity profile for a mixed solid–fluid state. An example optical image is shown in figure 2(a). Experimental details for PIV measurement are given in Appendix D. The velocity profile was analysed during horizontal imbibition due to its convenience of imaging. This method remains valid as a YSF solidifies even in the absence of gravity, and the penetration curve exhibits a plateau similar to that of the vertical rise (figure 2b). The velocity profile depicted in figure 2(c) is indeed a plug, suggesting that the fluid undergoes complete solidification instead of existing in a state of mixture between solid and liquid.

Figure 2. (a–c) Slip demonstration of horizontal imbibition using 0.12 wt. $\%$ Carbopol. (a) Microscope image of 0.12 wt. $\%$ Carbopol inside a glass capillary. The solution is seeded with 10 $\mathrm {\mu }$m hollow glass spheres. The scale bar is 100 $\mathrm {\mu }$m. The scheme on the right describes the depletion layer at the YSF–glass interface. (b) Horizontal capillary imbibition curve. (c) The PIV analysis result. The solid line is the plug flow velocity profile from the bulk imbibition speed. (d,e) Stress-velocity plots during the vertical capillary imbibition for (d) 0.12 wt. $\%$ and (e) various YSFs. The solid line is the Newtonian slip model (2.4). The shading in (d) shows the error bar. Insets in (c,d) show the vertical capillary rise curve $h(t)$, and coloured markers denote corresponding regimes. Velocity in (e) is normalized with $\mu$ and $l_{s}$. The concentrations in the legend in (e) are given in wt. $\%$.

To confirm the slip effect, the wall shear stresses are examined from all the points after the plateau. Because the wall shear stress can be computed simply by the rise ($h$) using (1.1), it is independent of the constitutive relationship to the fluid viscosity. In figure 2(d), when the height of the vertical capillary rise is converted into the wall shear stress, it shows the wall shear stress increases linearly with the rise velocity. Even if the experimental stress at the regime with high velocity shows a higher value than the theoretical prediction, it lies in the theoretical model when the lower bound of the error bar is accounted for. Moreover, the same relation has been found even more clearly in other vertical capillary rises of xanthan solutions (grey markers in figure 2e), as well as other YSFs shown in Zhang et al. (2017). In other words, even if the apparent viscosity changes, YSFs share similar slip behaviour at boundaries.

To implement slip, a Newtonian slip model (Zhang et al. 2017) is added to (2.1). The slip layer has thickness $l_{s}$, in which polymers are depleted. When its viscosity is assumed as the solvent viscosity ($\mu$), the slip velocity can be written with wall shear stress proportionally as

(2.3)\begin{equation} \dot{h}_{slip}(t)=\frac{l_{s}\sigma_{Y}}{\mu\,\epsilon(t)}, \end{equation}

which is added to $\dot {h}_{HB}(t)$ as

(2.4)\begin{equation} \dot{h}(t)=\dot{h}_{slip}(t)+\dot{h}_{HB}(t). \end{equation}

Even though other slip models have been suggested in Damianou et al. (2014), it has been demonstrated that these various models converge to a linear model when the edge evaporation effect is removed (Zhang et al. 2017). Therefore, in our evaporation-controlled study, we consider the linear model, which still shows great accordance with experimental data within the error range.

Equation (2.4) is connected with (1.1), and the value of $l_{s}$ is determined by fitting the slope in the stress–velocity plot, the linear plot of figure 2(d). The slip gradually develops even before solidification when $h< h_{p}$, but the slip's effect on the rising speed is negligible as $\sigma _{Y} \ll \sigma _{w} (t)$. After the solidification, the slip itself takes on the role of the apparent speed as the YSF is not flowing, considering that the contribution from (2.1) is zero with $\sigma _{w} (t)<\sigma _{Y}$. Due to the discontinuity of (2.1) between regimes before and after $h_{p}$, the numerical calculation explicitly connects the two regimes. The final solution of the early regime, which includes $h$, $\dot {h}$ and $\sigma _{w}$, serves as the initial condition after $h_{p}$.

In figures 3(a,b), the solid lines represent the capillary rise prediction of (1.1) and (2.4) for different YSFs. The experimental data are plotted using colour markers to indicate the yield stress. The rheological properties, surface tension and $l_{s}$ are given in tables 1 and 2. Surprisingly, model curves are capable of effectively capturing the diverse dynamics of capillary rise. The transition from solidification to gliding flow is accurately described smoothly, as supported by experimental results.

Figure 3. (a,b) Proposed model prediction (solid lines) against experiments. From top to bottom, concentrations are 0.08 and 0.12 wt. $\%$ in (a), and 0.3, 0.7, 1.0 and 1.3 wt. $\%$ in (b). Marker colours match the colour map for $\sigma _{Y}$. (c,d) Normalized numerical solutions of (1.1) and (2.4) varying yield stress and slip conditions. From the bottom, (c) $\sigma _{Y} = 40$, 10 and 2 Pa, while (d) $(l_{s},\phi ) = (0.01\ {\rm nm},0.006)$, (0.2 nm, 0.1) and (16 nm, 9) nm. Dashed lines in (c,d) are shear thinning curves without slip. (e) Contour plot of $\varDelta$. The solid line shows $l_{s} = 10$ nm, and the dashed line is $l_{s}=0.1$ nm. The lower inset is the graphical definition of $\varDelta$, while the upper inset shows the same contour over $l_{s}$ versus $\sigma _{Y}$. For (c–e), rheological properties other than $\sigma _{Y}$, surface tension and DCA parameters follow the Carbopol 0.12 wt. $\%$ values shown in table 1 and table 2.

To better reflect the initial trajectory of 1.3 and 1.0 wt. $\%$ xanthan solutions, the yield stress was fitted empirically from the rise curve. This is due to the ambiguity in the yield stress determination for xanthan (Barnes & Walters 1985; Tang et al. 2018). The residual discrepancy can be explained with shear banding, as shown in Appendix E. By reducing the shear banding effect with salt (Tang et al. 2018), experimental data follow the proposed model very well. To be within the scope of this work, accurate prediction of shear banding (de Souza Mendes 2011) is not included here. Although the Carbopol 0.08 wt. $\%$ sample's yield stress is also fitted, a plateau curve appears and could be interpreted with our proposed model, when the surface is altered into a less slippery condition through surface silanization (see Appendix E and figure 6b). The lowest concentration for which Carbopol maintains yield stress characteristics is 0.08 wt. $\%$, as highlighted by Jaworski et al. (2022), suggesting potential shear thinning behaviour under rapid flow conditions. Under conditions of a slippery surface, initial penetration speed is notably quicker compared to scenarios involving a silanized surface, where the pace significantly diminishes. The yield stress values $\sigma _{Y}$, deduced from the $h(t)$ profiles, are used for the theoretical prediction in figures 3(a) and 3(b), with values listed in table 1. Concurrently, the table also presents $\sigma _{HB}$ values, derived from rheological analyses detailed in Appendix B.

The plateau shows diverse intensities and shapes in both our experiments and the theoretical model. It can be expected that a larger yield stress results in a stronger plateau. In figure 3(b), the plateau at 1.3 wt. $\%$ xanthan solution in the darkest colour marker, however, is less pronounced and nearly unnoticeable as compared to the lower concentrations. It can be inferred that the strength of a plateau might be determined not solely by the yield stress, but also by other variables such as slip and consistency parameter ($K$).

Based on the proposed model, the various intensities and shapes of the plateau are studied theoretically. Capillary rise curves were obtained in figures 3(c) and 3(d) by changing only the rheology or slip, respectively, in order to distinguish the effects from each factor. All axes are normalized, where $\hat {h}$ is the rise height divided by the final height ($h_{E}$), and $\hat {t}$ is the time divided by the characteristic time scale of capillary rise ($t_{v}$, which is defined without yield stress and given in table 2) (Kim et al. 2020).

Table 2. DCA model parameters, slip dominance and plateau intensity.

In figure 3(c), an increase in the yield stress slows down the rise speed, and the plateau appears earlier. Figure 3(d) depicts the variation of slip strength while assuming a constant yield stress. In order to assess the contribution of slip at the plateau, the velocity ratio ($\phi$) is obtained by dividing the slip contribution ($\dot {h}_{slip}$) by the non-slip contribution ($\dot {h}_{HB}=\dot {h}-\dot {h}_{slip}$). At the point where the slip and non-slip velocities experience the same stress around the plateau, ($\sigma _{w}=2\sigma _{Y}$), we have

(2.5)\begin{equation} \phi={\frac{\dot{h}_{{slip},\epsilon=0.5}}{\dot{h}_{{HB},\epsilon=0.5}}}= \frac{8\,g(n)\,l_{s}}{R\mu } {\left(\frac{K}{{\sigma_{Y}}^{1-n}} \right)}^{1/n}, \end{equation}

where $g(n)={(n+1)(2n+1)(3n+1)}/{n(7n^2+8n+2)}$. The two graphs (figures 3c,d) are completely different, suggesting that each has a unique impact on the capillary rise dynamics. The position of the plateau and the transitional rate of rise are determined mainly by the yield stress, whereas the persistence of solidification and the rate of later rise are affected by slip.

The time scale for the duration of the plateau ($t_{p}$) is derived via dividing plateau height by slip velocity, as

(2.6)\begin{equation} t_{p}=\frac{h_{p}}{\dot{h}_{{slip},\epsilon=1}} =\frac{h_{E} \mu}{l_{s}\sigma_{Y}\left(1+\dfrac{2\sigma_{Y}}{\rho gR}\right)}. \end{equation}

We define $\hat {t}_{p}$ as $t_{p}$ divided by $t_{v}$, shown by dots on the solid lines in figures 3(c) and 3(d). Solidification finishes faster for higher yield stress and more slippery cases. It is interesting to note that when $l_{s}$ is reduced below a single atomic size, $t_{p}$ extends beyond $10^3$ s. According to Zhang et al. (2017), this duration is sufficient for the evaporation-induced residual stress to build up at the edge. The evaporation at the edge slows down the capillary imbibition speed because it affects primarily the wetting friction of the meniscus (Zhao et al. 2019). In Young's spreading of YSFs on rough surfaces, a similar cessation occurs when the plateau phase persists long enough for the meniscus edge to fully dry out (Géraud et al. 2014). Our observation of additional flow and its interpretation by incorporating the slip model distinguish our findings from the previous models and experiments (Balmforth 2022), where lubrication theory is applied identically in Hele-Shaw cell geometry. On the other hand, the previous no-slip asymptotic rise curve can be derived from the present model by assuming an extremely short $l_{s}$ and minimizing the DCA effect.

For a general understanding of the plateau intensity, we evaluated the flexion of plateau ($\varDelta$) as

(2.7)\begin{equation} \varDelta = \int_{0}^{\infty}({\tilde{h}_0(\tilde{t})- \tilde{h}(\tilde{t})})\,{\rm d} \tilde{t}, \end{equation}

where $\tilde {t}$ and $\tilde {h}$ are logarithmic $\hat {t}$ and $\hat {h}$, and $\tilde {h}_0$ is the normalized logarithmic rise height with $\sigma _{Y}=0$. As shown in the bottom inset $h(t)$ curve of figure 3(e), $\varDelta$ is the integration of the difference in height between curves $\tilde {h}_0$ and $\tilde {h}$. The intensity is calculated by varying $\sigma _{Y}$ and $l_{s}$ as shown on the plane of the normalized plateau time and height. According to the presented phase diagram, it can be observed that the plateau strength is positively correlated with the magnitude of the yield stress, while exhibiting an inverse relationship with the degree of slip. Even if the solidification occurs early (low $\hat {h}_{p}$), the flow can closely resemble the dynamics of shear thinning fluids when they experience a large slip effect (low $\hat {t}_{p}$), as well as in high $\hat {h}_{p}$ cases. The enhancement of slip can be obtained in cases of the increasing thickness of the slip layer (e.g. the enhancing repulsive forces between the YSF element and the smooth wall).

3. Conclusions

This study revisits the capillary imbibition problem in yield stress fluids (YSFs), demonstrating that YSFs inevitably solidify during wetting or spreading flows as the velocity gradually diminishes. The liquid–solid transition appears as a plateau at the rise curve, at a glance, but its re-rise occurs by a slip between the YSF and the inner wall surface. By assuming a depleted slip layer, the entire capillary rise curve could be predicted with a simple rheology model. The MATLAB code utilized for numerical calculations is readily accessible. Even if our demonstration is performed under ideal conditions without evaporation, the suggested model including DCA theory can provide crucial information about how YSFs flow inside channels at lower stresses, which is difficult to measure accurately, typically requiring a complicated force sensor or a micro-PIV device (Graziano et al. 2021). Our capillary imbibition model with experiments can be helpful to obtain parameters such as surface tension, slip length and dynamic wetting coefficients, required for wetting, spreading, simple pipe flow and any other slippery flow. Since predicting penetration rate and fluid uptake in porous media shares the same phenomena with capillary rise, our system can also provide clues for modelling YSF flow in porous media.