Stability and breakup of liquid threads and annular layers in a corrugated tube

We consider the axisymmetric evolution of a viscous core liquid thread surrounded by another immiscible fluid inside a cylindrical tube of (dimensional) radius $b(z)$ that varies in the axial direction with corrugation. See figure 1. The interface is given by $r=S(z,t)$ with constant surface tension $\gamma$ and the core region is occupied by fluid 1 with density $\rho_{1}$ and viscosity $\mu_{1}$, and the annular region is filled by fluid 2 with $\rho_{2},\mu_{2}$. The flow fields are assumed to be axisymmetric with velocity components $(u,0,w)$ in terms of cylindrical coordinates $(r,\theta,z)$.

The dimensionless equations can be made if lengths are rescaled by the undisturbed core thread radius $a$, velocities by $\gamma/\mu_{1}$, time by $a\mu_{1}/\gamma$, pressure by $\gamma/a$. The subscripts $i=1,2$ are hereinafter used to denote core and annular fluids respectively. The difference from the problem studied in [74] is that the tube wall here has cross section variations, $d=d(z)=b(z)/a$ and for simplicity, the sinusoidal shape is considered in present work. Thus we have the governing equations

where subscripts denoting derivatives, $\chi_{1}=1,\chi_{2}=\chi=\rho_{2}/\rho_{1},R_{e}=\rho_{1}\gamma a/\mu^{2}_{1}$ and $\lambda_{1}=1,\lambda_{2}=\lambda=\mu_{2}/\mu_{1}$ , with boundary conditions at tube wall $r=d=b/a$ and at interface $r=S$ respectively

Notice that our definition of viscosity ratio is the same as in [36], [51], [60], [74], where the small $\lambda$ corresponds to a case with a very viscous core and a relatively less viscous annulus, while in some other studies, such as [69], the viscosity ratio may mean the opposite case. To ease the discussion, we stick with our definition.

The stress balances at the fluid interface are given as

where $[()]_{1}^{2}$ denotes the jump across the interface, that is, $()_{2}-()_{1}$. Finally, there is a kinematic boundary condition for the interface position $S$ that is given by

The aim here is to capture the physics with a set of simplified equations and we perform an asymptotic reduction of the above equations by assuming the axial length scale is much larger than the radial one with slenderness parameter $\epsilon=a/\mathcal{L}\ll 1$. Then we have the transformation $\partial_{z}\rightarrow\epsilon\partial_{z^{\prime}}$ where $z^{\prime}$ is the axial variable in long wave regime. In the rest of the article, the prime is dropped to simplify the notation. We shall also assume a weak topography with $\epsilon d_{z}\ll 1$ so that the dimensional slope of the wall is small. Following [58] and [74], we introduce the following expansion

In the fluid region, in order to take both core and annular fluids into account, we take the annular fluid to be less viscous, $\lambda=\epsilon^{2}\lambda_{0}$ (see also the derivation in [74]). This is different from the usual long wave model in [24], [58] because here we also want to include the effect of tube wall by capturing the interaction between the tube wall and annular fluids. When there is no outer wall confinement, a similar limit $\lambda\ll 1$, where the core fluid thread is much more viscous than the surrounding medium, has been considered in the asymptotic modeling in [51, 7]. We further ignore the inertial effect in the annular region, which amounts to assuming $\chi R_{e}<O(\epsilon^{2})$. This allows us to make analytic progress with the model. From now on, we drop the superscript $0$ but retain $1$ to emphasize the high order term. The axial momentum equation then gives

whose solutions are readily to be obtained, by using continuity equation

where $A,\ B,\ C$ are functions of $z$ and $t$.

In the core, $w_{1}=w_{1}(z,t)$ in the leading order, which is independent of radial coordinate. Using continuity equation, $u_{1}=-rw_{1z}/2$ as in the classical single fluid jet problem [24, 58]. Before considering stress balances at interface, we eliminate $B$ and $C$ by using no-slip and no-penetration boundary conditions at wall $d=d(z)$ to obtain

