2.1. Experimental materials, apparatus and procedure

The experiments involve the downward injection of a heavy viscoplastic fluid into a light in situ Newtonian one (see figure 1 and table 2 for the flow parameters). Motivated by P&A applications, in general, an eccentric configuration is considered by installing an eccentric inner pipe inside a vertical closed-end outer pipe (although a limited number of concentric experiments are also performed for comparison purposes). The inner pipe is connected to a gear pump (Ismatech 405A) through which the viscoplastic fluid is injected at a precise flow rate. For the Newtonian fluid filling the outer pipe, water–glycerol mixtures are used. For the injected viscoplastic fluid, water–glycerol mixtures with the addition of Carbopol (Carbomer 940, Making Cosmetics Co.) at various concentrations are employed. The preparation of our viscoplastic fluid follows the established methods in the field (Eslami & Taghavi Reference Eslami and Taghavi2017). For visualization purposes, the viscoplastic fluid is dyed with 800 (mg l$^{-1}$) of ink (Fountain Pen India black ink). The densities are measured using a high-accuracy density meter (Anton Paar DMA 35).

Figure 1. Schematic of the experimental apparatus, wherein a viscoplastic fluid is injected into a vertical closed-end pipe filled with a Newtonian fluid.

We perform the rheological measurements of our Carbopol solutions using a rheometer (DHR-3, TA Instruments) with a parallel-plate geometry, with a diameter of 40 (mm) and mean gap of 1 (mm). Fine sandpapers (as rough surfaces) are attached to the rheometer plates to eliminate any possible wall-slip effects in the measurements at low shear rates (Habibi et al. Reference Habibi, Dinkgreve, Paredes, Denn and Bonn2016; Roberts & Barnes Reference Roberts and Barnes2001). We obtain the steady-state flow curves of our Carbopol solutions in experiments with controlled shear rate ($\hat {\dot {\gamma }}$), where the samples are sheared in an upward ramp of $10^{-2}<\hat {\dot {\gamma }}< 10^{3}$ (s$^{-1}$), as shown in figure 2(a). The figure shows the rheometry data represented by the symbols. The flow curves (figure 2a) display the characteristics of the materials above the critical stress. At low shear rates, the shear stress plateaus to a finite value corresponding to the yield stress ($\hat \tau _y$). In addition to the yield stress, the samples feature shear-thinning behaviours. As shown, the data points can be well fitted with the lines corresponding to the HB model parameters (Jaworski et al. Reference Jaworski, Spychaj, Story and Story2021). The model, fitted by the experimental data, has three parameters, i.e. the fluid consistency $\hat {\kappa }$, the yield stress $\hat {\tau }_y$ and the power-law index $n$, defined as

(2.1)\begin{equation} \left. \begin{array}{ll@{}} \hat{\dot\gamma}=0, & |{\hat \tau}| \le {\hat \tau _y},\\ \hat \tau = {\hat \tau _y} + \hat{\kappa} {\hat{\dot \gamma} ^{n}}, & |{\hat \tau}| > {\hat \tau _y}. \end{array} \right\}\end{equation}

Figure 2. (a) Flow curves of shear stress $\hat \tau$ versus shear rate $\hat {\dot {\gamma }}$ from the shear-rate-controlled tests. The symbols mark Sample I (${\blacksquare }$), Sample II ($\blacktriangleright$), Sample III ($\bullet$), Sample IV ($\blacktriangle$), Sample V ($\bigstar$) and Sample VI ($\blacktriangleleft$). The green dashed lines correspond to the HB fits, (2.1), with the rheological parameters that are reported for each fluid in table 1. (b) The results of the oscillation amplitude sweep tests for Samples II, III, V and VI. The symbols and dotted lines correspond to the storage ($\hat G^{\prime }$) and loss ($\hat G^{\prime \prime }$) moduli, respectively. The vertical lines show the intersection of $\hat G^{\prime }$ and $\hat G^{\prime \prime }$ as an indication of $\hat \tau _y$. (c) Shear stress as a function of shear strain at constant shear rates for Sample II at $\hat {\dot {\gamma }}= 0.0001$, 0.001, 0.01 and 0.1 s$^{-1}$ indicated by square, diamond, triangle and circle symbols, respectively. The scattered data at very low shear strains ($\gamma <10^{-2}$) are due to instrument artifacts. (d) Creep tests showing shear rate as a function of time for Sample II at $\hat \tau = 1$, 7, 13 and 19 Pa, marked by square, diamond, triangle and circle symbols, respectively. The shear rate variations at early times are due to instrument artifacts. Note that, to make the presented data visually distinguishable, the results for some samples are depicted, while the other sample data (not shown) exhibit similar trends. The size and the colour intensity of the symbols represent the following: in (a,b) the Carbopol concentration (i.e. larger and darker symbols represent higher Carbopol concentrations), in (c) the power of the applied shear rate and in (d) the intensity of the applied stress.

The rheological parameters of our Carbopol solutions, obtained using the HB model, are given in table 1 for the fluids depicted in figure 2(a). As expected, at smaller (larger) Carbopol concentrations, the yield stress values are lower (higher), i.e. the solid-like behaviour of the material is less (more) pronounced, and the required stress to result in yielding is lower (higher).

Table 1. The rheological parameters of Carbopol solutions determined from fitting the shear stress–shear rate curves to the HB model: the yield stress ($\hat \tau _y$), the consistency index ($\hat \kappa$) and the power-law index ($n$). The Carbopol concentration decreases from Sample I to Sample VI.

Table 2. The dimensional parameters and their ranges in our experiments.

Although our rheological model choice is the classical HB model, there are also alternative procedures to fit viscoplastic flow curves, e.g. by taking the HB model as a base. An example is the three-component (TC) model (Caggioni, Trappe & Spicer Reference Caggioni, Trappe and Spicer2020), defined by $\hat \tau = {\hat \tau _y} + {\hat \tau _y}{( {\hat {\dot \gamma } /{{\hat {\dot \gamma }}_c}} )^{1/2}} + {\hat \eta _{bg}}\hat {\dot \gamma }$, where $\hat {\dot \gamma }_{c}$ and $\hat {\eta }_{bg}$ are the critical shear rate and the viscosity of the continuous phase, respectively. The TC model represents three regimes of deformation, i.e. the elastic, plastic and viscous regimes. At low shear rates in this model, elastic straining determines the stress and, when the shear rate exceeds a critical value ($\hat {\dot \gamma }_{c}$), a transition takes place from elastic to plastic. Since the yield stress is an utmost important parameter in our study, we have also used the TC model to extract its value from the rheological data (results omitted for brevity), finding that the yield stresses obtained by the TC model are on average less than 10 % lower than those fitted by the HB model. Such closeness of yield stress values using these two models has been also reported by Caggioni et al. (Reference Caggioni, Trappe and Spicer2020). That said, in this study, we still use the classical HB model, as a simple, common, reliable and accurate model (Frigaard Reference Frigaard2019) for describing the rheological behaviours of our Carbopol solutions, in particular the yield stress values.

For the sake of completeness of the rheological characterization, we also perform oscillatory rheometry tests, to characterize elastic and viscous behaviours of our samples at a fixed frequency of 1 (Hz), with applied shear stresses ($\hat {\tau }$) ranging from 0.01 to 200 (Pa). This allows us to evaluate the storage modulus ($\hat G^{\prime }$) and loss modulus ($\hat G^{\prime \prime }$), which represent the elastic and viscous characteristics, respectively. As shown in figure 2(b), initially (at small $\hat {\tau }$) $\hat G^{\prime }$ is found to be higher than $\hat G^{\prime \prime }$ for all the sample solutions. This indicates a solid-like behaviour before yielding. For each sample, when the stress exceeds a critical value, an increase in $\hat G^{\prime \prime }$ accompanied by a decrease in $\hat G^{\prime }$ is observed. Such behaviour has been reported to occur in soft glassy materials (Hyun et al. Reference Hyun, Wilhelm, Klein, Cho, Nam, Ahn, Lee, Ewoldt and McKinley2011), such as Carbopol gels (Ovarlez, Barral & Coussot Reference Ovarlez, Barral and Coussot2010). Beyond this critical stress, the viscous response of the material overcomes the elastic one, and both $\hat G^{\prime }$ and $\hat G^{\prime \prime }$ vary with $\hat {\tau }$. Here, $\hat G^{\prime }$ significantly and $\hat G^{\prime \prime }$ slightly decrease versus $\hat {\tau }$, in such a way that $\hat G^{\prime \prime }$ becomes greater than $\hat G^{\prime }$, indicating that the material is indeed flowing and implying negligible elastic effects. Although the elastic effects of Carbopol solution may be important in certain situations, in this study, we neglect these effects on the flow patterns and dynamics.

The vertical lines in figure 2(b) represent the stress values at which the $\hat G^{\prime }$ and $\hat G^{\prime \prime }$ curves cross over, implying a rough estimation to the yield stress, as an alternative method (however, note that, as argued by Fernandes et al. (Reference Fernandes, Andrade, Franco and Negrão2017), estimating the yield point using the crossover must be made with caution). One can also compare these yield stress values with those reported in table 1. The values obtained from the HB model fitting in figure 2(a) are slightly different from the ones computed at the crossover of the $\hat G^{\prime }$ and $\hat G^{\prime \prime }$ curves. This is in line with recent measurements of similar types of viscoplastic fluids (Dinkgreve et al. Reference Dinkgreve, Paredes, Denn and Bonn2016; Fernandes et al. Reference Fernandes, Andrade, Franco and Negrão2017; Jalaal, Kemper & Lohse Reference Jalaal, Kemper and Lohse2019).

To further complete our rheological characterization, figure 2(c) presents an example of our results in terms of the shear stress as a function of the shear strain for different applied shear rates. Each curve presents the average of three rheological tests. At low shear strains, the stress curves are close to one another, implying nearly constant properties, regardless of the applied shear rate, hinting at a linear viscoelastic response (Fernandes et al. Reference Fernandes, Andrade, Franco and Negrão2017). At higher shear strains, however, the shear stress in each curve rises up as a function of the shear strain, and then deviates from the nearly linear viscoelastic region, to reach a maximum value associated with yielding and rupturing (Bonn et al. Reference Bonn, Denn, Berthier, Divoux and Manneville2017; Fernandes et al. Reference Fernandes, Andrade, Franco and Negrão2017). The maximum shear stress, however, increases with an increase in the applied shear rate value. On the other hand, the yield stress obtained from the shear stress–shear rate flow curve (figure 2a, table 1) is slightly lower than the maximum stress at the lowest shear rate in figure 2(c) ($\hat {\dot {\gamma }}=10^{-4}$ ($\textrm {s}^{-1}$)).