where $d_{z}=0$ corresponds to the case in [74]. Combining with the continuous velocity condition along $r=S$, $w_{1}=w_{2}$ and $u_{1}=u_{2}$, after some algebra we arrive at the following relation between $A$ and $p_{2}$,

where $f(t)$ is the quasi one dimensional force in the liquid thread (see also [64, 65, 57, 58]). Therefore

On the interface, the tangential and normal stress balances in leading order, read as

where in normal stress balance, $\epsilon^{2}S_{zz}$ is retained as in [24], to obtain the short wave cutoff in the linear stability region. The higher order $w^{1}_{1r}$ is obtained from the axial momentum equation,

where $\bar{R}_{e}=R_{e}/\epsilon^{2}$ (so to ensure inertia does not enter the leading order in the annular region, $\chi\sim o(1)$ is needed; so we have assumed that the highly viscous fluid is also more dense than the less viscous one in our long wave theory). Substituting (16) into (17) and using (18) leads to

where $\kappa=1/S-\epsilon^{2}S_{zz}$. With inviscid gas or vacuum instead of the viscous annular fluids, $A=0$ as in the single jet model [24]. Using the above results and kinematic condition (8) on interface, we reach the following equations for interface position $S$ and axial velocity $w_{1}$, namely,

where the subscript $1$ in $w$ is dropped and

In general, $f(t)$ is unknown and determined by some global constraint.

When $\lambda_{0}=0$ and $f(t)=0$, the model recovers the one in [24] for a single fluid jet. In current study, we assume spatial periodic solution on $[-L,L]$, which gives

To ease the discussion, for the rest of the paper, we further impose an even thread interface profile with respect to $z=0$ (also for $d(z)$), then $w$ is an odd function on $[-L,L]$ from (21). This is consistent with the stability analysis in literature in the sense that the analysis starts with a sinusoidal perturbation on the thread interface. Subsequently $f(t)\equiv 0$, which is used in the rest of study here and we leave $f(t)\neq 0$ case to future exploration. In addition, we notice that with $f(t)=0$ in (22), the model equations (21) and (22) are formally the same as the equations in [74], even though the wall shape effect has been included in the present study. Meanwhile, neglecting the wall variation, the two phase slender jet model in [51] is recovered formally by taking $\bar{R}_{e}=0$ and rescaling $\lambda_{0}=\ln d\bar{\lambda}$ (namely, $\lambda=\epsilon^{2}\ln d\bar{\lambda}$) with $\bar{\lambda}$ order one quantity as $d\rightarrow\infty$.

Before examining the fully nonlinear model, we bring up one interesting limit that many authors have investigated, the thin annular film limit [33, 59, 75, 76]. Following [75], we write

where $\delta\ll 1$, $\phi(z)$ is the wall shape function, $h(z,t)$ stands for the annular film thickness and $\sigma^{\prime}=0$ corresponds to a uncorrugated wall or a straight tube case [33, 59]. At the fluid interface, substituting (25) into (23), one arrives at

Rescaling time as $\tau=\delta^{3}t$, the leading order equation for $h$ is obtained as a modified Hammand equation

where Hammond equation (see [33]) can be recovered by taking $\sigma^{\prime}=0$, $\epsilon=1$ and adsorb the factors $12\lambda_{0}$ into time. A further simplified linearized equation can be obtained by taking $h=1+\alpha\xi$ in (27). After some algebra, the first few terms in equation are obtained as

Upon rescaling in time variable, this equation is the same as the linear equation in [75] in the strong surface tension limit if $\alpha\ll 1$ (see their (5.3)) and recovers the weakly nonlinear one in [76] for $\alpha\sim O(1)$ (see their (4.2)). In the derivation in [75, 76], the viscosity ratio is $\lambda\sim O(1)$ while in our long wave model, $\lambda\sim O(\epsilon^{2})$. Thus it seems that the different viscosity contrast (at least for $\lambda\leq O(1)$) only affects the time scale factor in the thin film equation.

We will not repeat the results in [75] in detail, but merely use their results to validate our method in computing the spectrum later in the linear stability analysis. We only focus on the long wave model (21), (22) in the rest of the discussion.

It is easy to see $\bar{S}=1$ and $\bar{w}=0$ still serve as the base state in our problem without any base flow or external force field. For simplicity and based on the observation in [74] that the inertia contributes little to the dynamics, only the inertialess case is studied for the linear stability and we probe the effect of inertia later by direct numerical simulations. We set $S=1+S^{\prime}$ and $w=0+w^{\prime}$, where $|S^{\prime}|\ll 1$, $|w^{\prime}|\ll 1$. For brevity, the prime is suppressed in the following analysis. The linearized system is reduced as following in the Stokes limit ($\bar{R}_{e}=0$),

Meanwhile, we consider the wall’s variation as a sinusoidal profile, $d=d_{0}+\sigma\cos(k_{w}z)$ where $k_{w}$ is the wavenumber of the tube wall. Recall we have assumed a weak topography, $\epsilon d_{z}\ll 1$, which leads to the condition $\epsilon\sigma k_{w}\ll 1$, i.e. $k_{w}\ll\left(\epsilon\sigma\right)^{-1}$. So the dimensionless wall wavenumber need not to be small provided $\sigma$ is not large (see [67] for similar argument in a related problem). For a straight tube, $\sigma=0$, then $G(1,d)$ is a constant and the growth rate is obtained by taking Fourier transform on (29) directly as in [74]. The dispersion relation is then given by

When the wall has cross section variation, $\sigma\neq 0$, the coefficient in (29) is a periodic function, which requires Floquet-Bloch theory to analyze the stability problem (see related studies in [56, 75, 3, 5]).

We employ the FFH method [20] for the eigenvalue problem here, which is a more straightforward numerical method than the ones in [75]. We have first used this method to reproduce the results in [75] and the convergence can be achieved quickly by using around $10\sim 20$ Fourier modes while [75] reported that a relatively large linear system is needed to solve the eigenvalue problem accurately. The details of FFH for our problem can be found in the Appendix.

To calculate the spectrum numerically, we choose a cut-off $N$ on the number of Fourier modes of the eigenfunctions $\psi$ (see the Appendix), resulting in a linear system of dimension $2N+1$. As mentioned earlier, we reproduced the results in [75] with around $N=5\sim 10$ (convergence within one period is achieved thus) and similar in the current linear problem. One advantage of FFH method is that one can handle any periodic function (in coefficient) in a general way by taking the Fourier transform (numerically), rather than simple cosine or sine function that can be easily incorporated in analysis. In current paper, we only focus on the simple sinusoidal profile though, which already gives a fairly complicated coefficient function $g$ in (37).

We show typical results in figure 2, where the wall shape wavenumber varies as indicated in figure. For comparison, the dotted lines in figure represent the growth curve in the case of uncorrugated tube, $\sigma=0$. Unlike the presentation in [75], where the primary and secondary branches are identified, only the dominant growth rate or the largest eigenvalue is shown in figure. Those results are in qualitatively good agreement with [75]. When the wall shape has relatively long wave length variation i.e. small $k_{w}$, panel (a) shows that both the long and short waves are excited due to the presence of wall corrugation which leads to larger growth rate than the uncorrugated case for small and large $k$. For example, when $k=1,k_{w}=6$, the dominant growth rate will be associated with the first wall harmonic $|k-k_{w}|=5$, which grows faster than $k=1$ in the uncorrugated or straight tube case. For intermediate range of $k$, the growth rate almost coincides with the uncorrugated case. As $k_{w}$ becomes larger, the overlap portion becomes larger, also for long waves, which is different from panel (a), the relatively small $k_{w}$ case. Furthermore, part of the unstable modes overlap the uncorrugated tube case completely as seen in the panel (c) and (d). When the wall shape has sufficiently short wave modes, a stable band is observed for $1<\epsilon k<\epsilon k_{w}-1$ (see panel (c) and (d)), which appears periodically as well. One cutoff mode $\epsilon k_{c}=1$ is from the uncorrugated tube branch while the other is from the first wall harmonic $|\epsilon k_{c}-\epsilon k_{w}|$. The result is similar to the findings in [75], which is as expected, since, although the lubrication model is only valid for asymptotically thin annular film, it in principle could reflect the basic features about the interaction between the capillarity and the wall harmonics.