As a final rheological characterization method, figure 2(d) presents the results of the creep tests for a given sample, by imposing a constant shear stress and recording the shear rate as a function of time. As seen, for $\hat \tau = 1$ (Pa) and $\hat \tau = 7$ (Pa), the shear rate approaches zero at long time, implying that the applied stress is lower than the yield stress, i.e. $\hat \tau _y = 10.25$ (Pa) evaluated from the shear stress–shear rate flow curve (figure 2a, table 1). On the other hand, for $\hat \tau = 13$ (Pa) and $\hat \tau = 19$ (Pa), the shear rate tends to a constant non-zero value (Espinoza et al. Reference Espinoza, Varges, Rodrigues, Naccache and de Souza Mendes2022), implying that the applied stress is above the yield stress, i.e. an observation generally consistent with the previous measurements.

Note that we consider our experimental pair of fluids (i.e. Carbopol solution and water) to be fully miscible, which means that there are no surface tension effects in our work. This consideration is in line with recent studies demonstrating that the difference between the surface tension of Carbopol solutions and pure water is very small (less than 10 %), and this difference does not vary significantly with the Carbopol concentration (Boujlel & Coussot Reference Boujlel and Coussot2013; Jørgensen et al. Reference Jørgensen, Le Merrer, Delanoë-Ayari and Barentin2015; Jalaal et al. Reference Jalaal, Kemper and Lohse2019).

Before running each experiment, a movable piston is used to adjust the falling height ($\hat H$). Then, the inner and outer pipes are filled with the viscoplastic and Newtonian fluids, respectively, while a gate valve at the inlet initially separates the two fluids. Each experiment begins with opening the gate valve and injecting/extruding the viscoplastic fluid into the Newtonian fluid at a fixed flow rate; simultaneously, the flow images are obtained using a camera (Basler acA2040) at 20 frames per second and transmitted to a computer for subsequent post-processing (via MATLAB). As the injection flow behaviour is expected to be observed below the inlet area, our camera's field of view is focused on this region. A large number of experiments (${>}300$) are performed for different flow rates, falling heights, rheological parameters and density and viscosity ratios.

2.2. Experimental observations of three distinct flow regimes

In this section, we describe the general observations made in our experiments. As mentioned before, our experiments are focused on analysing the viscoplastic fluid flow below the inlet area. For each experiment, the flow rate is chosen and the falling height is fixed. Upon opening the gate valve, which initially separates the fluids, the heavy fluid enters the flow domain below the gate valve and penetrates into the in situ light fluid.

Via our experiments, we have identified the following three distinct flow regimes during the injection of the heavy viscoplastic fluid into the light Newtonian fluid (as demonstrated in figure 3):

1. (i) Breakup regime. A flow example in this regime is illustrated in figure 3(a $_1$). The breakup regime is typically observed at small injection rates (small $\hat {V}_0$), large buoyancy forces, small to moderate yield stresses of the viscoplastic fluid, small viscosities of the in situ fluid and large falling heights. In this regime, the viscoplastic core fluid is initially extruded as a stable cylinder from the inlet and it advances towards the pipe end. At short times and axial lengths of the core, the stress in the viscoplastic core fluid is below the yield stress, so that a long stable filament can advance downward, holding up against the gravitational forces that attempt to extend it more rapidly. While propagating downward, the core fluid axial length naturally increases and, at a certain point, the core fluid can no longer resist the progressively increasingly buoyancy, which overcomes the yield stress; the core fluid starts to yield and thin near the inlet and it is eventually broken/ruptured. Increasing the buoyancy force (e.g. by increasing the density difference) promotes this flow regime and eases the breakup of the core fluid. After the break of the filament, its upper part which is still connected to the injection inlet is restored fairly rapidly to form a new filament. If the injection continues for a long time in this regime, several filaments can be formed and subsequently break up; the broken parts fall down evenly and make a structure of broken filaments at the bottom end of the pipe.

2. (ii) Coiling regime. A flow example in this regime is illustrated in figure 3(a $_2$). The coiling regime is typically observed at moderate injection rates (moderate $\hat {V}_0$), moderate to large yield stresses of the viscoplastic fluid and relatively short falling heights. In typical situations in this regime, the viscoplastic core fluid is initially extruded, as a stable cylinder, with a diameter that remains typically unchanged as the filament front moves downward. In most cases, the core fluid front can reach the pipe end before the core fluid is broken; then, its front sits on the pipe end and starts coiling evenly, with a fixed centre and a constant coiling frequency (termed regular coiling). As the heavy viscoplastic fluid is injected continuously, the viscoplastic fluid layers are positioned on top of one another, with a thickness initially equal to the nozzle diameter and the coiling of the core fluid gradually continues up to the inlet area. At large velocities, secondary or irregular behaviour (termed irregular coiling) may be observed; at small buoyancy, in some cases the coiling may start to occur before the viscoplastic fluid front reaches the pipe end (termed free coiling). For simplicity, we collectively classify all these as the coiling regime. In general, the coiling behaviour observed in our experiments is reminiscent of the coiling of a viscous fluid filament or an elastic rope falling onto a solid surface, as studied by previous works (Habibi et al. Reference Habibi, Maleki, Golestanian, Ribe and Bonn2006; Ribe et al. Reference Ribe, Habibi and Bonn2012; Ribe Reference Ribe2017), explaining the coiling mechanism via a balance among viscous, gravitational, inertial and elastic stresses.

3. (iii) Buckling regime. A flow example in this regime is illustrated in figure 3(a $_3$). The buckling (bulging) regime is typically observed at large injection rates (large $\hat {V}_0$), small buoyancy forces, low yield stresses of the viscoplastic fluid and high viscosities of the in situ fluid. In this regime, as the heavy core fluid penetrates into the in situ light fluid, it undergoes axial compressions, due to the pressure gradient and shear stresses at the interface, applied by the resistive upward flow of the annular fluid, resulting in the yielding of the core fluid. As the yield stress is overcome, the core fluid buckles and it quickly expands radially. With time, the radial expansion continues until the injected fluid completely blocks the region near the inlet boundary, and impedes the removal of the light in situ fluid from the bottom of the pipe. Increasing the injection velocity and viscosifying the in situ fluid promote the viscous compressional stresses towards yielding and buckling of the core fluid. Note that, in the context of our results throughout the text, the buckling regime concerns a radial expansion and buckling of the viscoplastic core fluid and, therefore, can be also called the bulging regime.

Figure 3. (a) Injection of viscoplastic fluid into closed-end pipe filled with Newtonian fluid through an eccentric inner pipe ($\hat {E}=6.35$ mm): ($a_1$) breakup regime for $\hat V_0=5.3$ mm s$^{-1}$, $\hat \rho _a= 998$ kg m$^{-3}$ and $\hat {\mu }_a= 0.001$ $\textrm {Pa}\, \textrm {s}$; ($a_2$) coiling regime for $\hat V_0=25.5$ mm s$^{-1}$, $\hat \rho _a= 998$ kg m$^{-3}$ and $\hat {\mu }_a= 0.001$ $\textrm {Pa}\,\textrm {s}$; ($a_3$) buckling (bulging) regime for $\hat V_0=87.5$ mm s$^{-1}$, $\hat \rho _a=1007.5$ kg m$^{-3}$ and $\hat \mu _a=0.025$ $\textrm {Pa}\,\textrm {s}$. For all the snapshots, $\hat H=200$ mm and the field of view is $20\times 3.81$ cm$^{2}$. (b) Injection through a concentric inner pipe ($\hat {E}=0$ mm): ($b_1$) breakup regime for $\hat V_0=7.3$ mm s$^{-1}$; ($b_2$) coiling regime for $\hat V_0=40.4$ mm s$^{-1}$; ($b_3$) buckling (bulging) regime for $\hat V_0=117.5$ mm s$^{-1}$. For all the snapshots, $\hat H=100$ mm, $\hat \rho _a= 997.5$ kg m$^{-3}$ and $\hat {\mu }_a= 0.001$ $\textrm {Pa}\,\textrm {s}$ and the field of view is $10\times 3.81$ cm$^{2}$. In all panels, the injection fluid is Sample IV (see table 1), ${\rm \Delta} \hat t$ is the time span between the snapshots, the injection inlet is indicated by an arrow and the outlet area of the flow domain is marked by a red surface.

In general, the flow regimes observed and explained above can be also visualized in a concentric flow configuration. To illustrate this, we have conducted a number of experiments in a second apparatus, made of two concentric pipes, for which the results are shown in figure 3(b). As seen, the observed flow regimes in the injection of a viscoplastic fluid into the in situ fluid using the concentric inner pipe are similar to the regimes that occur in the eccentric case (compare top and bottom panel groups in figure 3).

Let us clarify the criterion for classifying our experimental observations into the three defined regimes. To be systematic, this is based on analysing the variation of the viscoplastic filament radius in the upper half of the domain, which renders consistent measurements. To do so, the variation of the ratio between the minimum and maximum radii of the filament, $\hat R_{min}/\hat R_{max}$, is quantified versus time, as exemplified in figure 4(a). If this ratio ($\hat R_{min}/\hat R_{max}$) tends to zero at long time, implying that the core fluid is yielded towards breakup, we quantify the flow as the breakup regime. When the core buckles (bulges), $\hat R_{min}/\hat R_{max}$ tends to a constant value, i.e. $\hat R_{min}/\hat R_{c}\approx 1/3$ at long times (marked by the dashed line in figure 4a). For the experiments in which $\hat R_{min}/\hat R_{max}$ does not reach zero or the constant value of $\sim 1/3$ at long times (and in fact $\hat R_{min}/\hat R_{max} \approx 1$), the type of flow regime is classified as the coiling regime. Our classification criterion/approach allows us to circumvent various flow complexities and secondary flow behaviours, e.g. in the coiling regime (which may exhibit regular, free and irregular coiling), providing a simple way to classify our complex flow in terms of the main behaviours.

Figure 4. Variation of the core fluid minimum radius identified by $\hat R_{min}/\hat R_{max}$ versus the dimensionless flow time ($t={\hat t \hat V_0}/{\hat R}$) for Sample V and $\hat {E}=6.35$ mm. (a) A typical breakup regime ($\square$) for $\hat V_0= 8.2$ mm s$^{-1}$, $\hat \rho _a= 998$ kg m$^{-3}$, $\hat {\mu }_a= 0.001$ $\textrm {Pa}\,\textrm {s}$ and $\hat H=400$ mm; a typical coiling regime ($\circ$) for $\hat V_0= 34.0$ mm s$^{-1}$, $\hat \rho _a= 998$ kg m$^{-3}$, $\hat {\mu }_a= 0.001$ $\textrm {Pa}\,\textrm {s}$ and $\hat H=400$ mm; a typical buckling regime ($\vartriangle$) for $\hat V_0=87.50$ mm s$^{-1}$, $\hat \rho _a= 1007.5$ kg m$^{-3}$, $\hat {\mu }_a= 0.0252$ $\textrm {Pa}\,\textrm {s}$ and $\hat H=260$ mm. (b) An example of coiling case where the viscoplastic core fluid starts to yield but, due to the finite-sized pipe length, its front reaches the pipe end and it is finally compressed; thus, $\hat R_{min}/\hat R_{max}$ initially decreases but recovers. The results correspond to $\hat V_0=34.0$ mm s$^{-1}$, $\hat \rho _a= 998$ kg m$^{-3}$, $\hat {\mu }_a= 0.001$ $\textrm {Pa}\,\textrm {s}$, $\hat H=100$ mm and $\hat {E}=6.35$ mm. The inset shows the initial yielding, the recovering and the eventual coiling of the core fluid, for $t=[3.68~4.91~5.35~5.80~7.84]$ from left to right.