The excitation mechanism persists even for sufficiently small $\sigma$ in our long wave model. Considering the small $\sigma$ limit, namely, $\sigma\rightarrow 0$ and setting $g=(4\lambda_{0}G)^{-1}$ with $\lambda_{0}=1$, an expansion of $g$ in $\sigma$ gives

Then the equation (37) becomes

The first two terms of the above equation form the linearized equation for a straight tube ($\sigma=0$) for the dispersion $\omega$. Rewriting the above equation leads to

Then it is clear that even very small corrugation amplitude would bring in a small correction term to the growth rate for straight tube case. It is just that the time requires to see different behavior (due to corrugation) from the straight tube case is proportional to $\sigma^{-1}$ roughly (if the coefficient of $\sigma$ happens to be zero, then the correction comes from higher order terms in $\sigma$), namely, longer time is needed for small $\sigma$ to show the effect of corrugation.

The effect of wall wavenumber $k_{w}$ on the maximum growth rate $\omega_{m}$ is illustrated in the left panel of figure 3 which shows that $\omega_{m}$ reaches a local minimum around $k_{w}\approx 5$ for $d_{0}=2$ that is close to the most unstable mode $k_{m}(\sigma=0,d_{0}=2)=k^{0}_{m}\approx 4.5$ in the straight tube case. This minimum value is also close to the uncorrugated case $\omega_{m}(\sigma=0)$ when $\sigma$ is small. Meanwhile, the growth rate also tends to have a local maximum when $k_{w}\approx 0\ {\rm or}\ 9(\approx 2k^{0}_{m})$. For a larger $\sigma$, this local minimum (maximum) has a smaller (larger) value. The parameter $d_{0}=2$ is chosen in figure 3 but other choices of $d_{0}$ showed qualitatively similar results. For relatively large $k_{w}$, $\omega_{m}$ varies little and $\omega_{m}$ is seen to reach a smaller value for a larger $\sigma$. As will be shown in the nonlinear simulations later, the relatively large $k_{w}$ leads to a tube with an effective radius $d_{0}-\sigma$. Therefore a larger $\sigma$ (with large $k_{w}$) means an effectively narrower tube which enhances the confinement and gives smaller growth rate.

Meanwhile, we investigate the variation of the associated wavenumber $k_{m}$ at which $\omega_{m}$ is obtained in the right panel of figure 3. For most cases, there are two locations where $\omega_{m}$ is reached within one $k_{w}$ period. We only highlight the one within the long wave range or the smaller one. The other one can be calculated easily by using $k_{w}-k_{m}$, since $|k_{w}-k_{m}-k_{w}|=k_{m}$ is the first wall harmonic of the perturbation. For relatively small $k_{w}$, $k_{m}$ is found to be shifting between a small value close to zero and $k_{w}/2$ approximately. Corresponding to the local minimum from the left panel of figure 3, $k_{m}=0.5$ is found in the right panel (for $\sigma=0.1,0.2$). Then the other mode $|k_{m}-k_{w}|=4.5$ is found to approach $k_{m}^{0}$. For relatively large $\sigma$ ($\sigma=0.4,0.6$ in figure 3), $k_{m}\approx 0$ so the other one is close to $5$ which indicates an obvious shift in the dominant modes. As $k_{w}>5$ increases, $k_{m}$ also increases until $k_{w}$ is sufficiently large, where $k_{m}$ oscillates around 5, which corresponds to the case when the unstable branch from the uncorrugated case has been recovered (see the similar case in the figure 2 for $d_{0}=5$). But there still exists small deviation in $\omega_{m}$ and $k_{m}$ from the exact value in the uncorrugated straight tube case. Similar deviation phenomenon in $\omega_{m}$ and $k_{m}$ was also pointed out in [75].