We have experimentally observed that, in some cases categorized within the coiling regime, as the flow develops the core fluid is actually yielded and starts thinning, while its radius decreases; however, the core fluid front eventually reaches the pipe end before any breakup can occur and starts coiling. A typical example of such a case is presented in figure 4(b). As can be seen, the value of the core radius recovers after the initial decrease. Note that increasing the falling height ($H$) in such cases would lead to the breakup of the core. However, since the breakup does not actually occur for the used experimental falling height, we still classify these flows as the coiling regime (while highlighting the corresponding results in the results section).

It may be worth mentioning that, in terms of the breakup regime, previous relevant studies (e.g. those considering the injection of a viscoplastic fluid into air) have reported a strong conical shape of the core at the breakup point, which also highly depends on the capillary effects and the core density (Balmforth et al. Reference Balmforth, Dubash and Slim2010b; Moschopoulos et al. Reference Moschopoulos, Syrakos, Dimakopoulos and Tsamopoulos2020; Geffrault et al. Reference Geffrault, Bessaies-Bey, Roussel and Coussot2021). However, in our miscible buoyant flow, the surface tension effects are neglected, the density differences are small and the core plasticity is more pronounced, resulting in a relatively more flat shape at the breakup point (e.g. figure 3a). This may be conceptually in line with the finding of Balmforth et al. (Reference Balmforth, Dubash and Slim2010b) who have shown that, when the yield stress effects dominate the surface tension ones, the breakup point shape changes to a flatter end cone.

Before we proceed, it may be worth mentioning that viscoplastic fluids are known to slip on smooth solid surfaces, for example glass plates (Jalaal, Balmforth & Stoeber Reference Jalaal, Balmforth and Stoeber2015) and hydrophobic surfaces (Rahmani & Taghavi Reference Rahmani and Taghavi2020). However, in this study, we do not expect the viscoplastic Carbopol flow in the pipe to be influenced by slippage on the surface of the outer pipe, since there is no direct contact between the injected viscoplastic fluid and the outer pipe wall. In fact, in the time scale of our interest, the in situ (Newtonian) fluid always exists and separates the injected viscoplastic fluid and the pipe wall. Considering this argument, we do not expect the pipe wall properties to significantly change the experimental observations, e.g. by causing slippage effects.

Based on the experimental observations, it is clear that the development of a multilayer viscoplastic model is appropriate to gain a deeper understanding of our results and provide predictions to the transition between the aforementioned flow regimes; this also allows an analysis of the effects of certain flow parameters (e.g. the eccentricity or the inner pipe radius) that cannot be easily changed in our experiments. Thus, in what follows, we develop such a model, which will include the key governing dimensionless parameters of table 3, to generalize our results.

Table 3. The key dimensionless parameters and their ranges in this study. Note that, expect for the Reynolds number, these flow dimensionless parameters also serve as inputs to the model developed in this study.

3.1. Analysis of the viscoplastic core fluid flow

Regarding the Newtonian annular fluid in (3.18) ($k=a$), the term involving ${{\tau _{a,ZZ}}}$ is of ${O(\delta )}$ and can therefore be ignored. However, for the viscoplastic fluid ($k=c$ in (3.18)), the term involving ${{\tau _{c,ZZ}}}$ is large and cannot be ignored, as the flow is dominated by the extensional/compressional dynamics. Therefore, in the leading order, the main equation (axial momentum) for the viscoplastic fluid is

(3.24)\begin{equation} 0 ={-} \frac{{\partial {P_c}}}{{\partial Z}} + \frac{1}{{{\hbar ^{2}}}}\frac{{\partial ( {\hbar {\tau _{c,\xi Z}}} )}}{{\partial \xi }} + \frac{1}{{{\hbar ^{2}}}}\frac{{\partial ( {\hbar {\tau _{c,\phi Z}}} )}}{{\partial \phi }} + \delta \frac{{\partial {\tau _{c,ZZ}}}}{{\partial Z}} + \chi. \end{equation}

However, in general ${\tau }_{c,ZZ}$ (and accordingly the pressure) in the core phase is an unknown function. This situation results from the fact that the stress is indeterminate in the viscoplastic fluid below the yield stress and, therefore, the motion equations are not sufficient to uniquely determine all the stress components, i.e. a complexity that is frequently encountered in similar problems (Balmforth et al. Reference Balmforth, Dubash and Slim2010a). To overcome this problem, and in the interest of simplicity, we assume that ${\tau }_{c,\xi \xi } = { \tau }_{c,\phi \phi }$, which, in the spirit of the regularization of the constitutive equation (Frigaard & Nouar Reference Frigaard and Nouar2005), allows us to use the continuity equation to find

(3.25)\begin{equation} {\tau _{c,ZZ}} ={-} 2{\tau _{c,\xi\xi}} ={-} 2{\tau _{c,\phi \phi}}. \end{equation}

Considering the dominance of ${\tau _{c,ZZ}}$ also implies that

(3.26)\begin{equation} {u_{c,Z}}( {T,\xi ,\phi ,Z} ) = {u_{c,Z}}( {T,Z} ) + O( {{\delta ^{2}}} ). \end{equation}

Let us for now assume that the solution to $\tau _{c,ZZ}$ may be found. The yield condition can be obtained in the leading order as

(3.27)\begin{equation} \frac{{\sqrt 3 }}{2}| {{\tau _{c,ZZ}}} | + O( \delta ) > B. \end{equation}

The simplified yield condition criterion above is consistent with some previous theoretical approaches for viscoplastic fluid flows downwardly injected in air, such as the work of Balmforth et al. (Reference Balmforth, Dubash and Slim2010b). However, there are more recent works in the literature, such as the works of Moschopoulos et al. (Reference Moschopoulos, Syrakos, Dimakopoulos and Tsamopoulos2020) and Geffrault et al. (Reference Geffrault, Bessaies-Bey, Roussel and Coussot2021), emphasizing the importance of including shear components in the analysis, particularly because, in certain cases, the flows of our consideration can lead to yielded conical shapes even in slender filaments. Nevertheless, for the simplicity of our analysis, we assume that such shear component terms are of higher order in $\delta$ and can be ignored. To simplify the analysis, it is straightforward to apply the normal stress balance at the interface (Leal Reference Leal2007):

(3.28)\begin{equation} ( {{\sigma _a} \boldsymbol{\cdot} {\boldsymbol{n}}} )\boldsymbol{\cdot} {\boldsymbol{n}} = ( {{\sigma _c}\boldsymbol{\cdot} {\boldsymbol{n}}} )\boldsymbol{\cdot} {\boldsymbol{n}}, \end{equation}

in which $\sigma _a$ and $\sigma _c$ are the stress tensors for the annular and core phases and $\boldsymbol n$ is the unit normal vector inward to the surface:

(3.29)\begin{equation} {\boldsymbol{n}} = \frac{{\boldsymbol{\nabla} F}}{{ | {\boldsymbol{\nabla} F} |}} = \frac{{ ( {1/\hbar ,0, - \delta \partial {\xi _c}/\partial Z} )}}{{\sqrt {{{ ( {1/\hbar } )}^{2}} + {{ ( {\delta \partial {\xi _c}/\partial Z} )}^{2}}} }}, \end{equation}

and, in the leading order, we eventually find

(3.30)\begin{equation} \left. \begin{array}{c@{}} {P_c} = {P_a} - \dfrac{{\delta {\tau _{c,ZZ}}}}{2},\\ {\tau _{c,\xi Z}} = {\tau _{a,{\xi _c}Z}} + \dfrac{{3\hbar \delta {\tau _{c,ZZ}}}}{2}\dfrac{{\partial {\xi _c}}}{{\partial Z}}, \end{array} \right\} \end{equation}

where ${\tau _{a,{\xi _c}Z}} = {{{\tau _{a,\xi Z}}} |_{\xi = {\xi _c}}}$. First of all, note that since in our study both fluids are miscible (at the limit of large $Pe$), the surface tension effects are neglected in the derivations of the pressure terms above. Also, note that, again, the term involving $\delta \tau _{c,ZZ}$ is large and cannot be removed compared to the shear stress terms. We now integrate the axial momentum over the cross-section of the core fluid (${\xi _c} < \xi \le \infty ; 0 \le \phi \le 2{\rm \pi}$), while using the scale factor as well as the stress balance conditions, to find after considerable algebra

(3.31)\begin{equation} 0 ={-} R_i^{2}\frac{{\partial {P_a}}}{{\partial Z}} - \frac{1}{{\rm \pi} }\int_0^{2{\rm \pi} } {\hbar {\tau _{a,\xi_c Z}}} \, {{\rm d}} \phi + \frac{3}{2}\delta \frac{{\partial ( {{\tau _{c,ZZ}}R_i^{2}} )}}{{\partial Z}} +\chi R_i^{2}, \end{equation}

using which $\tau _{c,ZZ}$ can be obtained as

(3.32)\begin{equation} {\tau _{c,ZZ}} = \frac{2}{{3\delta R_i^{2}}}\int_Z^{L} {\left[ { - R_i^{2}\frac{{\partial {P_a}}}{{\partial Z}} - \frac{1}{{\rm \pi} }\int_0^{2{\rm \pi} } {\hbar {\tau _{a,\xi_c Z}}} \, {{\rm d}} \phi + \chi R_i^{2}} \right]\, {{\rm d}}Z}, \end{equation}

in which, for simplicity, we have applied the free-end boundary condition, i.e.

(3.33)\begin{equation} {\tau _{c,ZZ}}( L ) = 0. \end{equation}

Here, $L=L( T )$ is the core layer length, i.e. an unknown of the problem, which must be calculated for an evolving core viscoplastic layer. The first and second terms in the integral above are unknown at this point, as they represent the dynamic of the annular fluid flow; these will be calculated in the next section, but for now let us assume these terms can be obtained.

The second invariant, $\dot {\gamma }(\boldsymbol {u})$, in the leading order can be eventually found as

(3.34)\begin{equation} {\left[ {\frac{1}{2}\sum_{i,j = 1}^{3} {{{[ {{{\dot \gamma }_{ij}}} ]}^{2}}} } \right]^{1/2}} = \left| {\sqrt 3 \delta \frac{{\partial {u_{c,Z}}}}{{\partial Z}}} \right| + O( {{\delta ^{2}}} ) . \end{equation}

Therefore, the HB constitutive model can be simplified to

(3.35)\begin{equation} \left. \begin{array}{c@{}} {\sqrt 3 \left| {\dfrac{{\partial {u_{c,Z}}}}{{\partial Z}}} \right| = 0 \Leftrightarrow \dfrac{{\sqrt 3 }}{2}| {{\tau _{c,ZZ}}} | \le B},\\ {\dfrac{{\partial {u_{c,Z}}}}{{\partial Z}} = \dfrac{{{\rm{sgn}}( {{\tau _{c,ZZ}}} )}}{{\delta \sqrt 3 }}{{\left( {\dfrac{{\sqrt 3 | {{\tau _{c,ZZ}}} |}}{2} - B} \right)}^{{1}/{n}}} \Leftrightarrow \dfrac{{\sqrt 3 }}{2}| {{\tau _{c,ZZ}}} | > B}. \end{array}\right\} \end{equation}

Thus, the velocity gradient can be written as

(3.36)\begin{equation} \frac{{\partial {u_{c,Z}}}}{{\partial Z}} = \frac{{{\rm{sgn}}( {{\tau _{c,ZZ}}} )}}{{\delta \sqrt 3 }}{\left[ {\max \left( {\frac{{\sqrt 3 | {{\tau _{c,ZZ}}} |}}{2} - B,0} \right)} \right]^{1/n}}. \end{equation}

Finally, the velocity can then be obtained as

(3.37)\begin{equation} {u_{c,Z}} = \int_0^{Z} {\frac{{{\rm{sgn}}( {{\tau _{c,ZZ}}} )}}{{\delta \sqrt 3 }}{{\left[ {\max \left( {\frac{{\sqrt 3 | {{\tau _{c,ZZ}}} |}}{2} - B,0} \right)} \right]}^{1/n}}} \, {{\rm d}}Z + {u_{c,Z0}}, \end{equation}

in which $u_{c,Z0}= 1$ comes from the injection boundary condition.

In the next subsection, we analyse the Newtonian annular flow dynamics, which in return affects the viscoplastic core fluid flow development, e.g. via imposing a pressure gradient as well as a non-negligible shear stress at the interface (caused by the motion of the annular fluid). This consideration brings the multiphase nature of the flow into the analysis, highlighting that our annular fluid is dynamically active. However, such an analysis would not be relevant for dynamically inactive in-place fluids, e.g. in the simpler case of the injection of a viscoplastic fluid into air, which has been more frequently studied in the literature (Coussot & Gaulard Reference Coussot and Gaulard2005; Balmforth et al. Reference Balmforth, Dubash and Slim2010a,Reference Balmforth, Dubash and Slimb; Geffrault et al. Reference Geffrault, Bessaies-Bey, Roussel and Coussot2021).

3.2. Analysis of the Newtonian annular flow

Let us move on to analysing the annular flow, needed to calculate the stress in (3.32), starting with the leading-order axial momentum equation for the annular phase:

(3.38)\begin{equation} {0 ={-} \frac{{\partial {P_a}}}{{\partial Z}} + \frac{1}{{{\hbar ^{2}}}}\frac{{\partial ( {\hbar {\tau _{a,\xi Z}}} )}}{{\partial \xi }} + \frac{1}{{{\hbar ^{2}}}}\frac{{\partial ( {\hbar {\tau _{a, \phi Z}}} )}}{{\partial \phi }}}, \end{equation}

which using the Newtonian constitutive relation is simplified to

(3.39)\begin{equation} {f_a} = \frac{{{{( {\cosh \xi - \cos \phi } )}^{2}}}}{{{{\sin }^{2}}{\xi _w}}}\left( {\frac{{{\partial ^{2}}{u_{a,Z}}}}{{\partial {\xi ^{2}}}} + \frac{{{\partial ^{2}}{u_{a,Z}}}}{{\partial {\phi ^{2}}}}} \right), \end{equation}

where ${f_a} = ({1}/{M})({{\partial {P_a}}}/{{\partial Z}})$. Also, $u_{a,Z}$ is the leading-order solution to the velocity field in the annular phase. Note that the dependence of $u_{a,Z}( {T,\xi,\phi,Z } )$ on $( {T,Z} )$ enters only via the interface satisfying the kinematic condition, introduced further below.

Here, we extend the approach proposed by Goldstein, Ullmann & Brauner (Reference Goldstein, Ullmann and Brauner2017) to derive our analytical solution for the annular fluid flow in contact with the core viscoplastic fluid, as the inner boundary. Expectedly, the general solution of the equation is made of two solution types, i.e. the particular solution ($u_{ap,Z}$) and the homogeneous solution ($u_{ap,Z}$). The general solution becomes

(3.40)\begin{equation} {u_{a,Z}} = {u_{ah,Z}} + {u_{ap,Z}}. \end{equation}

Considering the form of the motion equations and non-homogeneous terms, following a classical style, we choose the following particular solution:

(3.41)\begin{equation} {u_{ap,Z}} = \frac{{{f_a}{{\sinh }^{2}}{\xi _w}\cosh \xi }}{{2( {\cosh \xi - \cos \phi } )}}. \end{equation}

Considering the symmetry condition with respect to the pipe centreline, $\phi = 0, 2{\rm \pi}$, i.e.

(3.42)\begin{equation} {{{\left. {\frac{{\partial {u_{ah,Z}}( {T,\xi ,\phi,Z } )}}{{\partial \phi }}} \right|}_{\phi = 0,{\rm \pi} }} = 0,} \end{equation}

it is possible to propose the following homogeneous solution in the form of the following Fourier series:

(3.43)\begin{equation} {{u_{ah,Z}} = \sum_{m = 0}^{\infty} {[ {{\alpha _{am}}\cosh ( {m\xi } ) + {\beta _{am}}\sinh ( {m\xi } )} ]} \cos m\phi}. \end{equation}

Following the procedure detailed in Appendix A, the solutions to the annular fluid flow velocity and shear stress are derived as

(3.44)\begin{gather} {u_{a,Z}} = \frac{{{f_a}{{\sinh }^{2}}{\xi _w}\cosh \xi }}{{2( {\cosh \xi - \cos \phi } )}} + {\alpha _{a00 }} + {\beta _{a00 }} + \sum_{m = 1}^{\infty} {[ {{\alpha _{am}}\cosh ( {m\xi } ) + {\beta _{am}}\sinh ( {m\xi } )} ]} \cos m\phi, \end{gather}

(3.45)\begin{align} {\tau _{a,\xi Z}} & = \frac{M}{\hbar }\frac{{\partial {u_{a,Z}}}}{{\partial \xi }} = \frac{M}{\hbar }\left[ { - \frac{{{f_a}{{\sinh }^{2}}{\xi _w}\sinh \xi \cos \phi }}{{2{{( {\cosh \xi - \cos \phi } )}^{2}}}}} \right] \nonumber\\ & \quad - \frac{M}{\hbar }\frac{1}{{{\xi _c} - {\xi _w}}}\left[ {\frac{{{f_a}{{\sinh }^{2}}{\xi _w}\cosh {\xi _c}}}{{2\sinh {\xi _c}}} - \frac{{{f_a}\sinh {\xi _w}\cosh {\xi _w}}}{2} - {u_{c,Z}}} \right]\nonumber\\ &\quad + \frac{M}{\hbar }\sum_{m = 1}^{\infty} {m [ {{\alpha _{am}}\sinh ( {m\xi } ) + {\beta _{am}}\cosh ( {m\xi } )} ]} \cos ( {m\phi } ), \end{align}

for which the coefficients are given in Appendix A.

3.3. Solution procedure for interface evolution

The annular and core fluid flux functions can be obtained as

(3.46)\begin{gather} {q_a} = \frac{1}{{\rm \pi} }\int_{{\xi _w}}^{{\xi _c}} {\int_0^{2{\rm \pi} } {{u_{a,Z}}{\hbar ^{2}}\, {{\rm d}} \phi\, {{\rm d}} \xi }}, \end{gather}

(3.47)\begin{gather}{q_c} = {u_{c,Z}}R_i^{2}, \end{gather}

which must satisfy the zero net flux condition, i.e.

(3.48)\begin{equation} q_a+q_c=0. \end{equation}

In this condition, the main unknown parameter of the flow is ${f_a} = ({1}/{M})({{\partial {P_a}}}/{{\partial Z}})$ and ${({{\partial {\tau _{c,ZZ}}}}/{{\partial Z}})}$ (or $u_{c,Z}$), which can be computationally obtained using established nested iterations for multilayer viscoplastic flows (Taghavi et al. Reference Taghavi, Seon, Martinez and Frigaard2009; Taghavi Reference Taghavi2018). The solution of the problem can be computed and it must satisfy the kinematic condition (Leal Reference Leal2007), i.e.

(3.49)\begin{equation} \frac{{\partial F}}{{\partial T}} + {\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla} F = 0, \end{equation}

which in combination with the continuity equation and some algebra translates to

(3.50)\begin{equation} {\frac{{\partial R_i^{2}}}{{\partial T}} + \frac{{\partial q_c}}{{\partial Z}} = 0}. \end{equation}

To find the interface motion versus time and space, the kinematic condition above (coupled to the other equations) is solved. To do so, an initially sharp interface is given at $Z = 0$, typically following a sharp linear function $R_i^{2} ( {Z,T = 0} ) = - a({{R_c^{2}}}/{2})Z + {{R_c^{2}}}/{2}$ (with the coefficient $a=10^{3}$). Let us recall that, here, ${R_c} = {{{{\hat R}_c}}}/{{\hat R}}$ is the inner radius of the injection pipe.

The discretization of the kinematic condition is implemented in conservative form, following a scheme that is first-order explicit in time and second-order in space, which for robust integration is coupled to a (shock-capturing) Van Leer flux limiter scheme (Yee, Warming & Harten Reference Yee, Warming and Harten1985). To calculate $u_{c,Z}$, an integration scheme is employed.

3.4. Yielding length

Now, it is straightforward to obtain the critical yielding length, $L_y$, which is the length at which the viscoplastic core layer starts to yield at the inlet ($Z=0$, where the stress is the highest). To do so, using (3.32) and the critical yield condition, i.e.

(3.51)\begin{equation} \frac{{\sqrt 3 }}{2}| {{\tau _{c,ZZ}}} | = B, \end{equation}

and knowing that $R_i=R_c$ and $u_{c,Z}=1$ hold before yielding, we find

(3.52)\begin{equation} {L_y} = \frac{{\sqrt 3 \delta B}}{{\displaystyle \left| { - \frac{{\partial {P_a}}}{{\partial Z}} - \frac{1}{{{\rm \pi} R_c^{2}}}\int_0^{2{\rm \pi} } {\hbar {\tau _{a,\xi_c Z}}} \, {{\rm d}} \phi + \chi } \right|}}. \end{equation}