In figure 4, we compare the growth rates calculated directly from the long wave model (21), (22) with the results based on the lubrication model (28). In order to make the comparison, $\epsilon=1$ is taken and a thin annular layer is considered with $d_{0}=1.1$. Based on (25), $\delta=0.1$ and $\sigma=0.02\ (\sigma^{\prime}=0.2)$ are used in the full long wave model. Suppose $\tilde{\omega}$ is the growth rate obtained from (28), the true growth rate in the long wave model is calculated as $\delta^{3}\tilde{\omega}$ which is plotted in figure 4. Two typical wall wavenumbers are chosen in figure 4, where the agreement is good, although the lubrication model seems to overestimate the growth rate slightly for this value of $d_{0}(\delta)$ in the left panel, but in general has predicted the modes associated with the maximum growth rate. In the next section, we carry out direct numerical simulations to further validate the results we obtained for the linear stability.

In the results reported here, we consider the wall shape as $d=d_{0}+\sigma\cos(k_{w}z)$, with initial and (no axial flux) boundary conditions given by

To ensure both the wall and interface are periodic in the computational domain, the period length is taken as $L=\pi/\beta$ where $\beta$ is the common divisor such that $(k,k_{w})=(n_{k},n_{w})\beta$ with both $n_{k,w}$ integers (see also [76]). In addition, $\epsilon=0.1$, $\lambda_{0}=1$ and $A_{0}=0.05$ are fixed in simulation unless otherwise stated. Since $A_{0}=0.05$ is small, the early stage evolution is expected to follow the linear theory which will be confirmed numerically later. The nonlinear equations are solved by the solver EPDCOL [40] which has been successfully used for related problems ([16] e.g.), and our previous studies ([73, 74]). This solver uses the finite element discretization in space and advances in time through the Gear’s method. In the results reported here, about $1600\sim 3200$ space points are used within the computational domain. The simulation is stopped when the neck of the core thread $S_{min}<0.003$ or the vertical gap between the wall and interface is smaller than $0.0008$, which indicates a touching solution or annular film rupture rather than pinching of the core liquid thread in our problem.

Our numerical results by choosing an initially slightly perturbed thread interface show fair agreement with the linear theory for which we solved by using the FFH method. A few typical results are shown in figure 5. In the panel (a), the evolution is shown for $k=7.5,k_{w}=3,d_{0}=5$, in which case the growth rate is larger than the straight tube case based on the inset, where the star symbol indicates the parameter values used in simulation. In fact, due to corrugation, the disturbance is modulated by exp$(i(k\pm nk_{w}))$ ($n=1,2,\cdots$) and the eventual dominant growth rate is associated with the mode (first wall harmonic) $|k-k_{w}|=4.5$ (or $4.5\pm nk_{w}$ due to the periodicity). It is shown that the interface evolves first as if there were no wall roughness and the growth rate follows the value in the straight tube case, for $k=7.5,\sigma=0,d_{0}=5$ at early stage. After $t\approx 25$, the interface recognizes the wall topography effect and then follows the new dominant growth mode that is predicted by our FFH method. The agreement is shown by viewing the solid line from the direct simulation and the dashed line from linear theory agree with each other well, even though in this particular case, the time period showing the agreement is relatively short and the amplitude of the perturbation grows to finite for which the linear theory is not strictly valid. Our results do not indicate that the linear theory is able to predict finite amplitude growth. The agreement merely coincides with the fact that the linear theory usually works well even beyond its valid domain. More details of the interface evolution will be revealed in figure 6. The important message that our numerical results delivered is the existence of the transition (see also the simulations in [75]). It is as expected and can be explained using our (33), which indicates that the growth rate $\omega(\sigma=0)$ from the uncorrugated case dominates in short times and there exists a critical time when the correction term due to corrugation eventually becomes the same order as $\omega(\sigma=0)$. The transition exists for most of cases as long as the growth rate differs from $\omega(\sigma=0)$. But it will not appear if the thread experiences breakup before reaching the critical time.

In panel (b) and (c), we fix $k_{w}=13,d_{0}=5$ but choose two different wavenumbers for the initial condition: one corresponds to a mode that is roughly equal to the uncorrugated tube case, $k=4$ while the other one corresponds to a mode that would have been stable if the wall is uncorrugated, $k=11$. Then it is not surprising that the interface follows a single growth rate from the linear theory in panel (b). In panel (c), we again observe some transient growth first. At early stage, the relatively shortwave perturbation saturates (the growth rate is negative for $k=11$ in the straight tube case) and after $t\approx 43.4$, the evolution follows the growth rate from the linear theory for $\sigma\neq 0$. Panel (d) illustrated a similar case as panel (c) but with a smaller $d_{0}$, i.e. the wall is initially closer to the thread interface, which shows qualitatively the same behavior as panel (c). Such transient growth from saturated shortwave perturbations to unstable long wave ones is also observed in [75], in which the annular layer is asymptotically thin and the core dynamics is neglected.

We also carry out a close inspection in the case in panel (a) and (c) of figure 5 by showing the interface profiles in figure 6 and 7 respectively. Corresponding to the case in panel (a) of figure 5, the interface evolves first as the case with no corrugation before $t\approx 25$ which is shown in the panel (a) of figure 6, while the panel (b) shows the evolution after the transition, which, despite the fact that some short waves surf on the large wave crest, follows longer wave length perturbation compared to the initial relatively short wave length perturbation. In particular, the eventual dominant mode is $|k-k_{w}|=4.5$ instead of $7.5$ as discussed above. The wall shape is illustrated in the upper panels only in order to view the relative phase between the interface shape and the wall. For the case in panel (c) of figure 5, the transition from stabilizing to destabilizing is obvious. As indicated in figure, we plotted the thread profiles before and after the turning point $t\approx 43.4$. In the panel (a) of figure 7, the perturbation is indeed saturated as the amplitude of the perturbation decreases, while in the panel (b), the perturbation is amplified obviously and the core tends to pinch off eventually. This result is in good qualitative agreement with those in [75]. The corresponding evolution of $S_{min}$ in bottom panel (c) further confirmed this transient evolution, where the circle indicates $S_{min}$ reached a maximum as time advances.

The dashed lines in the panels (c) of figure 6 and 7 are the results when inertia terms are included in calculations. Similar to the findings in [74], the inertia only affects the transient evolution slightly while having little effect both in the early stage and near pinching. It is also consistent with the scaling analysis in [51] where inertia is shown to give way to viscous fluid-fluid interaction when the region is slender. Therefore, we focus mainly on the Stokes flow problem, i.e. $\bar{R}_{e}=0$, for the discussion here.

In this subsection, we investigate the drop formation toward the breakup of the liquid threads or layers. In figure 8, drop and satellite formations for various wall shapes are shown when the tube wall is initially far away from the thread interface, $d_{0}=5$, so that the tube wall is out of sight in all the panels (the same tube wall shapes will be shown in figure 9). In the top panels of figure 8, the tube wall is kept uncorrugated, $\sigma=0$, while $k_{w}=3,13,23$ with $\sigma=0.2$ for the following three ones from the top to the bottom. It is seen that the drop formation can be quite different when the wall has wavy structure. In the panels in the second row, $(k,k_{w})=(7,3)$ corresponds to a point that leads to a larger linear growth rate than the straight tube case (see panel (a) in figure 2), which causes the thread to pinch with fewer large or mother drops than the top panel. The satellite formations are illustrated in the right column of figure 8. In addition, further short-wave perturbation causes the surface of large drop to have dimples. When $k_{w}=13$, the panel (b) of figure 2 shows the linear growth rate dominated by the resonant one $|k-k_{w}|=6$ which is close to $k=7$, while $k_{w}=23$, one unstable branch completely covers the growth rate curve in the case of no wall corrugation. Therefore, it is expected that the eventual drop shapes in the third and bottom row are similar to the top panel. However, a close inspection reveals that, although the drop shape is similar, the detail of pinch-off differs. For $\sigma=0$, the satellites are obtained simultaneously since the actual period is $2\pi/7\approx 0.9$ in the top row of figure 8. For $\sigma=0.2$ and $k_{w}=13$ in the third row, the middle satellite is obtained via the first breakup while the rest possible satellites can only appear from the subsequent breaks up (since the necking region has already formed, it is likely they will be the satellite drops; see the related work in [70] which presented a ’self-repeating’ mechanism on a series of breakups to obtain satellites). In the bottom row, the left (right) most satellite together with the middle satellite breakups earlier than the ones next to them. The calculation after the first pinch off may be interesting, but this is beyond the scope of this paper. The difference in droplet sizes (together with the cases having different values of $d_{0}$) will be summarized later in figure 12.