Accordingly, the yielding time is $T_y=L$. Also, $f_a$ can be simply found by satisfying

(3.53)\begin{equation} \frac{1}{{\rm \pi} }\int_{{\xi _w}}^{{\xi _c}} {\int_0^{2{\rm \pi} } {{u_{a,Z}}{\hbar ^{2}}\, {{\rm d}} \phi\, {{\rm d}} \xi } } + R_c^{2} = 0, \end{equation}

and, then, the shear stress ${{\tau _{a,\xi _c Z}}}$ can be accordingly found using (3.45). Note that the sign of the denominator in (3.52) determines if the yielding occurs due to an extension or compression, i.e.

(3.54)\begin{equation} \left.\begin{array}{c@{}} \displaystyle {\rm sgn} \left( { - \dfrac{{\partial {P_a}}}{{\partial Z}} - \dfrac{1}{{{\rm \pi} R_c^{2}}}\int_0^{2{\rm \pi} } {\hbar {\tau _{a,\xi_c Z}}} \, {{\rm d}} \phi + \chi } \right) > 0\text{: yielding due to extension} \Rightarrow {\rm breakup};\\ \displaystyle {\rm sgn} \left( { - \dfrac{{\partial {P_a}}}{{\partial Z}} - \dfrac{1}{{{\rm \pi} R_c^{2}}}\int_0^{2{\rm \pi} } {\hbar {\tau _{a,\xi_c Z}}} \, {{\rm d}} \phi + \chi } \right) < 0\text{: yielding due to compression} \Rightarrow {\rm buckling}. \end{array}\right\} \end{equation}

Also, the critical transition between the buckling (bulging) and breakup regimes occurs when the denominator in (3.52) becomes zero, i.e.

(3.55)\begin{equation} - \frac{{\partial {P_a}}}{{\partial Z}} - \frac{1}{{{\rm \pi} R_c^{2}}}\int_0^{2{\rm \pi} } {\hbar {\tau _{a,\xi_c Z}}} \, {{\rm d}} \phi + \chi = 0, \end{equation}

which corresponds to an infinitely long viscoplastic core layer, which is in a perfect balance of forces with no yielding.

The model presented in this section so far delivers the transition conditions from an unyielded to a yielded viscoplastic core fluid, for infinitely long core fluid (e.g. infinite falling height). The yielding can be due to the extensional or compressional stresses, eventually leading to the breakup or buckling (bulging) of the viscoplastic core, respectively. However, if the applied stresses on the core fluid perfectly counterbalance one another, the core fluid does not buckle (bulge) nor break up, and it remains unyielded. On the other hand, considering a finite falling height ($H$), i.e. the case in our experiments, (3.52) delivers a yielding length ($L_y$) that can be compared with the falling height. Provided that the falling height is larger than the yielding length ($H>L_y/\delta$), the theory would suggest that the core fluid is yielded (due to extensional or compressional stresses depending on the parameters), leading to the breakup or buckling regimes. Conversely, if the falling height is smaller than the yielding length ($H\le L_y/\delta$), the theory would suggest that the core fluid remains radially unyielded over the falling height (implying that the core fluid front theoretically can reach the pipe end without bucking/bulging or breakup), which would roughly translate to the coiling regime in the context of our experiments.

Considering its limit case, our theoretical result for the yielding length is consistent with those obtained from previous theoretical approaches for viscoplastic fluid flows injected in air (i.e. a dynamically inactive fluid). In this case, of course, the yielding would be always due to elongation not compression. If the annular fluid was ignored (as would be the case for the injection into air), the yielding length of (3.52) would become

(3.56)\begin{equation} \frac{{{L_y}}}{\delta } = {\ell _y} = \frac{{{{\hat \ell }_y}}}{{\hat R}} = \frac{{\sqrt 3 B}}{\chi } = \frac{{\sqrt 3 {{\hat \tau }_y}}}{{( {{{\hat \rho }_c} - {{\hat \rho }_a}} )\hat g\hat R}}, \end{equation}

which for air (${{\hat \rho }_a} \ll {{\hat \rho }_c}$) becomes

(3.57)\begin{equation} {\hat{\ell} _y} \approx \frac{{\sqrt 3 {{\hat \tau }_y}}}{{{{\hat \rho }_c}\hat g}}. \end{equation}

The relation above is fully consistent with the yielding length obtained by Balmforth et al. (Reference Balmforth, Dubash and Slim2010b) and Coussot & Gaulard (Reference Coussot and Gaulard2005). In particular, the yielding length from Balmforth et al. (Reference Balmforth, Dubash and Slim2010b) can be transferred to our scaling as

(3.58)\begin{equation} {{\hat \ell }_y} = \frac{{\sqrt 3 {{\hat \tau }_y}}}{{{{\hat \rho }_c}\hat g}} + \frac{{\hat \varsigma }}{{{{\hat \rho }_c}\hat g{{\hat R}_c}}}, \end{equation}

in which $\hat \varsigma$ denotes the surface tension. Note that the latter term that includes the surface tension effects is neglected in our miscible flow (which is also at the limit of large $Pe$).

3.5. Specific case of concentric injection flows

For concentric injection flows ($E=0$), the bipolar coordinates $( {\xi,\phi,Z} )$ simplify to the cylindrical coordinates $(r,\theta,Z)$. In this case, fully analytical expressions can be obtained for the quantities of interest, such as the pressure gradient and the interfacial shear stress, respectively:

(3.59)\begin{gather} {f_a} = \dfrac{1}{M}\dfrac{{\partial {P_a}}}{{\partial Z}} = \dfrac{{ - 4{u_{c,Z}}}}{{\ln ( {{R_i}} )R_i^{2} - R_i^{2} + \ln ( {{R_i}} ) + 1}}, \end{gather}

(3.60)\begin{gather}{\tau _{a,rZ}} = \dfrac{{{u_{c,Z}}M( {1 - R_i^{2}} )}}{{( {\ln ( {{R_i}} )R_i^{2} - R_i^{2} + \ln ( {{R_i}} ) + 1} ){R_i}}}, \end{gather}

where ${\tau _{a,rZ}} = {{{\tau _{a,rZ}}} |_{r = {R_i}}}$. The critical yielding length, $L_y$, can be also obtained straightforwardly:

(3.61)\begin{equation} {L_y} = \frac{{\sqrt 3 \delta B}}{{\left| {\dfrac{{2M\left( {2 + \left( {\dfrac{1}{{R_c^{2}}} - 1} \right)} \right)}}{{ ( {\ln ( {{R_c}} )R_c^{2} - R_c^{2} + \ln ( {{R_c}} ) + 1} )}} + \chi } \right|}}, \end{equation}

in which the denominator determines if the yielding occurs due to an extension or compression, i.e.

(3.62)\begin{equation} \left.\begin{array}{c} {\rm sgn} \left({\dfrac{{2M\left( {2 + \left( {\dfrac{1}{{R_c^{2}}} - 1} \right)} \right)}}{{( {\ln ( {{R_c}} )R_c^{2} - R_c^{2} + \ln ( {{R_c}} ) + 1} )}} + \chi } \right) > 0\text{: yielding due to extension} \Rightarrow {\rm breakup};\\ {\rm sgn} \left({\dfrac{{2M\left( {2 + \left( {\dfrac{1}{{R_c^{2}}} - 1} \right)} \right)}}{{( {\ln ( {{R_c}} )R_c^{2} - R_c^{2} + \ln ( {{R_c}} ) + 1} )}} + \chi } \right) < 0\text{: yielding due to compression} \Rightarrow {\rm buckling}. \end{array} \right\} \end{equation}

4.1. Flow regimes

Let us begin with the transition boundaries between the breakup, coiling and buckling regimes. First of all, note that (3.52) gives a simple prediction to the onset of yielding of the viscoplastic core fluid, as a function of $R_c, \chi, B$, $M$, $\delta \equiv 1/H$ and $E$ (note that the onset of yielding is independent of the power-law index, $n$). Conservative critical boundaries between the regimes can be simply obtained by considering the critical yielding length to be equal to the falling height ($L_y/\delta =H$), which via rearranging (3.52) gives

(4.1)\begin{equation} \frac{{MH}}{B} = \frac{{\sqrt 3 }}{{\displaystyle \left| { - \frac{1}{M}\frac{{\partial {P_a}}}{{\partial Z}} - \frac{1}{{{\rm \pi} R_c^{2}M}}\int_0^{2{\rm \pi} } {\hbar {\tau _{a,\xi_c Z}}} \, {{\rm d}} \phi + \frac{\chi }{M}} \right|}},\end{equation}

where the first and second terms in the dominator, corresponding to the annular flow dynamics (i.e. the pressure gradient and shear stress terms), are directly calculated using the model. Thus, for given values of $E$ and $R_c$, the two relevant dimensionless groups of the flow are

(4.2)\begin{gather} \frac{\chi }{M} = \frac{{( {{{\hat \rho }_c} - {{\hat \rho }_a}} )\hat g{{\hat R}^{2}}}}{{{{\hat \mu }_a}{{\hat V}_0}}}, \end{gather}

(4.3)\begin{gather}\frac{{MH}}{B} = \frac{{{{\hat \mu }_a}{{\hat V}_0}\hat H}}{{{{\hat \tau }_y}{{\hat R}^{2}}}}. \end{gather}

Note that, for concentric injection flows ($E=0$), (4.1) reduces to a simplified analytical relation:

(4.4)\begin{equation} \frac{{MH}}{B} = \frac{{\sqrt 3 }}{{ \displaystyle \left| {\frac{{2\left( {2 + \left( {\frac{1}{{R_c^{2}}} - 1} \right)} \right)}}{{( {\ln ( {{R_c}} )R_c^{2} - R_c^{2} + \ln ( {{R_c}} ) + 1} )}} + \frac{\chi }{M}} \right|}}. \end{equation}

Before we proceed, let us clarify the discussion above on how to calculate the regime transition boundaries using the model. First of all, note that the model considers the yielding/non-yielding condition at the inlet of the core liquid since, prior to the yielding, the stresses are theoretically maximum at this location. From a theoretical perspective, this implies that, provided the core fluid does not yield at this location, it will not yield anywhere else in the pipe. Therefore, using the model, we can first calculate the transition boundaries between the core fluid flows which theoretically yield radially before reaching the pipe end and those which do not (using (4.1)). Once these transitions are identified, we can look at the sign of the denominator of (4.1) to theoretically categorize the yielding flows as breakup or buckling.

Figure 6 plots the theoretical critical boundaries via (4.1) in the plane of $MH/B$ and $\chi /M$ (for fixed $E=0$ and $R_c=0.5$); the breakup, coiling and buckling regimes are marked. Each regime corresponds to a different balance among the governing forces that control the motion of the viscoplastic core. The breakup regime, where the viscoplastic core is yielded due to the extensional stress overcoming the yield stress, is found at large values of $\chi / M$, corresponding to large buoyancy, low injection velocities and low viscosities of the annular fluid. The buckling regime, on the other hand, occurs at $\chi / M \lessapprox O(10^{2})$ and $M H/B\gtrapprox O(10^{-2})$, representing small buoyancy, small yield stresses as well as large injection velocities and high viscosities of the annular fluid. In this regime, the yielding of the viscoplastic core towards buckling is due to the compressional stress dominating the yield stress. The area between the critical boundaries of the breakup and buckling regimes marks the coiling regime, where forces counterbalance one another in a way such that, theoretically, the viscoplastic core fluid neither breaks up nor buckles but, instead, shows a coiling behaviour. It is evident, from the regime map, that large yield stresses (large $B$) withstand extensional or compressional stresses, resulting in the coiling regime.

Figure 6. Regime classification in terms of the breakup, buckling (bulging) and coiling regimes in the plane of $MH/B$ and $\chi / M$, for $E=0$ and $R_c=0.5$, where different regimes are marked by different colours.

There exist alternative approaches to define the appropriate dimensionless groups, especially when dealing with viscoplastic fluids. In particular, if instead of $( {{{\hat \mu }_c}{{\hat V}_0}} )/\hat R$ we use ${\hat \tau _y} + ( {{{\hat \mu }_c}{{\hat V}_0}} )/\hat R$ as the characteristic stress (Caliman, Soares & Thompson Reference Caliman, Soares and Thompson2017), the main dimensionless quantities change to

(4.5)\begin{equation} \left. \begin{array}{c@{}} {B^{*}} = B\dfrac{{( {{{\hat \mu }_c}{{\hat V}_0}} )/\hat R}}{{{{\hat \tau }_y} + ( {{{\hat \mu }_c}{{\hat V}_0}} )/\hat R}},\\ {M^{*}} = M\dfrac{{( {{{\hat \mu }_c}{{\hat V}_0}} )/\hat R}}{{{{\hat \tau }_y} + ( {{{\hat \mu }_c}{{\hat V}_0}} )/\hat R}},\\ {\chi ^{*}} = \chi \dfrac{{( {{{\hat \mu }_c}{{\hat V}_0}} )/\hat R}}{{{{\hat \tau }_y} + ( {{{\hat \mu }_c}{{\hat V}_0}} )/\hat R}}. \end{array} \right\} \end{equation}

Using these new definitions, the main quantities employed to present the main results (e.g. the regime classification of figure 6) would not change, i.e.

(4.6)\begin{gather} \frac{{{\chi ^{*}}}}{{{M^{*}}}} \equiv \frac{\chi }{M}, \end{gather}

(4.7)\begin{gather}\frac{{{M^{*}}H}}{{{B^{*}}}} \equiv \frac{{MH}}{B}. \end{gather}

Although the above-mentioned way of addressing the dimensionless analysis can be more aligned with the recent findings in yield stress materials, e.g. Kordalis et al. (Reference Kordalis, Varchanis, Ioannou, Dimakopoulos and Tsamopoulos2021) highlighting the advantages of using the plastic number ($B^{*}$) over the Bingham number ($B$), it does not change the outcomes of our work.

We have so far shown the existence of three main flow regimes, i.e. the breakup, coiling and buckling regimes, when a viscoplastic fluid is injected into a miscible medium (at the limit of large $Pe$), and we have accordingly delineated conservative transition boundaries between the flow regimes through the mathematical model, via the flow regime map in the plane of $MH/B$ and $\chi / M$ (figure 6). Now, it is worth mentioning that, in terms of (4.1), there exist at least two obvious limit cases, i.e. ${\chi }/{M}= {{( {{{\hat \rho }_c} - {{\hat \rho }_a}} )\hat g{{\hat R}^{2}}}}/{{{{\hat \mu }_a}{{\hat V}_0}}}\gg 1$ and ${\chi }/{M}={{( {{{\hat \rho }_c} - {{\hat \rho }_a}} )\hat g{{\hat R}^{2}}}}/{{{{\hat \mu }_a}{{\hat V}_0}}}\ll 1$, corresponding to the situations where the ambient fluid flow dynamics is either completely negligible or completely dominant, respectively. Regarding the former (${\chi }/{M}\gg 1$), the viscoplastic fluid flow can only transition between the breakup and coiling regimes and the critical transition boundary can be simply obtained as

(4.8)\begin{equation} \frac{H}{B} \approx \frac{{\sqrt 3 }}{\chi },\end{equation}

which relates, for example, the critical transition falling height ($H$) to the Bingham number ($B$) and the buoyancy number ($\chi$). In this scenario, as $B$ increases or $\chi$ decreases, the flow transitions from the breakup regime (where the viscoplastic core fluid breaks up) to the coiling regime (where the viscoplastic core fluid does not break up), while the ambient fluid flow dynamics does not play much of a role. From this perspective, therefore, the theory and the critical relation obtained can be directly used for much simpler flow configurations, e.g. the extrusion of a typical viscoplastic fluid into air (for which ${\chi }/{M}= {{( {{{\hat \rho }_c} - {{\hat \rho }_a}} )\hat g{{\hat R}^{2}}}}/{{{{\hat \mu }_a}{{\hat V}_0}}}\gg 1$, since ${\hat \rho }_a$ and ${\hat \mu }_a$ are both very small, implying that air is practically a dynamically inactive fluid), albeit with ignoring the surface tension effects. This finding is consistent with the spirit of previous works (e.g. Coussot & Gaulard Reference Coussot and Gaulard2005; Balmforth et al. Reference Balmforth, Dubash and Slim2010a,Reference Balmforth, Dubash and Slimb; Geffrault et al. Reference Geffrault, Bessaies-Bey, Roussel and Coussot2021), which have mainly focused on the breakup of extruded viscoplastic fluids into air, although they have not typically considered/investigated a coiling regime due to unboundedness of their flow domains. On the other hand, regarding the other limit in (4.1), i.e. ${\chi }/{M}\ll 1$, the viscoplastic fluid flow can only transition between the buckling and coiling regimes and the critical transition boundary can be simplified as

(4.9)\begin{equation} \frac{{MH}}{B} \approx c,\end{equation}

in which $c$ is a constant, i.e. only dependent on $E$ and $R_c$. Therefore, when the ambient fluid flow dynamics is dominant, the flow almost always belongs to the buckling regime, unless $B$ is very large or $H$ is very small.

4.2. Comparison with the experimental data

Figure 7 compares the three regimes based on the model and experiments. Despite the simplifying assumptions made in the model as well as complex experimental behaviours, the agreement between the model predictions and experiments is reasonable. However, some discrepancies are observed, e.g. the experimental coiling regime region is larger than that predicted by the model, a feature which may be justified. First, note that the critical condition of (4.1) is quite conservative and delivers only a prediction to the onset of yielding, which in the case of an extensional flow would not necessarily be equivalent to the breakup; in fact, when the core fluid yields near the inlet due to an extension, as the injection continues, the core fluid front can still evolve and reach the pipe end before the full breakup can occur. In this case, the coiling may occur for a core fluid that is partially yielded but not fully broken. Therefore, the experimental coiling (breakup) regime region is slightly larger (smaller) than that predicted by the model. On the other hand, yielding due to a compression almost certainly results in the buckling (bulging) regime because, in this case, the evolution of the core fluid front is significantly slowed down as the core fluid radially expands rapidly and cannot coil. This may justify why the model prediction to the transition between the buckling and coiling regimes may be more satisfactory.

Figure 7. Comparison between the regime classifications based on the model (solid lines) and experiments (symbols, representing $>300$ experiments) for $R_c=E=1/3$, in the plane of $MH/B$ and $\chi / M$, for $R_c=E=1/3$: the breakup regime ($\square$); the coiling regime ($\circ$); the buckling regime ($\vartriangle$). The filled circle symbols mark the cases where the core fluid yields and thins but its front eventually reaches the pipe end before breakup, and starts coiling.

As mentioned above, there are certain cases categorized within the coiling regime for which the core fluid is yielded but, due to the finiteness of the domain, its front reaches the pipe end before breaking up and starts coiling (see figure 4b). We have quantified and marked such flows, using the filled symbols, in the regime map of figure 7. Interestingly, such flows are all within the breakup regime zone predicted by the model, implying that the model can provide a good prediction to the yielding behaviour and these experiments would indeed exhibit a breakup for a longer falling height. We should also note that the coiling data points in figure 7 include different sorts of observed coiling behaviours, i.e. regular, irregular and free coiling (as discussed in § 2).

Before we proceed, let us summarize some possible sources/reasons behind the experimental–modelling discrepancy observed in our work (see figure 7). (i) One of the possible reasons for the deviation of the experimental data from the model predictions (e.g. near the boundary of the coiling and breakup regimes) could be due to the yield stress values (obtained based on the HB model) used in our analysis. As shown earlier in § 2, the yield stresses obtained from the HB model fitting (figure 2(a) and table 1) are slightly different from those obtained using other measurement methods (e.g. compare figures 2a, 2b, 2c and 2d). This difference highlights the fact that viscoplastic fluids are complex materials in terms of their yielding behaviour and that the yield stress measured using steady-state shear rheometry may not be the best choice across the board in different flow confirmations. (ii) Relevant to (i), a possible source of mismatch may be due to differences between extensional and viscometric yielding (which is the most widely used method to probe viscoplastic fluid rheology). For instance, Sica, de Souza Mendes & Thompson (Reference Sica, de Souza Mendes and Thompson2020) and Thompson & de Souza Mendes (Reference Thompson and de Souza Mendes2020) have considered extensional and compressional states of viscoplastic fluids extruded into air and they have quantified the corresponding yield stresses (namely static and dynamic yield stress values), which can be different from those delivered by steady-state shear rheometry. (iii) Another source of discrepancy between the model and experiments may be related to considering only normal components in the viscoplastic fluid flow and ignoring shear components in our asymptotic scaling analysis. Although in our work we may expect such an omission to have only negligible effects, there also are works in the literature (e.g. Moschopoulos et al. Reference Moschopoulos, Syrakos, Dimakopoulos and Tsamopoulos2020; Geffrault et al. Reference Geffrault, Bessaies-Bey, Roussel and Coussot2021) arguing that neglecting shear components of the strain rate tensor may lead to inaccuracies of slender filament approximations. (iv) Regarding the breakup regime, the discrepancy can also originate from the notion of yielding and rupturing in viscoplastic fluids (as also mentioned further above). Although the model may provide a reasonable understanding and prediction to the onset of yielding, the final rupturing of our viscoplastic fluid is hard to model and fully analyse, especially considering the multiphase nature of the flow. However, there are studies in the literature for simpler flow configurations (Hu et al. Reference Hu, Bian, Grotberg, Filoche, White, Takayama and Grotberg2015; Tarcha et al. Reference Tarcha, Forte, Soares and Thompson2015; Thompson & de Souza Mendes Reference Thompson and de Souza Mendes2020; Bahrani et al. Reference Bahrani, Hamidouche, Moazzen, Seck, Duc, Muradoglu, Grotberg and Romanò2021) which have investigated the rupturing of viscoplastic fluids but in different conditions; for instance, Tarcha et al. (Reference Tarcha, Forte, Soares and Thompson2015) have analysed the yielding and rupturing of viscoplastic fluids using rheological tests, revealing how such fluids experience a destruction process and an eventual rupturing after the yield stress is overcome. (v) Finally, another discrepancy source is related to the nature of the model and experiment in our work. While the former is highly simplified, the latter is a complex flow, many features of which are ignored in the model for simplicity. For instance, the model does not obviously consider the coiling regime details (including regular, free and irregular coiling behaviours). In free coiling, in particular, the coiling starts before the core front touches the pipe end, i.e. a feature that is not included in the model and must be studied in future works.