Since the tube wall would appear far away from the local pinching region (assuming the tube cross section variation is not too large), the pinching behavior would be expected to resemble the case with no wall corrugation (see the study on the pinching dynamics in [74]). Therefore this is not pursued further in current work.

When the annular layer is relatively thin, namely, the averaged tube radius $d_{0}$ is small, we show drop formation for $d_{0}=1.5,1.3$ in figure 9, where the wall shapes are chosen to be the same as in figure 8 but they are visible in figure now. In the top two panels, the tube is straight, $\sigma=0$ and we obtain the so called plug-with-collar drop formation (see [32] and [74]) where the fluids are trapped between the tube wall and the large drops along with pinch off. When the tube wall has variation in shape, the thread needs to cope with the wall constraint and larger drops can be expected at first breakup even starting with the same initial perturbation wavenumber. The drop shapes resemble the results in [55] in the full axisymmetric simulations when the wall variation is present. The panels in the second row of figure 9 show larger drops in the middle of the thread compared to the top panels where the tube wall is straight, because a longer wave length perturbation is excited by the wall harmonics, $|k-k_{w}|=4<k=7$ (see also the argument in [75]). In the third row, pinching occurs near $z=0$, leaving a long thread both at $z>0$ and $z<0$. Additionally, the left one in the third row would perhaps experience a second breakup due to the thin neck around $z\approx\pm 1$, which leads to multiple drop formations. Meanwhile, in some portion of the tube, the annular layer tends to breakup as well, since as pinching occurs, the film drainage regime is also approached slowly. Due to the wall topography, it is seen that more fluids compared to the case $\sigma=0$ can be trapped between the wall and drops. In the last row of figure 9, drop formation is seen to differ significantly. In the left panel, a ultimate multiple drop formation is expected, while in the right panel, a large size drop is obtained in the middle which is connected by drops in comparable sizes. As shown later, this case in the right bottom panel stands near a transition from pinching solutions to film rupture solutions.

The associated evolutions of the corresponding minimum thread radius and the minimum thread-tube gap are plotted in figure 10, where the left panel for the case $d_{0}=1.5$, shows as $S_{min}\rightarrow 0$, i.e. the core thread pinches, the minimum gap still remain significantly larger compared to $S_{min}$ or at least, $S_{min}$ decays faster than the gap value. Therefore pinching solutions are obtained (at least up to $k_{w}=43$ here). As shown in the thin annulus limit, owing to the different time scales (recall $\tau=\delta^{3}t$ with $\delta\ll 1$ in (27) and (28)), pinching is expected to occur in advance of film drainage that happens in longer time dynamics, and perhaps even after the first pinch off, as shown in [32] in the straight uncorrugated tube case. In addition, we observe from the left panel of figure 10 that the change in minimum gap along pinching with $k_{w}$ is not monotonic. It seems that the gap remains larger for $k_{w}=43$ than $k_{w}=33$. The right panel of figure 10 shows similar results to the left panel, except that $d_{0}=1.3$ and for sufficiently large $k_{w}$ ($k_{w}>23$), the neck of core thread $S_{min}$ remains order one while the gaps seem to reach a very small value first. In those cases, we provided the numerical evidence that a tube wall with large wavenumber may induce the film drainage regime, meanwhile suppress pinching even for a case where the annulus is not asymptotically thin. In the case without wall corrugation, $\sigma=0$, to suppress pinching, sufficiently thin annulus is required (the numerical work by [54] gave a threshold value $d_{0}\approx 1.19$ for the uncorrugated tube in Stokes flow); in contrast, we found $d_{0}=1.3$ (namely, the average film thickness $0.3$) is sufficiently small to produce this phenomenon. In general, this should also depend on the size of wall variation $\sigma$ (as well as the initial conditions on the thread profiles). For the discussion here, $\sigma=0.2$ is fixed all the time and we merely provide the evidence of possible pinching suppression for an annular layer whose thickness is not asymptotically small. Figure 11 shows the thread interface shape at the final stage of our simulations, at $t=257.1$ and $t=417.3$ respectively, where it seems that the effective tube radius is reduced ($\sim d_{0}-\sigma$) due to the rapid variation of tube cross section, which makes the fluids inside the wall undulation effectively rigid. Notice that in figure 11, the core fluid does not intrude the inside of the wall undulation part, while in figure 9 (the third and fourth row in particular), the core fluid still can wrap around the undulation ’tips’. As the liquid layer keeps thinning, locally, one still expects the lubrication model holds as well as the scaling laws based on it (see [33] and [50], where $h\sim t^{-1}$ for a collar next to a collar, $h\sim t^{-1/2}$ for a lobe next to a collar, with $h$ the minimum film thickness to the substrate and $t$ the time). Unfortunately our numerical results fail to show the scalings even the gap is small and in the order of $O(10^{-4})$. It is possible that the thin film regime has not been fully reached yet as the fluid interface is only touching the wall undulation tip and the fluid region adjacent to the thinning part is different from the collar or lobe that is seen in [50].