A closer look at (3.52) reveals that, at the limit of a zero-injection velocity, the dimensional length of the viscoplastic core fluid at the yielding onset from the model ($\ell _y$) can be approximated as

(4.10)\begin{equation} {\left. {\frac{{{{\hat \ell }_y}}}{{\hat R}}} \right|_{{{\hat V}_0} \to 0}} \approx \frac{{\sqrt 3 B}}{\chi } = \frac{{\sqrt 3 {{\hat \tau }_y}}}{{( {{{\hat \rho }_c} - {{\hat \rho }_a}} )\hat g\hat R}} \Rightarrow { {{{\hat \ell }_y}} |_{{{\hat V}_0} \to 0}} \approx \frac{{\sqrt 3 {{\hat \tau }_y}}}{{( {{{\hat \rho }_c} - {{\hat \rho }_a}} )\hat g}}, \end{equation}

since the terms regarding the pressure gradient and the shear stress applied by the annular fluid flow vanish in (3.52) as $\hat {V}_0 \to 0$. Equation (4.10) can also provide a way to approximate the yield stress, provided that the yielding length is known at the limit of zero-injection rates. This approach has been previously demonstrated and used by Sica et al. (Reference Sica, de Souza Mendes and Thompson2020) and Thompson & de Souza Mendes (Reference Thompson and de Souza Mendes2020) in conceptually similar contexts. Figure 8(a) plots the length of the viscoplastic core at the yielding onset from the experiments ($\hat {\ell }^{experiment}_y$) for small injection velocities. In this figure, the symbols mark Sample III (circles) and Sample IV (diamonds). As can be seen, $\hat {\ell }^{experiment}_y$ for Sample III is larger than that for Sample IV. Following Sica et al. (Reference Sica, de Souza Mendes and Thompson2020) and Thompson & de Souza Mendes (Reference Thompson and de Souza Mendes2020), it is possible to use a linear least-squares fit as a crude way to estimate the yielding length at the limit of the zero-injection velocity, namely where the fitted line crosses the vertical axis, i.e. ${{\hat \ell _y^{\textrm {{experiment}}}} |_{{{\hat V}_0} \to 0}}$. Accordingly, using (4.10), the associated yield stresses for Samples III and IV can be approximated as 6.47 and 4.31 Pa, respectively, which in comparison with the rheometry results (table 1) show a difference of less than 7 %.

Figure 8. (a) Viscoplastic core length at the yielding onset $\hat {\ell }^{experiment}_y$ versus the injection velocity $\hat V_0$ for an injection flow in an eccentric pipe ($\hat {E}=6.35$ mm). The symbols mark Sample III (circle) and Sample IV (diamond). The lines mark the least-squares fitting of the data points. Length ${{\hat \ell _y^{\textrm {{experiment}}}} |_{{{\hat V}_0} \to 0}}$ is 0.115 and 0.039 m for Samples III and IV, respectively. (b) Dimensionless elapsed time between the yielding and rupturing $t_{r}-t_y$ versus $\chi /M$ for the breakup experiments for $E=1/3$. The symbol colours represent the values of the Bingham number.

For the sake of completeness of the result presentation, let us show the dimensionless elapsed time between the yielding onset (${t_y} = {{{{\hat t}_y}}}/{{( {\hat R/{{\hat V}_0}} )}}$) and the eventual of breakup or rupturing (${t_r} = {{{{\hat t}_r}}}/{{( {\hat R/{{\hat V}_0}} )}}$) of the core viscoplastic fluid, as shown in figure 8(b) for our experimental results. The yielding and rupturing times are extracted from the experimental data using an image processing technique. This figure shows $t_r-t_y$ versus $\chi /M$, while the values of the Bingham number are superimposed by the symbol colours. The results show that, by increasing $\chi /M$ and $B$, the elapsed time decreases, i.e. the core fluid ruptures faster. The results for the elapsed time between yielding and rupturing, such as those presented in figure 8(b), could be perhaps further analysed in future works to distinguish between two critical stresses, i.e. associated with yielding and rupturing.

4.3. Effect of $R_c$

It is insightful to study the effects of the variation of $R_c$ on the regime boundaries, as shown in figure 9. As seen, varying $R_c$ can significantly change the regime boundary positions, albeit in a non-monotonic way; by increasing $R_c$ up to 0.3, the critical boundaries are shifted to the left but, for higher values of $R_c$, they are moved to the right. In the latter case, the buckling regime region becomes significantly wider.

Figure 9. Critical regime transition boundaries in the plane of $MH/B$ and $\chi / M$, for $E=0$ and $R_c=0.1$ (dashed line), 0.3 (dash-dotted line), 0.5 (red solid line), 0.7 (dotted line) and 0.9 (black solid line).

The model predicts a non-monotonic effect of $R_c$ on the critical regime boundaries, especially in terms of the transition between the buckling and breakup regimes. In an attempt to explain this non-monotonic behaviour, let us consider a long core fluid, for which $R_c$ increases from small values, while the other parameters are fixed. On the one hand, note that the buoyancy force is the only force that pulls the core fluid downward, while the other forces (including the yield stress force and the shear stress force at the interface) resist against this force and, therefore, against the breakup. On the other hand, the buoyancy force per unit volume is constant (proportional to $\chi$), while the interfacial shear force per unit volume is proportional to $({1}/{{{R_c}}})( {{{{\tau _{a,rZ}}}}/{M}} )$. Therefore, initially at small $R_c$, the interfacial shear force per unit volume decreases with increasing $R_c$. However, as $R_c$ keeps increasing and the annular flow cross-section area decreases, the shear stress applied at the interface increases due to the flux constraint. This results in the interfacial shear force per unit volume eventually increasing at larger $R_c$, implying that a larger $\chi$ would be needed to overcome this force and lead to the transition from the buckling regime to the breakup regime.

To better illustrate the discussion above, figure 10 presents the variation of three parameters versus the core fluid radius. To simplify the analysis, the condition before yielding ($R_i=R_c$ and $u_{c,Z}=1$) for a concentric flow ($E=0$) is considered. The solid black line marks the non-monotonic variation of $| {({1}/{{{R_c}}})( {{{{\tau _{a,rZ}}}}/{M}} )} |$ versus $R_c$, i.e. initially decreasing and then increasing as $R_c$ grows. The variation of $({1}/{M})({{\partial {P_a}}}/{{\partial Z}})$ versus $R_c$ (the red dashed line) shows that, as may be expected, the pressure gradient increases monotonically as $R_c$ increases. Also, at large $R_c$, the value of $({1}/{M})({{\partial {P_a}}}/{{\partial Z}})$ progressively increases and dominates that of $| {({1}/{{{R_c}}})( {{{{\tau _{a,rZ}}}}/{M}} )} |$. Finally, the variation of ${( {{\chi }/{M}} )_{{{critical}}}}$ versus $R_c$ is also shown. For an infinitely long viscoplastic core layer, ${({{\chi }/{M}})_{{{critical}}}}$ delineates the critical transition between the buckling and breakup regimes (i.e. simply found by setting the denominator of (3.62) equal to zero). A final note inferred from figure 10 is that, as $R_c$ approaches the value of zero or one, $| {({1}/{{{R_c}})}( {{{{\tau _{a,rZ}}}}/{M}} )} |$ and ${( {{\chi }/{M}})_{{{critical}}}}$ approach large values, implying that the model may not be valid as $R_c$ approaches its extreme values.

Figure 10. Variation of $| {({1}/{{{R_c}}})( {{{{\tau _{a,rZ}}}}/{M}} )} |$ (solid black line), $({1}/{M})({{\partial {P_a}}}/{{\partial Z}})$ (red dashed line) and ${( {{\chi }/{M}} )_{{{critical}}}}$ (blue dash-dotted line) versus $R_c$ for $E=0$.

Figure 11(a) shows the regime classification in the plane of $R_c$ and $\chi /M$ for an infinitely long pipe. The flow configuration is concentric (i.e. $E=0$). The black solid line delineates the transition between the breakup and buckling regimes; there is, of course, no coiling regime and the line itself represents a core fluid layer that never breaks up or buckles in an infinitely long pipe. Figure 11(b) shows the regime transition boundaries for different values of ${{M{H}}}/{{B}}$. For ${{M{H}}}/{{B}}=\infty$, there is no coiling regime but, as ${{M{H}}}/{{B}}$ decreases, the coiling regime region appears and expands in the area sandwiched between the buckling and breakup regimes (corresponding to the dashed and dash-dotted lines). Moreover, for the ranges observed, the transition boundary has a minimum value of $\chi /M$ for a core with $R_c \approx 0.3$. Figure 11(c) depicts the regime classification for a typical finite value of ${{M{H}}}/{{B}}=0.01$. As seen, for a given value of $R_c$, on increasing $\chi /M$, the flow transitions from the buckling to coiling and eventually to the breakup regime. Also, interestingly, for $0.10 \lessapprox R_c \lessapprox 0.64$ and all the values of $\chi /M\ge 1$, the viscoplastic core fluid does not buckle.

Figure 11. Regime classification for an infinitely long pipe for $E=0$. (a) The regime classification of breakup and buckling for an infinitely long pipe. The thick solid line marks the critical transition. (b) Contour values of ${{M{H}}}/{{B}}$ for $E=0$. The thick solid line marks ${{M{H}}}/{{B}}=\infty$. The dash-dotted lines correspond to yielding (breakup) due to extension, where ${{M{H}}}/{{B}} = 0.1$, 0.03 and 0.01, from left to right. The dashed lines correspond to yielding (buckling) due to compression, where ${{M{H}}}/{{B}} = 0.1$, 0.03 and 0.01, from right to left. (c) Regime classification of the breakup, coiling and buckling regimes, in a pipe with the typical length of $H = (1/100)(B/M)$.

4.4. Effect of $E$

Figure 12(a–d) shows the regime classification in the plane of $R_c$ and $\chi /M$, for an infinitely long pipe and different eccentricities. The black solid line delineates the transition between the breakup regime (on the right-hand side) and the buckling regime (on the left-hand side). In an infinitely long pipe, there is no coiling regime and each line for different $E$ represents a core fluid layer that never breaks up or buckles. It should be noted that, with increasing eccentricity, the variation of $R_c$ becomes limited. It can be seen that, with increasing $E$, the buckling regime region is expanded. Figures 12(e) and 12(f) show observations similar to those in figures 11(b) and 11(c), respectively, but for an eccentric flow ($E=0.4$).

Figure 12. Regime classification for an infinitely long pipe for various eccentricities: (a) $E=0$, (b) $E=0.2$, (c) $E=0.4$ and (d) $E=0.8$. (e) Contour values of ${{M{H}}}/{{B}}$ for $E=0.4$. The thin solid line marks ${{M{H}}}/{{B}}=\infty$. The dash-dotted lines correspond to yielding (breakup) due to extension, where ${{M{H}}}/{{B}} = 0.1$, 0.03 and 0.01, from left to right. The dashed lines correspond to yielding (buckling) due to compression, where ${{M{H}}}/{{B}} = 0.1$, 0.03 and 0.01, from right to left. (f) Regime classification of the breakup, coiling and buckling regimes in a pipe with typical length of $H = (1/{20})({B}/{M})$.

Figure 13 illustrates a comparison between the results obtained for the two different values of $E$. The lines and data points represent the modelling and experimental results, respectively. For $R_c=1/3$, for the value of which the experimental results are also superimposed, the coiling regime is seen for smaller values of $\chi /M$ when $E=0$. In addition, the comparison between the experimental and modelling results shows reasonable agreement.

Figure 13. Regime classification based on the model (lines) and experiments (symbols) for (a) $E=1/3$, $H=(1/20)({B}/{M})$ and (b) $E=0$, $H=(3/100)({B}/{M})$. The bounded regions via lines show the breakup, coiling and buckling regimes from right to left, respectively. Different symbol shapes represent the breakup regime ($\square$), the coiling regime ($\circ$) and the buckling regime ($\vartriangle$) in the experiments.

4.5. Interface evolution and velocity field

Figure 14 depicts typical examples of the evaluation of a core fluid interface, at a given time $T$, in the absence and presence of the annular flow dynamics, for a concentric flow configuration. As seen, when the core fluids are yielded, the stress contours exceed $B$ everywhere, except for the core fluid's lowest part. Figures 14(a) and 14(c) show that, without the annular flow, the extensional stress (positive $\tau _{c,ZZ}$) dominates the yield stress in the thinning region, as the core progresses towards the breakup. When the annular flow dynamics is considered, and depending on the values of $\chi$, $B$ and $M$ (for fixed $\delta$, $E$, $R_c$ and $n$), figures 14(b) and 14(d) reveal that the flow can transition between the breakup and buckling regimes. For instance, for small $M$, the extensional stress (positive $\tau _{c,ZZ}$) is dominant, while for large $M$ the compressional stress (negative $\tau _{c,ZZ}$) becomes dominant.

Figure 14. Interface evolution for $\delta =0.05$, $n=0.5$, $R_c=0.5$ and $E=0$ (a,c) without and (b,d) with considering the annular flow dynamics: (a,b) $B=\chi =1$ and $M=0.005$; (c,d) $B=\chi =0.1$ and $M=0.02$; (a) $T=0.19$, (b) $T=0.5$, (c) $T=0.6$ and (d) $T=0.24$. The superimposed colours and the colour bars represent the values of $\tau _{c,ZZ}$.

Figure 15 presents the interface shape, along with the dominant stress contours (figure 15a,c) and the axial velocity contours (figure 15b,d), for compressional and extensional flows of an eccentric viscoplastic core flow ($E=1/3$), at a given time $T$. Figures 15(a,b) and 15(c,d) correspond to the breakup and buckling regimes, respectively. In figures 15(a,c) and 15(b,d), the values of $\tau _{c,ZZ}$ and $u_{c,Z}$ are superimposed on the interface shape, respectively. Additionally, in figure 15(b,d), the annular flow velocity field near the necking/buckling region is also superimposed. As seen, increasing the viscosity of the annular fluid causes the extensional stress (positive $\tau _{c,ZZ}$) to become a compressional stress (negative $\tau _{c,ZZ}$); this leads to the transition of the flow dynamics from the breakup to the buckling regime. In both regimes, however, the unyielded region ($| {{\tau _{c,ZZ}}} |\le B$) is located near the core end, with lower velocity for the buckled core. Figure 15 also shows that, in the buckling regime, the core fluid layer grows much slower (note that all the results are at $T=0.3$). Also, for the breakup flow, the core velocity has the lowest value at the inlet, while for the buckling flow the core fluid velocity is highest near the inlet. In addition, the local velocity of the annular fluid is larger in the wide side of the eccentric core than that in the narrow side. As the eccentricity promotes the flow along the wide side, it affects on the dynamics of the injected viscoplastic core fluid, through changing the interfacial stresses due to the upward motion of the annular fluid. Finally, our results show that the predictions of the interface evolution using the model are in qualitative agreement with the early-time experimental interface shapes in the breakup and bucking regimes; see figure 1($b_1$–$b_3$).

Figure 15. Model results of an eccentric injection, for $\delta =0.05$, $n=0.5$, $R_c=E=1/3$ and $T=0.3$: (a,b) $B=\chi =1$ and $M=0.005$; (c,d) $B=\chi =0.1$ and $M=0.02$. (a,c) Interface shape, along with $\tau _{c,ZZ}$ superimposed as colours. (b,d) Interface shape, along with axial velocity ($u_{Z}$) superimposed as coluors. The core fluid velocity is given throughout the viscoplastic layer length while the annular velocity is given at $Z\approx L/3$.

In order to gain more insight into the flow, here, let us further analyse the core fluid behaviour without and with considering the annular fluid flow dynamics, as shown in figures 16 and 17, respectively. The former (i.e. when the annular fluid flow dynamics is ignored) would resemble the simpler case of the injection of a viscoplastic fluid into a dynamically inactive fluid, such as air (Coussot & Gaulard Reference Coussot and Gaulard2005; Balmforth et al. Reference Balmforth, Dubash and Slim2010a,Reference Balmforth, Dubash and Slimb; Geffrault et al. Reference Geffrault, Bessaies-Bey, Roussel and Coussot2021), albeit with ignoring surface tension effects. Figure 16(a) shows the variation of the interface radius along the length of the core versus the time. At long times, the core fluid yields after advancing some distance with respect to the injection inlet (when, loosely speaking, a sufficient weight of the core fluid is hung below the injection inlet). Meanwhile, as time progresses, the core fluid becomes thinner at the yielding point. The core is yielded, except for a small region towards its front. Figure 16(b) shows the growth of the core fluid length ($L$) versus time ($T$), which occurs due to the continuous injection of the viscoplastic fluid. Initially, before yielding, the layer length varies linearly with time. After yielding, however, $L(T)$ deviates from the linear line of $L(T)=T$, and grows more rapidly as expected. A sequence of profiles illustrating $\tau _{c,ZZ}$ along the core length are shown in figure 16(c), for the corresponding times in figures 16(a) and 16(b). The horizontal dash-dotted line in figure 16(c) represents the value of $B$. At early times, the magnitude of $\tau _{c,ZZ}$ along the core length shows a gradual increase with time. Before yielding ($\tau _{c,ZZ} \le B$), the maximum stress is at $Z=0$, which represents the injection inlet. Nevertheless, after yielding ($\tau _{c,ZZ} > B$), the point where the stress is maximum moves downward as time increases. Finally, considering the observations made in figure 16, very similar behaviours are observed when buoyancy is large while the annular flow dynamics is considered.

Figure 16. Model results for the dynamics of a breakup flow, for a case without considering the annular flow dynamics. The flow parameters are $\delta =0.05$, $n=0.5$, $B=\chi =1$, $R_c=0.5$ and $T_{end}=0.19$. (a) Interface radius $R_i$ versus $Z$ and $T$. The results are shown for $T=0,T_{end}/10,\ldots, 9T_{end}/10,T_{end}$. The last profile is marked by the red line, for which the shaded area marks the unyielded region. (b) The viscoplastic fluid layer length $L(T)$ versus time $T$. The dashed line is for comparison purposes and marks $L(T)=T$. (c) Plot of $\tau _{c,ZZ}$ versus $Z$ and $T$. The solid arrow shows the increase in time $T$. The results are shown for $T=0, T_{end}/10,\ldots, 9T_{end}/10, T_{end}$. The horizontal dash-dotted line marks $\tau _{c,ZZ}=B$.

Figure 17. Model results for the dynamics of a buckling flow, while considering the annular flow dynamics. The flow parameters are $\delta =0.05$, $n=0.5$, $B=0.1$, $\chi =0$, $R_c=0.5$, $M=0.02$, $E=0$ and $T_{end}=0.22$. (a) Interface radius $R_i$ versus $Z$ and $T$. The results are shown for $T=0, T_{end}/10,\ldots, 9T_{end}/10, T_{end}$. The last profile is marked by the red line, for which the shaded area marks the unyielded region. (b) The viscoplastic fluid layer length $L(T)$ versus time $T$. The dashed line is for comparison and shows $L(T)=T$. (c) Plot of $\tau _{c,ZZ}$ versus $Z$ and $T$. The solid arrow shows the increase in time $T$. The results are shown for $T=0, T_{end}/10,\ldots, 9T_{end}/10, T_{end}$. The horizontal dash-dotted line marks $\tau _{c,ZZ}=-B$.

Figure 17(a) shows the variation of the interface radius along the core length for different times in the presence of the annular flow. It can be observed that the inclusion of the annular flow creates sufficient viscous stresses for yielding and buckling of the core fluid, as a yielded region appears adjacent to the point of injection, where the compressive stress overcomes the yield stress. While the core fluid buckles radially, we can observe that the core fluid front section remains unyielded. As seen in figure 17(b), at early times, the growth rate of the viscoplastic layer length is constant, until the time at which the yielding and buckling of the core occur. After yielding, however, $L(T)$ deviates from the linear line of $L(T)=T$, and grows more slowly as expected. The variation of the stresses along the core in each time step is depicted in figure 17(c), showing a gradual increase in the stress magnitude with time. Since $B$ is small, $\tau _{c,ZZ}$ quickly overcomes $B$ and core fluid yields radially. Comparing figures 17(c) and 16(c) shows that the signs of $\tau _{c,ZZ}$ are opposite in these figures, which is trivial. However, comparing figures 16(c) and 17(c) reveals that the maximum value of the stress remains always close to $Z=0$, even after yielding (at least for the case considered).

The results presented in this section on the interface dynamics, obtained using the model, reveal that our theoretical predictions of the interface shape and its evolution are in agreement, albeit rather very qualitative, with the shape of the thinning region in the extensional flow or the buckling region in the compressional flow, e.g. as observed in our experiments; see e.g. figure 3.