Finally, for pinching solutions, we show the effect of the wall wavenumber $k_{w}$ on the breakup times and drop sizes, with the latter characterized by an effective drop radius $R_{eff}$ that is defined as

where $z_{1},z_{2}$ are two points associated with sufficiently thin thread neck that we consider to be the pinching points. In the left panel of figure 12, the breakup time $t_{b}$ first decreases and then increases as $k_{w}$ varies between $0$ and $10$. This is expected from our linear theory (also qualitatively similar to the results in figure 3). When $k-k_{w}$ corresponds to a mode with larger growth rate, the breakup time is smaller. For example, when $d_{0}=1.5$, the dominant mode is $k=4$ for $k_{w}=3$ and $k=7$ or $1$ for $k_{w}=8$. Accordingly, the growth rate is $0.06$ and $0.049$ respectively, which indicates the former case breakups earlier than the latter one. As $k_{w}$ increases, $t_{b}$ increases and the increasing is more profound when $d_{0}$ is smaller. For $d_{0}=2,5$, the breakup time eventually approaches the value in the straight tube case for large $k_{w}$ approximately, because the initial wavenumber is already a long wave one, and one of the unstable branches in the linear theory for $\sigma>0$ has almost covered the straight tube case (see figure 4). While for the relatively small $d_{0}$, the annular film also tends to reach the drainage regime which occurs in a longer time scale and tends to suppress the pinching. As given in figure 11, the annular film touching solution is obtained instead of pinching for $k_{w}=33$ and the interface profile that is shown is at $t=257.1$. The right panel of figure 12 illustrates the distribution of a large and satellite drops. For clarity here, we only calculate the satellite size from the first breakup (typically the one satellite in the middle) together with the adjacent large drops, where we estimated a second pinching point, as our current code can only track down to the first breakup. But for $k_{w}=3,d_{0}=5$ in figure 8, there are clearly three large drops and two small satellites, which are all taken into account. This has also been done for the small $d_{0}$ cases. For $d_{0}=5$ in figure 8, multiple satellite formation is expected eventually, which is similar for $d_{0}=2$ (data not shown). As $k_{w}$ increases, it is seen that the size of both the large and satellite drops increases to a local maximum for some $k_{w}$ within the range $0<k_{w}<10$, or $0<\epsilon k_{w}<1$, which depends on the value of $d_{0}$. The drop sizes are similar to the straight tube case when $k_{w}$ is relatively large except when $d_{0}$ is sufficiently small. For $d_{0}=1.3$, the size of the satellite drop keeps decreasing as $k_{w}$ increases. Based on our numerical results in figure 10, this indicates a transition from the pinching to touching solution.