2.1. Scalings and governing equations in dimensionless form

We now rescale all relevant variables. Specifically, we rescale lengths by $a$ , velocities by $a\varOmega$ , stresses by $(\mu _s + \mu _p)\varOmega$ , and time by $1/\varOmega$ :

(2.11) \begin{align} &Z = \frac {z}{a}, \quad R = \frac {r}{a},\quad H = \frac {h}{a},\quad \boldsymbol {U} = \frac {\boldsymbol {u}}{a\varOmega }, \quad T = t\varOmega ,\notag \\ &\boldsymbol{\boldsymbol{E}} = \frac {\unicode{x1D640}}{\varOmega },\quad P =\frac {p}{(\mu _s+\mu _p)\varOmega }, \quad \boldsymbol{\varSigma }' = \frac {\boldsymbol{\sigma }}{(\mu _s+\mu _p)\varOmega }, \quad \boldsymbol{\varSigma }_p = \frac {\boldsymbol{\sigma }_p}{(\mu _s+\mu _p)\varOmega }. \end{align}

These dimensionless variables are accompanied by dimensionless parameters for viscosity ratios

(2.12) \begin{equation} \beta _p = \frac {\mu _p}{\mu _s + \mu _p} \quad \text{and} \quad \beta _s = 1 - \beta _p, \end{equation}

the Weissenberg number $\textit{Wi}$ , the gravity parameter ${\mathcal{G}}$ , and the capillary number $\textit{Ca}$ :

(2.13) \begin{equation} Wi = \lambda \varOmega , \quad \mathcal{G} = \frac {\rho g a}{(\mu _s + \mu _p)\varOmega }, \quad Ca = \frac {(\mu _s+\mu _p)a\varOmega }{\gamma }. \end{equation}

The Weissenberg number $\textit{Wi}$ , defined as the product of the fluid’s relaxation time and the characteristic shear rate, is commonly used to assess the viscoelastic properties of a material, indicating whether its behaviour is more elastic or more viscous. In experiments, $\textit{Wi}$ typically ranges from 1 to 10 or higher. In this work, however, we restrict attention to sufficiently small $\textit{Wi}$ , allowing us to employ a regular perturbation expansion in $\textit{Wi}$ . Although it is difficult to measure the interface height experimentally when ${\textit{Wi}} \lt 1$ , we anticipate that our analysis may provide insights applicable to larger $\textit{Wi}$ . Another key parameter is the Deborah number $De$ , which represents the ratio of the material’s response time to the characteristic time of the flow. In this problem, $\textit{Wi}$ and $De$ are equivalent, so no distinction is made between the two. Here, we focus on the case where ${\textit{Wi}} \ll 1$ , corresponding to weak viscoelastic effects. The ranges of ${\mathcal{G}}$ and $\textit{Ca}$ will depend on the value of $\textit{Wi}$ , so that we can derive analytical expressions for the small deformation of the interface shape $H$ , as will be discussed in §3.1.

Using the scalings in (2.11), the governing equations (2.1), (2.2) and (2.4) become

(2.14a) \begin{equation}{ \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {U} = 0, \quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\varSigma }' = \mathcal{G} \boldsymbol {e}_z,} \end{equation}

(2.14b) \begin{equation}{ \boldsymbol{\varSigma }' = -P \boldsymbol {I} + 2\beta _s \boldsymbol{E} + \boldsymbol{\varSigma }_p,} \end{equation}

(2.14c) \begin{equation}{ \boldsymbol{\varSigma }_p + Wi\left (\frac {\partial \boldsymbol{\varSigma }_p}{\partial T} + \boldsymbol {U} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\varSigma }_p - \boldsymbol{\varSigma }_p \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol {U} - (\boldsymbol{\nabla }\boldsymbol {U})^{\textrm {T}}\boldsymbol{\cdot }\boldsymbol{\varSigma }_p\right ) + \frac {\alpha\, Wi}{\beta _p}\boldsymbol{\varSigma }_p \boldsymbol{\cdot }\boldsymbol{\varSigma }_p = 2\beta _p\boldsymbol{E} . } \end{equation}

Similarly, the corresponding boundary and initial conditions (2.5)–(2.10) become

(2.15a) \begin{equation}{ \text{No-slip boundary condition:}\quad \boldsymbol {U} = \boldsymbol {e}_{\theta } \quad \text{at} \quad R =1, } \end{equation}

(2.15b) \begin{equation}{ \text{Contact angle boundary condition:}\quad \frac {\partial H}{\partial R} = 0 \quad \text{at} \quad R =1, } \end{equation}

(2.15c) \begin{equation}{ \text{Stress boundary condition:} {\quad \boldsymbol {n} \boldsymbol{\cdot }\boldsymbol{\varSigma }' = -\frac {1}{\textit{Ca}}(\boldsymbol{\nabla} _s\boldsymbol{\cdot }\boldsymbol {n})\boldsymbol {n}} \quad \text{at} \quad Z = H(R,T), } \end{equation}

(2.15d) \begin{equation}{ \text{Kinematic boundary condition:} \quad \frac {\partial H}{\partial T} = U_z - U_r\frac {\partial H}{\partial R} \quad \text{at} \quad Z = H(R,T), } \end{equation}

(2.15e) \begin{equation}{ \text{Far field:} \quad H \rightarrow 0, \quad \boldsymbol {U} \rightarrow \boldsymbol {0}, \quad { \boldsymbol{\varSigma }_p \rightarrow \boldsymbol {0}} \quad \text{as} \quad R \rightarrow \infty,} \end{equation}

(2.15f) \begin{equation}{ \text{Initial condition:} \quad H = 0, \quad { \boldsymbol{\varSigma }_p = \boldsymbol {0}} \quad \text{at} \quad T = 0. } \end{equation}

To further simplify the momentum equation (2.14a ), it is useful to introduce a modified pressure variable above hydrostatic pressure and redefine the stress tensor to incorporate the conservative gravitational force:

(2.16) \begin{equation} \boldsymbol{{\varSigma }} = {\boldsymbol{\varSigma }}' - {\mathcal{G}}Z{\boldsymbol {I}}, \quad {\mathcal{P}} = P + {\mathcal{G}}Z. \end{equation}

With these reformulated variables, the total stress equation in the fluid (2.14b ) and the momentum equation (2.14a ) simplify to

(2.17) \begin{equation} \boldsymbol{{\varSigma }} = -\mathcal{P} \boldsymbol {I} + 2\beta _s \boldsymbol{E} + \boldsymbol{\varSigma }_p, \quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{{\varSigma }} = \boldsymbol {0}. \end{equation}

To account for the new variables, the stress boundary condition (2.15c ) becomes

(2.18) \begin{equation} \boldsymbol {n} \boldsymbol{\cdot }\boldsymbol{{\varSigma }} = -\left (\mathcal{G}Z +\frac {1}{\textit{Ca}}(\boldsymbol{\nabla} _s\boldsymbol{\cdot }\boldsymbol {n})\right ) \boldsymbol {n} \quad \text{at} \quad Z = H(R,T). \end{equation}

3.1. Distinguished limit

We focus on three key non-dimensional parameters, i.e. the Weissenberg number $\textit{Wi}$ , the gravity parameter ${\mathcal{G}}$ , and the capillary number $\textit{Ca}$ ; we assume that the viscosity ratio $\beta _p$ is given. It is natural to examine the asymptotic limit where ${\textit{Wi}} \ll 1$ in which viscoelasticity appears as a correction to the Newtonian flow. However, careful consideration of the ranges of ${\mathcal{G}}$ and $\textit{Ca}$ is required, as different values of these parameters will lead to different asymptotic regimes depending on $\textit{Wi}$ . Before advancing our analysis, we will investigate how key quantities, e.g. velocity, pressure, polymeric stress and interface height, depend on these non-dimensional parameters, to gain insight into the distinguished limit that forms the foundation of our study.

First, we consider the relationship between the gravity parameter and the capillary number. A natural dimensionless quantity that characterises the ratio of gravitational to capillary forces is the Bond number, defined as

(3.1) \begin{equation} B = \frac {\rho g a^2}{\gamma } = \textit{Ca}\,\mathcal{G}. \end{equation}

In the following sections, we aim to incorporate both effects into our analysis. To this end, we assume that the Bond number is of order unity, i.e. $B = O(1)$ . This allows us to write the capillary number in terms of the Bond number and the gravity parameter, $Ca = B/\mathcal{G}$ .

Next, we consider the relationship between the Weissenberg number and the gravity parameter. Initially, we assume that the polymers are in equilibrium. As the rod begins rotating at a constant angular velocity, the fluid forms a vortex around the rod, causing the suspended polymers to stretch. Once the polymers extend beyond their equilibrium state, they induce additional polymeric stress of order $O(Wi)$ :

(3.2) \begin{equation} \boldsymbol{{\varSigma }}_{\textit{p}} \thicksim \textit{Wi}. \end{equation}

When viscoelasticity is absent from the flow, the Stokes flow only yields a flat interface, i.e. in the absence of inertial effects. Therefore, we expect the normal stress near the interface to arise from the polymeric stress. Using the normal stress boundary condition (2.18), we find that the normal stress near the interface (deformation $H$ ), in the limit of small deformation, behaves as

(3.3) \begin{equation} \boldsymbol{{\varSigma }} \thicksim \boldsymbol{{\varSigma }}_{\textit{p}} \thicksim \mathcal{G}H \quad (\text{near the interface}). \end{equation}

Here, we have assumed that the gravitational and capillary effects are comparable when $B = O(1)$ . By combining (3.2) and (3.3), we anticipate that the typical interface height will depend on both viscoelasticity and gravity as

(3.4) \begin{equation} H \thicksim \frac {\textit{Wi}}{\mathcal{G}}. \end{equation}

As discussed earlier, the dynamics of the interface is driven by the (weak) viscoelasticity of the fluid. Therefore, we expect the characteristic time scale for temporal changes in the interface height to be determined by the polymer relaxation time, rather than the flow time scale itself. With this in mind, we use the kinematic boundary condition (2.15d ) and the typical interface height (3.4) to estimate the typical magnitude of the upward fluid velocity near the interface as

(3.5) \begin{equation} U_{\textit{interface}} \thicksim \frac {H}{\textit{Wi}} \thicksim \frac {1}{\mathcal{G}}. \end{equation}

Near the interface, we observe from (3.2) that the polymeric stress scales as $\textit{Wi}$ , whereas from (3.5), the viscous stress scales like $\mathcal{G}^{-1}$ . When $1 \gg \textit{Wi} \gg \mathcal{G}^{-1}$ , the normal stress at the interface comes mainly from the polymeric stress, as we had assumed to obtain (3.3).

However, when $\mathcal{G} = {\textit{Wi}}^{-1}$ , the upward flow near the interface induces a viscous stress comparable to the polymeric stress, adding further complexity to the problem. The distinguished limit $\mathcal{G} = {\textit{Wi}}^{-1}$ , where our qualitative analysis (3.2)–(3.5) begins to break down, motivates us to explore the gravity parameter ${\mathcal{G}}$ within the range $\mathcal{G} \gtrsim {\textit{Wi}}^{-1}$ . To simplify the complexity of the double asymptotic expansion into a single expansion in terms of $\textit{Wi}$ , based on the distinguished limit $\mathcal{G} = {\textit{Wi}}^{-1}$ , we introduce a rescaled gravity parameter ${G}$ by defining

(3.6) \begin{equation} {G} = {\textit{Wi}}\mathcal{G}, \end{equation}

where ${G}$ is now taken to be of order $O(1)$ . Intuitively, when gravitational effects are large ( ${G} \gg 1$ ), the viscoelastic fluid has difficulty climbing the rotating rod, which leads to only small deformations on the interface shape. This is precisely the small-deformation situation that we aim to address in this study.

3.2. Perturbation expansion in Weissenberg number

In the following sections, we expand the velocity, pressure and the polymeric stress as regular perturbations in ${\textit{Wi}} \ll 1$ ,

(3.7a) \begin{equation}{{P} = \mathcal{P}^{(0)} + {\textit{Wi}} \mathcal{P}^{(1)} + \cdots , } \end{equation}

(3.7b) \begin{equation}{ \boldsymbol{E} = \boldsymbol{E}^{(0)} + {\textit{Wi}} \boldsymbol{E}^{(1)} + \cdots , } \end{equation}

(3.7c) \begin{equation}{\boldsymbol{\varSigma }_p = \boldsymbol{\varSigma }_p^{(0)}+ {\textit{Wi}} \boldsymbol{\varSigma }_p^{(1)} + \cdots , } \end{equation}

where the superscript indicates the corresponding order in $\textit{Wi}$ . Along with the quantities expressed as expansions in (3.7), we introduce an expansion of the interface shape in the situation where the interface is slightly deformed:

(3.8) \begin{equation} H(R,T) = {\textit{Wi}}^2 H^{(2)} + {\textit{Wi}}^3 H^{(3)} + \cdots. \end{equation}

Although it is more natural to start the expansion of $H$ at order $O(\textit{Wi})$ , since the interface remains flat in the absence of viscoelasticity, we have opted to express it as in (3.8) for two reasons. First, we observe from our preliminary scaling argument (3.4) that $H \sim {\textit{Wi}}\mathcal{G}^{-1} \sim {\textit{Wi}}^2\, {G}^{-1} = O (\textit{Wi}^2 )$ . Second, as we will see in §3.3, expanding the free-surface boundary condition to be consistent with the perturbation expansion in the Weissenberg number involves extensive calculations unless each term is scaled correctly. By keeping the non-zero leading-order expansion of $H$ , the derivations for the boundary conditions at each order become more organised. In this work, we focus on obtaining an analytical formula for the leading-order term of the interface shape $H^{(2)}$ , given in (3.8).

3.3. Domain perturbation method for the interface boundary condition

In §2.1, we introduced the normal stress boundary condition (2.18) and the kinematic boundary condition (2.15d ), in which both boundary conditions are evaluated at the interface $Z = H$ . However, when expanding the interface shape $H$ in terms of $\textit{Wi}$ (see (3.8)), it is challenging to evaluate the boundary conditions using the full series expansion. To resolve this issue, we take advantage of the assumption that the interface deformation is small to simplify the boundary conditions using a domain perturbation expansion. This step enables us to systematically write down the corresponding interface boundary conditions for the velocity, pressure and polymeric stress fields at each order in $\textit{Wi}$ .

We start by expressing the normal stress boundary condition (2.18) at each order in $\textit{Wi}$ . The assumption of small deformation of the interface provides two simplifications. First, we can approximate the normal vector $\boldsymbol {n}$ at the interface as $\boldsymbol {n} \approx \boldsymbol {e}_z - ({\partial H}/{\partial R})\boldsymbol {e}_r = \boldsymbol {e}_z + O(\textit{Wi}^2)\boldsymbol {e}_r$ . Second, when evaluating the boundary condition at $Z = H$ , the small value of $H$ allows us to perform a Taylor expansion around $Z = 0$ . From these two remarks, we may simplify the left-hand side of (2.18) as

(3.9) \begin{align} \boldsymbol {n} \boldsymbol{\cdot }\boldsymbol{{\varSigma }}\bigg|_{Z = H} \bigg.\vphantom{\frac {\partial (\boldsymbol {n} \boldsymbol{\cdot }\boldsymbol{{\varSigma }})}{\partial Z}}&= \boldsymbol {n} \boldsymbol{\cdot }\boldsymbol{{\varSigma }}\bigg|_{Z = 0} \bigg. + H \frac {\partial (\boldsymbol {n} \boldsymbol{\cdot }\boldsymbol{{\varSigma }})}{\partial Z} \bigg|_{Z = 0} \bigg. + \left .\frac {1}{2}H^2 \frac {\partial ^2 (\boldsymbol {n} \boldsymbol{\cdot }\boldsymbol{{\varSigma }})}{\partial Z^2} \right|_{Z = 0} + \cdots \nonumber \\ &= \boldsymbol {e}_z \boldsymbol{\cdot }\left . \vphantom{\frac {\partial (\boldsymbol {n} \boldsymbol{\cdot }\boldsymbol{{\varSigma }})}{\partial Z}}\left (-\mathcal{P}^{(0)}\boldsymbol {I} + 2\beta _s \boldsymbol{E}^{(0)} + {\boldsymbol{\varSigma }_p^{(0)}}\right )\right|_{Z = 0} \notag \\ &\quad {}+ \textit{Wi} \left (\boldsymbol {e}_z \boldsymbol{\cdot }\left (-\mathcal{P}^{(1)}\boldsymbol {I} + 2\beta _s \boldsymbol{E}^{(1)} + {\boldsymbol{\varSigma }_p^{(1)}}\right )\bigg|_{Z = 0}\right ) + O(\textit{Wi}^2). \end{align}

On the other hand, the small-deformation assumption allows us to simplify the curvature of the interface ( $\boldsymbol{\nabla} _s \boldsymbol{\cdot }\boldsymbol {n}$ ) so that the right-hand side of (2.18) becomes

(3.10) \begin{align} -\left (\mathcal{G}Z + \left . \vphantom{\frac {\partial (\boldsymbol {n} \boldsymbol{\cdot }\boldsymbol{{\varSigma }})}{\partial Z}}\frac {1}{\textit{Ca}}(\boldsymbol{\nabla} _s\boldsymbol{\cdot }\boldsymbol {n})\right ) \boldsymbol {n}\right|_{Z = H} &= -\left ( \mathcal{G}H-\frac {1}{\textit{Ca}}\frac {1}{R}\,\frac {\partial }{\partial R}\left (R\frac {\partial H}{\partial R}\right )\right )\boldsymbol {e}_z + O(\textit{Wi}^2) \notag \\ &= -\textit{Wi}\, {G}\left (H^{(2)}-\frac {1}{\textit{BR}}\,\frac {\partial }{\partial R}\left (R\frac {\partial H^{(2)}}{\partial R}\right )\right )\boldsymbol {e}_z + O(\textit{Wi}^2). \end{align}

Combining (3.9) and (3.10), we conclude that the stress boundary condition (2.18) at $O(1)$ and $O(\textit{Wi})$ becomes

(3.11a) \begin{equation}{ \kern-.8pc \boldsymbol {e}_z \boldsymbol{\cdot }\left . \vphantom{\frac {\partial (\boldsymbol {n} \boldsymbol{\cdot }\boldsymbol{{\varSigma }})}{\partial Z}}\left (-\mathcal{P}^{(0)}\boldsymbol {I} + 2\beta _s \boldsymbol{E}^{(0)} + {\boldsymbol{\varSigma }_p^{(0)}}\right )\right |_{Z = 0} = \boldsymbol {0}, } \end{equation}

(3.11b) \begin{equation}{ \kern-.8pc \boldsymbol {e}_z \boldsymbol{\cdot }\left . \vphantom{\frac {\partial (\boldsymbol {n} \boldsymbol{\cdot }\boldsymbol{{\varSigma }})}{\partial Z}}\left (-\mathcal{P}^{(1)}\boldsymbol {I} + 2\beta _s \boldsymbol{E}^{(1)} + {\boldsymbol{\varSigma }_p^{(1)}}\right )\right |_{Z = 0} = {-{G}\left (H^{(2)}-\frac {1}{\textit{BR}}\,\frac {\partial }{\partial R}\left (R\frac {\partial H^{(2)}}{\partial R}\right )\right )\boldsymbol {e}_z.} } \end{equation}

At this stage, we choose not to express the kinematic boundary condition (2.15d ) for each order in $\textit{Wi}$ . As mentioned in §3.1, the temporal evolution of the interface height is expected to be governed by the polymer relaxation time, rather than the flow time scale. This implies that, starting from a flat interface, the interface will rise rapidly on a time scale of order $O(\textit{Wi})$ before approaching its steady state. This behaviour suggests a boundary-layer-like structure in time, prompting the examination of two time scales: $T = O(1)$ and $T = O(\textit{Wi})$ . In §4, we focus on $T = O(1)$ , representing the outer layer approximation, to recover the steady-state interface shape. In §5, we analyse the shorter time scale $T = O(\textit{Wi})$ using a stretched time variable, corresponding to the inner layer approximation, to study the transient interface shape. We will revisit the kinematic boundary condition (2.15d ) separately for each rescaled time in §§4 and5.

4.1. Leading-order solution

In this subsection, we examine the interface dynamics for $T = O(1)$ , representing the state when time has progressed beyond the initial start-up of the flow. To this point, we have provided only a general qualitative outline explaining why a distinction in time scales is necessary. Now we examine the evolution of the polymeric stress by writing (2.14c ) in non-dimensional form as

(4.1) \begin{equation} \textit{Wi}\left (\frac {\partial \boldsymbol{\varSigma }_p}{\partial T} + \boldsymbol {U} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\varSigma }_p - \boldsymbol{\varSigma }_p \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol {U} - (\boldsymbol{\nabla }\boldsymbol {U})^{\textrm {T}}\boldsymbol{\cdot }\boldsymbol{\varSigma }_p + \frac {\alpha }{\beta _p}\boldsymbol{\varSigma }_p \boldsymbol{\cdot }\boldsymbol{\varSigma }_p \right )= -\boldsymbol{\varSigma }_p + 2\beta _p\boldsymbol{E} . \end{equation}

Equation (4.1) features a first-order time derivative, with the highest-order term $ {\partial \boldsymbol{\varSigma }_p}/{\partial T}$ being multiplied by a small parameter $\textit{Wi}$ . This structure frequently arises in boundary layer problems, where singular perturbation analysis is required to understand the solution behaviour. Typically, the analysis is split into two regions, one where (4.1) is used as is, and another where time is rescaled to maintain the dominance of the highest-order term; in this subsection, we focus on the former case.

Before proceeding further, we note that using (3.7), the right-hand side of (4.1) enables us to express the leading-order polymeric stress in terms of the leading-order velocity field:

(4.2) \begin{equation} \boldsymbol{\varSigma }_p^{(0)} = 2\beta _p\boldsymbol{E}^{(0)}. \end{equation}

In addition, when $T = O(1)$ and $H = O (\textit{Wi}^2 )$ , we can derive the kinematic boundary conditions (2.15d ) up to first-order terms, by employing a domain perturbation technique similar to that used in §3.3, as

(4.3) \begin{equation} U_Z^{(0)}\Big |_{Z = 0} = U_Z^{(1)}\Big |_{Z = 0} = 0. \end{equation}

To proceed with the problem solution, we substitute the leading-order polymeric stress (4.2) into the leading-order momentum equation (2.17) to obtain

(4.4) \begin{equation} \boldsymbol {0} = \boldsymbol{\nabla }\boldsymbol{\cdot }\left (-\mathcal{P}^{(0)} \boldsymbol {I} + 2\beta _s \boldsymbol{E}^{(0)} + {\boldsymbol{\varSigma }_p^{(0)}}\right ) = -\boldsymbol{\nabla }\mathcal{P}^{(0)} + {\nabla} ^2 \boldsymbol {U}^{(0)}, \end{equation}

where $\beta _s+\beta _p=1$ , and from (4.2), we have $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\varSigma }_p^{(0)}= \beta _p\, {\nabla} ^2 \boldsymbol {U}^{(0)}$ . Further, substituting (4.2) into the leading-order stress boundary condition (3.11a ) yields

(4.5) \begin{equation} \boldsymbol {e}_z \boldsymbol{\cdot }\left .\left (-\mathcal{P}^{(0)}\boldsymbol {I} + 2 \boldsymbol{E}^{(0)} \right )\right |_{Z = 0} = \boldsymbol {0}. \end{equation}

The leading-order governing equation (4.4) is solved subject to the boundary conditions (2.15a ), (2.15e ), (4.3) and (4.5), which no longer depend on the viscoelastic effects. At this order, the problem simplifies to the rotation of a rod in a Newtonian fluid at low Reynolds number, which is known to produce an azimuthal vortex flow profile described by

(4.6) \begin{equation} \boldsymbol {U}^{(0)} = \frac {1}{R}\boldsymbol {e}_{\theta }, \quad \mathcal{P}^{(0)} = 0. \end{equation}

With this leading-order velocity profile, and (4.2), we can write down the leading-order polymeric stress:

(4.7) \begin{equation} {\boldsymbol{\varSigma }_p^{(0)} = 2\beta _p\boldsymbol{E}^{(0)}= -\frac {2\beta _p}{R^2}(\boldsymbol {e}_{r}\boldsymbol {e}_{\theta } + \boldsymbol {e}_{\theta }\boldsymbol {e}_{r})}. \end{equation}

Notice that (4.7) does not satisfy the condition $\boldsymbol{\varSigma }_p^{(0)} = \boldsymbol {0}$ at the initial time; we will resolve this point in §5.

4.2. First-order correction to the solution

We now proceed to determine the first-order correction to the shape of the interface. From (4.7), we observe that the leading-order polymeric stress appears as the rate-of-strain tensor from the leading-order vortical flow. This tensor has zero diagonal elements, resulting in no normal stress differences, which are typically responsible for the rod-climbing effect. To reveal the ‘hoop stress’ that leads to the rod-climbing effect, we solve for the first-order correction to the polymeric stress. Considering again the evolution of the polymeric stress equation (4.1), and using (4.6) and (4.7), we obtain

(4.8) \begin{align} \boldsymbol{\varSigma }_p^{(1)}&= {-\boldsymbol {U}^{(0)}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\varSigma }_p^{(0)}} + { \boldsymbol{\varSigma }_p^{(0)}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol {U}^{(0)}+(\boldsymbol{\nabla }\boldsymbol {U}^{(0)})^{\textrm T}\boldsymbol{\cdot }\boldsymbol{\varSigma }_p^{(0)} }- {\frac {\alpha }{\beta _p}\boldsymbol{\varSigma }_p^{(0)}\boldsymbol{\cdot }\boldsymbol{\varSigma }_p^{(0)}}+2\beta _p\boldsymbol{E}^{(1)} \notag \\ &= {-\frac {4\beta _p}{R^4}(\boldsymbol {e}_{r}\boldsymbol {e}_{r} - \boldsymbol {e}_{\theta }\boldsymbol {e}_{\theta })} + {\frac {4\beta _p}{R^4}(\boldsymbol {e}_{r}\boldsymbol {e}_{r} + \boldsymbol {e}_{\theta }\boldsymbol {e}_{\theta })} -{\frac {4\beta _p \alpha }{R^4}(\boldsymbol {e}_{r}\boldsymbol {e}_{r} + \boldsymbol {e}_{\theta }\boldsymbol {e}_{\theta })}+ 2\beta _p\boldsymbol{E}^{(1)} \notag \\ &= \frac {4\beta _p}{R^4}\left ((2-\alpha )\boldsymbol {e}_{\theta }\boldsymbol {e}_{\theta } - \alpha \boldsymbol {e}_r \boldsymbol {e}_r\right )+ 2\beta _p\boldsymbol{E}^{(1)}. \end{align}

As in §4.1, we substitute (4.8) into the momentum equation (2.17) at first order, to find

(4.9) \begin{align} \boldsymbol {0} &= -\boldsymbol{\nabla }\mathcal{P}^{(1)} + \beta _s \, {\nabla} ^2 \boldsymbol {U}^{(1)} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\varSigma }_p^{(1)} \notag \\ &= -\boldsymbol{\nabla }\mathcal{P}^{(1)} + {\nabla} ^2 \boldsymbol {U}^{(1)} + 4\beta _p\,\boldsymbol{\nabla }\boldsymbol{\cdot }\left (\frac {2-\alpha }{R^4}\boldsymbol {e}_{\theta }\boldsymbol {e}_{\theta } - \frac {\alpha }{R^4} \boldsymbol {e}_r \boldsymbol {e}_r\right )\notag \\ &= -\boldsymbol{\nabla }\mathcal{P}^{(1)} + {\nabla} ^2 \boldsymbol {U}^{(1)} -\frac {8\beta _p(1-2\alpha )}{R^5}\boldsymbol {e}_r. \end{align}

Although the anisotropic hydrodynamic drag in the Giesekus model contributes to both shear-thinning behaviour and a non-zero second normal stress difference, we observe from (4.8) and (4.9) that only the normal stress is responsible for the rod climbing at this order. In particular, the second normal stress difference, represented by the terms containing the mobility factor $\alpha$ , decreases the magnitude of the hoop stress in the final term of (4.9), thereby reducing the rise height, while the shear-thinning behaviour of the Giesekus fluid affects the interface shape at higher-order corrections.

Observe from (4.8) that $\boldsymbol {e}_z \boldsymbol{\cdot }\boldsymbol{\varSigma }_p^{(1)}=2\beta _p\boldsymbol {e}_z \boldsymbol{\cdot }\boldsymbol{E}^{(1)}$ . Thus we can rewrite the stress boundary condition (3.11b ) at first order in $\textit{Wi}$ as

(4.10) \begin{equation} \boldsymbol {e}_z \boldsymbol{\cdot }\left . \left (-\mathcal{P}^{(1)}\boldsymbol {I} + 2 \boldsymbol{E}^{(1)} \right )\right |_{Z = 0} = {-{G}\left (H^{(2)}-\frac {1}{\textit{BR}}\,\frac {{\textrm d}}{{\textrm d} R}\left (R\frac {{\textrm d}H^{(2)}}{{\textrm d} R}\right )\right )\boldsymbol {e}_z.} \end{equation}

An important observation in solving the problem at first order is that all boundary conditions for the velocity field are constrained to be zero. Specifically, the no-slip boundary condition (2.15a ) imposes $\boldsymbol {U}^{(1)} = \boldsymbol {0}$ at $R = 1$ , the far-field boundary condition (2.15e ) enforces $\boldsymbol {U}^{(1)} = \boldsymbol {0}$ as $R \rightarrow \infty$ , and the kinematic boundary condition (4.3) gives $U_Z^{(1)} \big|_{Z = 0} \big. = 0$ . These features suggest $\boldsymbol {U}^{(1)} = \boldsymbol {0}$ , reducing the problem to solving for the pressure and interface height using (4.9) and (4.10) with the corresponding boundary conditions. Specifically, (4.9) implies that $\mathcal{P}^{(1)}(R)$ , allowing us to simplify (4.10) to $\mathcal{P}^{(1)} = {G} (H^{(2)}-({1}/{\textit{BR}})( {\textrm d}/{{\textrm d} R}) (R\,{{\textrm d} H^{(2)}}/{{\textrm d} R} ) )$ . The far-field boundary condition (2.15e ) then implies that the pressure decays to zero as the interface flattens, i.e. $\mathcal{P}^{(1)}(R \rightarrow \infty ) = 0$ .

We can now integrate the momentum equation (4.9) once with respect to $R$ to find the pressure:

(4.11) \begin{equation} \mathcal{P}^{(1)} (R)= \frac {2\beta _p(1-2\alpha )}{R^4}. \end{equation}

On the other hand, the free-surface shape $H^{(2)}$ is governed by the differential equation

(4.12) \begin{equation} H^{(2)} - \frac {1}{\textit{BR}} \frac {\textrm d}{{\textrm d} R} \left (R \frac {{\textrm d}H^{(2)}}{{\textrm d} R} \right ) = \frac {2\beta _p(1-2\alpha )}{{G}R^4}. \end{equation}

Equation (4.12) is subject to the contact angle boundary condition (2.15b ) and the far-field condition (2.15e ), yielding the solution

(4.13) \begin{eqnarray} H^{(2)}_{\textit{steady}}(R) \!\!&=&\!\! \frac {2\beta _p(1-2\alpha )B}{{G}}\left(I_1 \big(\sqrt {B} \big)K_1^{-1}\big(\sqrt {B} \big) \left(\int _{1}^{\infty }s^{-3}K_0 \big(\sqrt {B}s \big) {\textrm{d}}s \right)K_0 \big(\sqrt {B}R \big) \notag \right. \nonumber\\ &&\!\! \left. +\, I_0 \big(\sqrt {B}R \big)\int _{R}^{\infty }s^{-3}K_0 \big(\sqrt {B}s \big) {\textrm d}s + K_0 \big(\sqrt {B}R \big)\int _{1}^{R}s^{-3}I_0 \big(\sqrt {B}s \big) {\textrm{d}}s \right), \nonumber\\ \end{eqnarray}

where $I_\nu$ and $K_\nu$ are the first- and second-kind modified Bessel functions of order $\nu$ , respectively. In figure 2, we present the scaled steady interface shape $H(R)\,\mathcal{G}/[2\beta _p(1-2\alpha )\,\textit{Wi}]$ around an infinitely long rotating rod in a Giesekus fluid for different Bond numbers.

Figure 2. Steady interface shape around an infinitely long rotating rod in a Giesekus fluid, with inertial effects absent, displayed for different Bond numbers $B = 1, 5, 25$ .

Although it is not straightforward to visualise the height profile in (4.13), we can infer from (4.12) that for $R \gg 1$ , the dominant balance is between the gravitational term and the viscoelastic term:

(4.14) \begin{equation} H_{\textit{steady}}(R) \approx \frac {2\beta _p(1-2\alpha )\,\textit{Wi}}{\mathcal{G}R^4} \quad (R \gg 1). \end{equation}

A result similar to (4.14) was obtained previously for the second-order fluid model (Joseph & Fosdick Reference Joseph and Fosdick1973; More et al. Reference More, Patterson, Pashkovski and McKinley2023), which is known to agree with the Giesekus model at steady state for small Weissenberg numbers with $\alpha = 0$ . However, in earlier studies, domain perturbation methods were developed using dimensional rotation speeds, without explicitly addressing the impact of individual dimensionless parameters. In addition, instead of providing solutions in explicit form, a method of successive approximations was used to estimate the final steady interface height. In §§5 and6, we increase the complexity by introducing time-dependent and inertia terms, while maintaining the same dimensionless analytical structure established here.

5.1. Leading-order solution

In §4, we examined the problem at an intermediate time scale $T = O(1)$ , and observed that the leading-order interface shape (4.13) is time-independent. This result occurs because at $T = O(1)$ , substituting the polymeric stress expansion (3.7c ) into (4.1) has corrections (4.7) and (4.8) that do not satisfy the initial condition $\boldsymbol{\varSigma }_p = \boldsymbol {0}$ at $T = 0$ . Physically, these features suggest that no transient changes in the polymeric stress are observed at intermediate times, as the dynamics of these changes occurs on a much shorter time scale. Therefore, we introduce a stretched time variable to balance the time derivative term on the left-hand side of (4.1), ${\textit{Wi}}\, ({\partial \boldsymbol{\varSigma} _p}/{\partial T})$ , with the right-hand side of (4.1), $-\boldsymbol{\varSigma }_p + 2\beta _p\boldsymbol{E}$ , by defining

(5.1) \begin{equation} \tau = \frac {T}{\textit{Wi}}. \end{equation}

In terms of physical variables, time has now been scaled by the relaxation time. More precisely, we rescale time by interpreting $T = O(\textit{Wi})$ as $\tau = O(1)$ .

The evolution equation for the polymeric stress (4.1) in terms of the stretched time variable $\tau$ now becomes

(5.2) \begin{equation} \textit{Wi}\left ( \boldsymbol {U} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\varSigma }_p - \boldsymbol{\varSigma }_p \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol {U} - (\boldsymbol{\nabla }\boldsymbol {U})^{\textrm {T}}\boldsymbol{\cdot }\boldsymbol{\varSigma }_p + \frac {\alpha }{\beta _p}\boldsymbol{\varSigma }_p \boldsymbol{\cdot }\boldsymbol{\varSigma }_p \right )= -\left (\frac {\partial \boldsymbol{\varSigma }_p}{\partial \tau } +\boldsymbol{\varSigma }_p - 2\beta _p\boldsymbol{E} \right ), \end{equation}

where we distinguish the polymeric stress in the inner layer with $\boldsymbol{\varSigma }_p(R,Z,\tau )$ .

Following the approach in §4.1, we express the leading-order polymeric stress in terms of the leading-order velocity profile using (5.2) as

(5.3) \begin{equation} \frac {\partial \boldsymbol{\varSigma }_p^{(0)}}{\partial \tau } + \boldsymbol{\varSigma }_p^{(0)} = 2\beta _p\boldsymbol{E}^{(0)} \quad \Longrightarrow \quad \boldsymbol{\varSigma }_p^{(0)} = 2\beta _p\, {\rm e}^{-\tau }\int _{0}^{\tau }\boldsymbol{E}^{(0)}\, {\rm e}^{s}\, {\textrm d}s, \end{equation}

where we have used the condition $\boldsymbol{\varSigma }_p^{(0)} = \boldsymbol {0}$ at $\tau = 0$ . We then substitute (5.3) into the leading-order momentum equation (2.17) to arrive at

(5.4) \begin{equation} \boldsymbol {0} = -\boldsymbol{\nabla }\mathcal{P}^{(0)} + \beta _s\, {\nabla} ^2 \boldsymbol {U}^{(0)} + \beta _p\, {\rm e}^{-\tau }\int _0^{\tau } {\nabla} ^2 \boldsymbol {U}^{(0)}\, {\rm e}^s\, {\textrm d}s. \end{equation}

Equation (5.4) is solved subject to the no-slip boundary condition (2.15a ), the far-field boundary condition (2.15e ), the stress boundary condition (3.11a ), and the kinematic boundary condition (2.15d ).

It is helpful to introduce a linear operator acting on the time-dependent vector field that resembles the right-hand side of (5.4) by defining

(5.5) \begin{equation} {\mathcal{L}}(\boldsymbol {F}(\boldsymbol {X},\tau )) = \beta _s\, \boldsymbol {F}(\boldsymbol {X},\tau ) + \beta _p\, {\rm e}^{-\tau } \int _{0}^{\tau } \boldsymbol {F}(\boldsymbol {X},s)\, {\rm e}^{s}\, {\textrm d}s, \end{equation}

so that (5.4) becomes $\boldsymbol {0} = -\boldsymbol{\nabla }\mathcal{P}^{(0)} + {\mathcal{L}} ({\nabla} ^2 \boldsymbol {U}^{(0)} )$ . With this notation, we simplify the stress boundary condition (3.11a ) by substituting (5.3) to obtain

(5.6) \begin{equation} 2{\mathcal{L}}\left (\boldsymbol {e}_z \boldsymbol{\cdot }\boldsymbol{E}^{(0)}\right|_{Z = 0}\left . \vphantom{\boldsymbol{\textit{E}}^{(0)}}\right ) = \left . \vphantom{\boldsymbol{\textit{E}}^{(0)}}\mathcal{P}^{(0)}\right|_{Z = 0} \boldsymbol {e}_{z}. \end{equation}

The kinematic boundary condition (2.15d ), however, requires more careful examination, as it involves a time derivative that is crucial to the analysis in this section. Nevertheless, at leading order, we still find $U_Z^{(0)} \big |_{Z = 0} = 0$ . Since the flow begins at rest with a flat interface, this absence of upward velocity suggests that the interface will remain flat throughout. Therefore, we aim to find a solution where $ U_Z^{(0)}$ is identically zero, and the velocity field $ \boldsymbol {U}^{(0)} = \boldsymbol {U}^{(0)}(R, \tau )$ is independent of $Z$ . The continuity equation (2.14a ), expressed in cylindrical coordinates as $ ({1}/{R}) ({\partial (R U_R^{(0)})}/{\partial R}) + ({\partial U_Z^{(0)}}/{\partial Z}) = 0$ , then forces $ U_R^{(0)} = c(\tau )/{R}$ for some function $ c=c(\tau )$ . Applying the no-slip boundary condition (2.15a ) implies $c(\tau ) = 0$ , resulting in a purely azimuthal velocity field, $\boldsymbol {U}^{(0)} = U_{\theta }^{(0)}(R, \tau )\, \boldsymbol {e}_\theta$ .

Continuing this approach, the momentum equation (2.17), at leading order in the $r$ - and $z$ -directions, yields ${\partial \mathcal{P}^{(0)}}/{\partial R} = {\partial \mathcal{P}^{(0)}}/{\partial Z} = 0$ , so that $ \mathcal{P}^{(0)}(\tau )$ only, while the stress boundary condition (5.6) in the $z$ -direction further implies $\mathcal{P}^{(0)} \equiv 0$ . Thus under the flat interface assumption, we have deduced that the pressure is absent from the leading-order solution. The momentum equation (2.17) and the stress boundary condition (5.6) simplify to

(5.7) \begin{equation} {\mathcal{L}}\left ({\nabla} ^{2}\boldsymbol {U}^{(0)}\right ) = {\mathcal{L}}\left (\boldsymbol {e}_z \boldsymbol{\cdot }\boldsymbol{E}^{(0)}\right|_{Z = 0}\left . \vphantom{\boldsymbol{\textit{E}}^{(0)}}\right ) = \boldsymbol {0} \quad \Longrightarrow \quad {\nabla} ^{2}\boldsymbol {U}^{(0)} = \left .\boldsymbol {e}_z \boldsymbol{\cdot }\boldsymbol{E}^{(0)}\right|_{Z = 0} = \boldsymbol {0}, \end{equation}

where we have used the fact that the kernel of the linear operator ${\mathcal{L}}$ is null (see Appendix A).

We observe that (5.7) has the same structure as the leading-order solution analysed in §4.1; in the weak viscoelastic limit, the leading-order solution is effectively Newtonian. Although the momentum equation (5.4) contains a time-dependent integral, the solution remains time-independent. Specifically, the velocity field retains the form of a vortex solution as in (4.6):

(5.8) \begin{equation} \boldsymbol {U}^{(0)} = \frac {1}{R}\boldsymbol {e}_{\theta }. \end{equation}

With this leading-order velocity profile, we can express the leading-order polymeric stress, using (5.3), as

(5.9) \begin{equation} \boldsymbol{\varSigma }_p^{(0)} = 2\beta _p\, {\rm e}^{-\tau }\int _{0}^{\tau }\boldsymbol{E}^{(0)}\, {\rm e}^{s}\, {\textrm d}s = 2\beta _p(1-{\rm e}^{-\tau })\boldsymbol{E}^{(0)} = -\frac {2\beta_p(1- {\rm e}^{-\tau })}{R^2}(\boldsymbol {e}_{r}\boldsymbol {e}_{\theta } + \boldsymbol {e}_{\theta }\boldsymbol {e}_{r}). \end{equation}

Starting from zero at $\tau =0$ , $ \boldsymbol{\varSigma }_p^{(0)}$ increases monotonically towards the steady state (4.7), approaching it with an exponential decay rate on the time scale of $ \tau$ . This behaviour confirms our remark at the start of this subsection that transient changes in the polymeric stress occur on a shorter time scale when $ T = O(\textit{Wi})$ .

5.2. First-order correction solution

Hitherto, we have derived that for an inner layer in time where $ T =O(\textit{Wi})$ , the leading-order flow remains the same vortex solution as in the outer layer over times $ T = O(1)$ . Although the velocity fields are the same at the leading order, the leading-order polymeric stresses differ: in the inner layer, $\boldsymbol{\varSigma }_p^{(0)}$ varies with time, whereas it does not in the outer layer. Similarly, we anticipate the polymeric stress arising from $\boldsymbol{\varSigma }_p^{(1)}$ , which defined the interface shape in §4.2, to also be time-dependent, enabling us to capture the transient interface shape. We can determine the time evolution of $ \boldsymbol{\varSigma }_p^{(1)}$ using calculations similar to those in (4.8):

(5.10) \begin{align} \frac {\partial \boldsymbol{\varSigma }_p^{(1)}}{\partial \tau } + \boldsymbol{\varSigma }_p^{(1)}&= {-\boldsymbol {U}^{(0)}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\varSigma }_p^{(0)}} + { \boldsymbol{\varSigma }_p^{(0)}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol {U}^{(0)}+(\boldsymbol{\nabla }\boldsymbol {U}^{(0)})^{\textrm T}\boldsymbol{\cdot }\boldsymbol{\varSigma }_p^{(0)} }\nonumber\\ & \quad -\, {\frac {\alpha }{\beta _p}\boldsymbol{\varSigma }_p^{(0)}\boldsymbol{\cdot }\boldsymbol{\varSigma }_p^{(0)}}+2\beta _p\boldsymbol{E}^{(1)} \notag \\ &= {-\frac {4\beta _p(1-{\rm e}^{-\tau })}{R^4}(\boldsymbol {e}_{r}\boldsymbol {e}_{r} - \boldsymbol {e}_{\theta }\boldsymbol {e}_{\theta })} + {\frac {4\beta _p(1-{\rm e}^{-\tau })}{R^4}(\boldsymbol {e}_{r}\boldsymbol {e}_{r} + \boldsymbol {e}_{\theta }\boldsymbol {e}_{\theta })} \notag \\ &\quad {}-{\frac {4\beta _p \alpha (1-{\rm e}^{-\tau })^2 }{R^4}(\boldsymbol {e}_{r}\boldsymbol {e}_{r} + \boldsymbol {e}_{\theta }\boldsymbol {e}_{\theta })}+ 2\beta _p\boldsymbol{E}^{(1)} \notag \\ &= \frac {4\beta _p}{R^4}\big (\big (2(1-{\rm e}^{-\tau })-\alpha (1-{\rm e}^{-\tau })^2\big )\boldsymbol {e}_{\theta }\boldsymbol {e}_{\theta } - \alpha (1-{\rm e}^{-\tau })^2\boldsymbol {e}_r \boldsymbol {e}_r \big ) \!+\! 2\beta _p\boldsymbol{E}^{(1)}. \end{align}

Multiplying both sides by the integrating factor ${\rm e}^{\tau }$ , we find $\boldsymbol{\varSigma }_p^{(1)}$ to be

(5.11) \begin{align} \boldsymbol{\varSigma }_p^{(1)} &= 2\beta _p\, {\rm e}^{-\tau }\int _{0}^{\tau }\boldsymbol{E}^{(1)}\, {\rm e}^s\, {\textrm d}s + \frac {4\beta _p\left(2-\alpha -2\, {\rm e}^{-\tau }-(2-2\alpha )\tau {\rm e}^{-\tau }+\alpha {\rm e}^{-2\tau}\right)}{R^4}\boldsymbol {e}_{\theta }\boldsymbol {e}_{\theta } \nonumber \\ &\quad {}- \frac {4\beta _p\alpha \left(1-2\tau {\rm e}^{-\tau }-{\rm e}^{-2\tau}\right)}{R^4}\boldsymbol {e}_r \boldsymbol {e}_r, \end{align}

where we have used the initial condition $\boldsymbol{\varSigma }_p^{(1)} = \boldsymbol {0}$ at $\tau = 0$ to write the lower limit of the integral. We then substitute (5.11) into the first-order momentum equation (2.17) to obtain

(5.12) \begin{align} \boldsymbol {0} &= -\boldsymbol{\nabla }\mathcal{P}^{(1)} + \beta _s \, {\nabla} ^2 \boldsymbol {U}^{(1)} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\varSigma }_p^{(1)} \notag \\&= -\boldsymbol{\nabla }\mathcal{P}^{(1)} + {\mathcal{L}}\big ({\nabla} ^2 \boldsymbol {U}^{(1)}\big ) -\frac {8\beta _p\left(1-2\alpha -{\rm e}^{-\tau }-(1-4\alpha )\tau {\rm e}^{-\tau }+2\alpha {\rm e}^{-2\tau}\right)}{R^5}\boldsymbol {e}_r. \end{align}

From (5.12), we observe that as the suspended polymers begin to stretch from their relaxed state, the elastic response generates a time-dependent force per volume that acts purely in the radial direction and is zero initially. This hoop stress is expected to drive the fluid radially inwards towards the rod, which then, by conservation of mass, would lead the interface to rise, resulting in the rod-climbing effect. Notably, the force induced by the hoop stress is conservative. Thus it is convenient to define a new pressure variable that includes the polymeric stress, by setting

(5.13) \begin{equation} {\textrm{P}^{(1)} = \mathcal{P}^{(1)} - \frac {2\beta _p\left(1-2\alpha -\textrm{e}^{-\tau }-(1-4\alpha )\tau {\rm e}^{-\tau }+2\alpha {\rm e}^{-2\tau }\right)}{R^4},} \end{equation}

which allows (5.12) to be rewritten as

(5.14) \begin{equation} 0 = -\boldsymbol{\nabla }{{\rm P}}^{(1)} + {\mathcal{L}}\left ({\nabla} ^2 \boldsymbol {U}^{(1)}\right ). \end{equation}

In addition to the usual no-slip and far-field boundary conditions, we anticipate that the free-surface boundary conditions at the interface are particularly crucial for this first-order solution. We substitute (5.11) into the first-order normal stress boundary condition (3.11b ), and express the equations using the newly defined pressure variable, leading to the result

(5.15) \begin{align} \kern-4pt -{{\rm P}}^{(1)}\Big|_{Z = 0} \Big. \boldsymbol {e}_z+ 2{\mathcal{L}}\vphantom{\boldsymbol{\textit{E}}^{(1)}} \left (\boldsymbol {e}_z \boldsymbol{\cdot }\boldsymbol{E}^{(1)}\right|_{\textit{Z} = 0}\left . \vphantom{\boldsymbol{\textit{E}}^{(1)}}\!\right ) & = \left(\frac {2\beta _p\left(1-2\alpha -{\rm e}^{-\tau }-(1-4\alpha )\tau {\rm e}^{-\tau }+2\alpha {\rm e}^{-2\tau}\right)}{R^4} \notag \right. \\ & \left.\ \ \ -\, {G}\left (H^{(2)}-\frac {1}{\textit{BR}}\,\frac {\partial }{\partial R}\left (R\frac {\partial H^{(2)}}{\partial R}\right )\right )\right)\boldsymbol {e}_z. \end{align}

On the other hand, recall from (3.8) that $H =O(\textit{Wi}^2)$ . This implies that when $T = O(\textit{Wi})$ , we have $({\partial H}/{\partial T}) = O(\textit{Wi})$ . In other words, temporal variations in the interface shape will emerge in this first-order solution. Unlike for the steady interface shape (4.3), the corresponding kinematic boundary condition (2.15d ) at this order becomes

(5.16) \begin{equation} U_Z^{(1)}\big {|}_{Z=0} = \frac {\partial H^{(2)}}{\partial \tau }. \end{equation}

It is difficult to obtain exact analytical expressions for the velocity field that satisfies (5.14), (5.15) and (5.16) for arbitrary values of ${G} = {O}(1)$ . However, when ${G} \gg 1$ , the method of dominant balance enables progress in determining the interface shape, with the dominant terms in this limit being the terms on the right-hand side of (5.15),

(5.17) \begin{equation} H^{(2)}_{\textit{transient}}(R,\tau ) \sim \left (\frac {1-2\alpha -{\rm e}^{-\tau }-(1-4\alpha )\tau {\rm e}^{-\tau }+2\alpha {\rm e}^{-2\tau }}{1-2\alpha }\right )H^{(2)}_{\textit{steady}}(R) = O\left ({G}^{-1}\right ), \end{equation}

with $H^{(2)}_{\textit{steady}}(R)$ given in (4.13). Note that (5.17) is regular at $\alpha = 1/2$ because the steady-state solution (4.13) is also proportional to $1-2\alpha .$ To confirm the consistency of this dominant balance, we need to assess the relative magnitudes of the pressure and velocities in terms of the parameter ${G}$ . From (5.17), the kinematic boundary condition (5.16) implies that $U_Z^{(1)} = O ({G}^{-1} )$ . According to the continuity equation (2.14a ), $U_R^{(1)}$ must match the order of $U_Z^{(1)}$ , giving $\boldsymbol {U}^{(1)} = O ({G}^{-1} )$ . This, in turn, implies from the momentum equation (5.14) that ${P}^{(1)}$ is also of order $O ({G}^{-1} )$ . Consequently, both the pressure and the rate of strain on the left-hand side of (5.15) are of order $O ({G}^{-1} )$ , while the right-hand side remains of order $O(1)$ . Therefore, as ${G} \gg 1$ , so that ${G}^{-1}\ll 1$ , we conclude that the right-hand side forms a consistent dominant balance as desired. In §4, we have derived a steady interface shape (4.13) for $T = O(1)$ , whereas in this section, we obtain a transient interface shape (5.17) for $T =O(\textit{Wi})$ . We note that

(5.18) \begin{equation} \lim _{\tau \rightarrow \infty } H_{\textit{transient}}(R,\tau ) = H_{\textit{steady}}(R), \end{equation}

thus confirming that the interface dynamics occurs on a short time scale, with $T = O(\textit{Wi})$ , and that the transient interface shape approaches the steady interface shape at long times.

Figure 3(a,b) show the time evolutions of the scaled interface shape $H(R,\tau )\,\mathcal{G}/[2\beta _p(1-2\alpha )\,\textit{Wi}]$ and the scaled climbing height of a Giesekus fluid on the rotating rod, $H(R=1,\tau )\,\mathcal{G}/[2\beta _p(1-2\alpha )\,\textit{Wi}]$ , given in (5.17), with $B = 5$ and $\alpha = 0.25$ . We observe that the transient interface height reaches a steady state at $\tau \approx 6$ , with the interface shape becoming nearly flat for $R \gtrapprox 3$ .

Figure 3. Time-varying interface shape for a Giesekus fluid. ( $a$ ) Time evolution of the interface shape around an infinitely long rotating rod in a Giesekus fluid, with inertial effects absent, shown at times $\tau = 0, 1, 3, 10$ . ( $b$ ) Time evolution of the climbing height of a Giesekus fluid on the rotating rod, with inertial effects absent, evaluated at $R = 1$ . All calculations were performed using $B = 5$ and $\alpha = 0.25$ .

6.1. Leading-order solution

As derived in §5, the transient interface shape (5.17) remains positive, owing to the hoop stress in the viscoelastic fluids, which drives the fluid radially inwards and is responsible for rod climbing. However, introducing a small inertial effect may change the dynamics, as the resulting centrifugal force would push the fluid radially outwards. This competition between forces leads us to investigate how finite, small inertia may impact the interface height when both effects are present. In this section, we continue to use the stretched time scale $\tau$ to examine the transient interface dynamics.

Incorporating inertial effects into our analysis, the momentum equation (2.17) is

(6.1) \begin{equation} {Re} \left (\frac {1}{\textit{Wi}}\frac {\partial \boldsymbol {U}}{\partial \tau } + \boldsymbol {U}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol {U}\right ) = \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{{\varSigma }}. \end{equation}

Rather than focusing our analysis on the Reynolds number, we find it more useful to base our analysis on the elasticity number $El$ , which represents the ratio of elastic stress to inertial stress (Denn & Porteous Reference Denn and Porteous1971):

(6.2) \begin{equation} El = \frac {\textit{Wi}}{{Re}}. \end{equation}

Consequently, the momentum equation (6.1) can be expressed as

(6.3) \begin{equation} El^{-1} \left (\frac {\partial \boldsymbol {U}}{\partial \tau } + {\textit{Wi}} \boldsymbol {U}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol {U}\right ) = \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{{\varSigma }}. \end{equation}

Previously, we derived the interface shape for $El \rightarrow \infty$ , and now aim to extend our analysis to cases where $El$ is large. To isolate the nonlinear term $\boldsymbol {U} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol {U}$ from the leading-order velocity profile, we assume that $El$ is at least $O(1)$ , indicating small but finite inertial effects with ${Re} \lesssim \textit{Wi}$ . When $El$ is finite, the time-derivative $ {\partial \boldsymbol {U}}/{\partial \tau }$ of the velocity field reappears, in contrast to the low-Reynolds-number approximation in §§4 and5. To ensure that the problem is well-posed, we specify the velocity field at the initial time:

(6.4) \begin{equation} \boldsymbol {U} = \boldsymbol {0} \quad \text{at} \quad \tau = 0. \end{equation}

Besides the additional terms on the left-hand side of the momentum equation (6.3) and the new initial condition (6.4), the remaining constitutive equations and boundary conditions are identical to those in §5. Specifically, the leading-order polymeric stress $\boldsymbol{\varSigma }_p^{(0)}$ can be expressed in terms of the leading-order rate-of-strain tensor, as shown in (5.3). This expression, combined with the momentum equation (6.3), allows us to derive the differential equation for the leading-order velocity field:

(6.5) \begin{equation} El^{-1}\,\frac {\partial \boldsymbol {U}^{(0)}}{\partial \tau } = -\boldsymbol{\nabla }\mathcal{P}^{(0)} + {\mathcal{L}}\left ({\nabla} ^2 \boldsymbol {U}^{(0)}\right ). \end{equation}

Following the approach in §5.1, we seek a solution where the leading-order interface is flat. Using the same reasoning, we deduce that $\mathcal{P}^{(0)} = 0$ and the velocity field is purely azimuthal, $ \boldsymbol {U}^{(0)} = U^{(0)}(R,\tau )\, \boldsymbol {e}_{\theta }$ . This allows us to reduce the vector differential equation (6.5) to a partial differential equation for an unknown scalar variable:

(6.6) \begin{equation} El^{-1}\,\frac {\partial U^{(0)}}{\partial \tau } = {\mathcal{L}}\left (\frac {1}{R}\, \frac {\partial }{\partial R} \left (R\frac {\partial U^{(0)}}{\partial R}\right ) - \frac {U^{(0)}}{R^2}\right ). \end{equation}

Equation (6.6) is formulated on the semi-infinite domain $R \in [1, \infty )$ subject to the no-slip boundary condition (2.15a ), the far-field boundary condition (2.15e ), and the initial condition (6.4), i.e.

(6.7) \begin{equation} U^{(0)}(1,\tau ) = 1, \quad U^{(0)}(\infty ,\tau ) = 0, \quad U^{(0)}(R,0) = 0. \end{equation}

Since the vortex solution ${1}/{R}$ represents the steady-state solution of (6.6) along with the specified no-slip and far-field boundary conditions, we decompose $U^{(0)}$ into a steady-state component and a transient component:

(6.8) \begin{equation} U^{(0)}(R,\tau ) = \frac {1}{R} + \tilde {U}^{(0)}(R,\tau ). \end{equation}

As (6.6) is linear, $\tilde {U}^{(0)}$ also satisfies (6.6) but with homogeneous spatial boundary conditions:

(6.9) \begin{equation} \tilde {U}^{(0)}(1,\tau ) = 0, \quad \tilde {U}^{(0)}(\infty ,\tau ) = 0, \quad \tilde {U}^{(0)}(R,0) = -\frac {1}{R}. \end{equation}

We approach this problem using a Weber transform $\mathcal{W}$ of order one (Watson Reference Watson1966; Piessens Reference Piessens2000), an analogue of the Hankel transform on the semi-infinite domain $[1, \infty )$ , by introducing

(6.10) \begin{equation} V[\tau ; k] = \mathcal{W}\big \{\tilde {U}^{(0)}(R, \tau )\big \} = \int _{1}^{\infty } R\, \tilde {U}^{(0)}(R, \tau ) \, \textit{Z}_{1}(k, R) \, {\textrm d}R, \end{equation}

with the kernel $\textit{Z}_{1}(k, R)$ given by

(6.11) \begin{equation} \textit{Z}_{1}(k, R) = \textit{J}_{1}(kR)\, \textit{Y}_{1}(k) - \textit{Y}_{1}(kR)\, \textit{J}_{1}(k), \end{equation}

where $\textit{J}_1$ and $\textit{Y}_1$ are the first- and second-kind Bessel functions of order one, respectively. This transformation enables us to simplify the Bessel differential operator of order one acting on $\tilde {U}^{(0)}$ on the right-hand side of (6.6), provided that $\tilde {U}^{(0)}$ decays to zero at infinity:

(6.12)

Note that the operator $\mathcal{W}$ acts on the $R$ variable, while ${\mathcal{L}}$ acts on $\tau$ , which implies that the two operators commute: $\mathcal{W} \circ {\mathcal{L}} = {\mathcal{L}} \circ \mathcal{W}$ . We apply the Weber transform to (6.6), yielding

(6.13) \begin{equation} El^{-1}\,\frac {\partial V[\tau ; k]}{\partial \tau } = -k^2\,{\mathcal{L}}\left (V[\tau ;k]\right ) = -k^2\left (\beta _s\, V[t;k] + \beta _p\, {\rm e}^{-\tau }\int _{0}^{\tau } V[s;k]\, {\rm e}^s \, {\textrm d}s\right ). \end{equation}

In this transformed variable, the initial condition becomes

(6.14) \begin{align} V[0;k] = \mathcal{W}\left \{\tilde {U}^{(0)}(R,0)\right \} = -\int _{1}^{\infty } \textit{Z}_{1}(k,R) \, {\textrm d}R \notag \\ &= \textit{J}_{1}(k)\int _{1}^{\infty } \textit{Y}_{1}(kR) \, {\textrm d}R - \textit{Y}_{1}(k)\int _{1}^{\infty }\textit{J}_{1}(kR) \, {\textrm d}R \notag \\ &= \frac {1}{k}\left (\textit{J}_{1}(k)\,\textit{Y}_{0}(k) - \textit{J}_0(k)\,\textit{Y}_{1}(k)\right ) = \frac {2}{\unicode{x03C0} k^2}, \end{align}

where the last equality follows from the Wronskian of the zero-order Bessel differential equation. To solve (6.13), we first multiply both sides by ${\rm e}^{\tau }$ and then differentiate with respect to $\tau$ , to obtain

(6.15) \begin{equation} El^{-1}\left (\frac {\partial ^2 V}{\partial \tau ^2} + \frac {\partial V}{\partial \tau }\right ) = -k^2\left (\beta _s \frac {\partial V}{\partial \tau } + V\right ), \end{equation}

where we have used the relation $\beta _s + \beta _p = 1$ . Equation (6.15) is a constant-coefficient second-order differential equation in time, which has solutions of the form

(6.16) \begin{equation} V[\tau ;k] = C_1(k)\, {\rm e}^{\lambda _{+} \tau } + C_2(k)\, {\rm e}^{\lambda _{-} \tau }, \end{equation}

where $\lambda _{+}$ and $\lambda _{-}$ are roots of the equations

(6.17) \begin{equation} \lambda ^2 + \left (1 + k^2 El \beta _s\right )\lambda + k^2 El = 0. \end{equation}

In particular,

(6.18) \begin{equation} \lambda _{\pm } = \frac {-1-k^2 El \beta _s \pm \sqrt {(1+k^2 El \beta _s)^2 - 4k^2 El}}{2}. \end{equation}

Although $\lambda _{\pm }$ can be complex depending on the value of $k$ , their real parts are always negative, except at $k=0$ , where one root is zero and the other is $-1$ . In addition to the initial condition (6.14), evaluating (6.13) at $\tau = 0$ provides an initial condition for ${\partial V}/{\partial \tau }$ :

(6.19) \begin{equation} \frac {\partial V}{\partial \tau }[0;k] = -k^2 El \beta _s\,V[0;k] = -\frac {2 El\beta _s}{\unicode{x03C0} }. \end{equation}

Using (6.14), (6.16) and (6.19), and simplifying with (6.17), we find

(6.20) \begin{equation} V[\tau ;k] = \frac {2}{\unicode{x03C0} }\left (\frac {\left (1+\lambda _{+}\right ){\rm e}^{\lambda _{+} \tau } - \left (1 + \lambda _{-}\right ) {\rm e}^{\lambda _{-}\tau }}{k^2(\lambda _{+} - \lambda _{-})}\right ), \end{equation}

where $\lambda _{+}$ and $\lambda _{-}$ are functions of $k$ , $El$ and $\beta _s$ . Finally, we take an inverse Weber transform to obtain

(6.21) \begin{align} \tilde {U}^{(0)}(R, \tau ) &= \int _{0}^{\infty } \frac {k\,V[\tau ;k]\,\textit{Z}_1(k,R)}{{J}_1^2(k) + {Y}_1^2(k)}\, {\textrm d}k \notag \\ &= \frac {2}{\unicode{x03C0} }\int _{0}^{\infty } \frac {\left (\left (1+\lambda _{+}\right ){\rm e}^{\lambda _{+} \tau } - \left (1 + \lambda _{-}\right ) {\rm e}^{\lambda _{-} \tau }\right ) \textit{Z}_1(k,R)}{ k(\lambda _{+} - \lambda _{-})\left ({J}_1^2(k) + {Y}_1^2(k)\right )}\, {\textrm d}k. \end{align}

Although proceeding with (6.21) analytically is difficult, it is straightforward to handle numerically. In figure 4, we show the time evolution of the leading-order velocity profile $U^{(0)}(R,\tau )$ for various values of $El$ . We observe that larger $El$ leads to a faster approach to a steady state in the transient velocity profile. Intuitively, when $El$ is large, the elastic effect ( $\textit{Wi}$ ) is larger than the inertial effect ( ${Re}$ ). A larger viscoelasticity accelerates the flow’s return to equilibrium by emphasising the elastic ‘spring-like’ properties over inertial resistance, thus pushing the flow to a steady state faster.

Figure 4. Time evolution of the leading-order velocity field $U^{(0)}(R, \tau )$ around an infinitely long rotating rod in a Giesekus fluid, incorporating small but finite inertial effects, shown at times $\tau = 0.1, 1, 10$ , for ( $a$ ) $El=0.2$ , ( $b$ ) $El=1$ , and ( $c$ ) $El=5$ . All calculations were performed using $\beta _p = 0.5$ .

6.2. First-order correction solution

In §§4.2 and5.2, we used the leading-order velocity field $\boldsymbol {U}^{(0)}$ to calculate the polymeric force $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\varSigma }_p^{(1)}$ , which is then incorporated into the pressure gradient $\boldsymbol{\nabla }\mathcal{P}^{(1)}$ to determine the transient interface shape. However, the leading-order velocity field (6.21) is challenging to work with directly. To address this, we first express $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\varSigma }_p^{(1)}$ in terms of $\boldsymbol {U}^{(0)}$ , calculate the interface shape, and then transition to numerical computations as a final step. We find (see Appendix B for more detailed calculations)

(6.22) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\varSigma }_p^{(1)} = \beta _p\, {\rm e}^{-\tau }\int _{0}^{\tau }{\nabla} ^2 \boldsymbol {U}^{(1)}\, {\rm e}^{s}\, {\textrm d}s - f(R,\tau )\, \boldsymbol {e}_r, \end{equation}

where

(6.23) \begin{align} f(R,\tau ) &= 2\beta _p\, {\rm e}^{-\tau }\int _{0}^{\tau }R\frac {\partial }{\partial R}\left (\frac {U^{(0)}}{R}\right )\left (\int _{0}^{S}\frac {\partial }{\partial R}\left (\frac {U^{(0)}}{R}\right ){\rm e}^{s}\, {\textrm d}s\right ) {\textrm d}S \notag \\ &\quad {}+\alpha \beta _p\, {\rm e}^{-\tau }\frac {\partial }{\partial R}\left ( \int _{0}^{\tau } {\rm e}^{-S}\, R^2 \left (\int _{0}^{S}\frac {\partial }{\partial R}\left (\frac {U^{(0)}}{R}\right ){\rm e}^{s}\, {\textrm d}s\right )^2 {\textrm d}S\right ). \end{align}

At first order, the momentum equation (6.3) becomes

(6.24) \begin{equation} El^{-1}\left (\frac {\partial \boldsymbol {U}^{(1)}}{\partial \tau } + \boldsymbol {U}^{(0)}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol {U}^{(0)}\right ) = -\boldsymbol{\nabla }\mathcal{P}^{(1)} + \beta _s \,{\nabla} ^2 \boldsymbol {U}^{(1)} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\varSigma }_p^{(1)}. \end{equation}

The second term on the left-hand side of (6.24) represents the centrifugal force generated by the (start-up) vortex:

(6.25) \begin{equation} \boldsymbol {U}^{(0)}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol {U}^{(0)} = -\frac {(U^{(0)})^2}{R}\boldsymbol {e}_r. \end{equation}

Observe that both the centrifugal force (6.25) and the polymeric force induced by the leading-order flow (second term of (6.22)) are directed radially. To simplify, we introduce a modified pressure variable that combines these two conservative forces:

(6.26) \begin{align} {{\rm P}}^{(1)}(R,\tau ) = \mathcal{P}^{(1)}(R,\tau ) + \int _{R}^{\infty }\left ( \frac {(U^{(0)}(s,\tau ))^2}{El\, s} - f(s,\tau )\right ) {\textrm d}s. \end{align}

Combining (6.24) and (6.26) allows us to express the momentum equation as

(6.27) \begin{equation} El^{-1}\,\frac {\partial \boldsymbol {U}^{(1)}}{\partial \tau } = -\boldsymbol{\nabla }{{\rm P}}^{(1)} + {\mathcal{L}}\left ({\nabla} ^2 \boldsymbol {U}^{(1)}\right ). \end{equation}

Following an argument similar to that in §5.2, in the limit where ${G} \gg 1$ , we argue through the dominant balance that the interface shapes are determined by the pressure contributions from the inertial and polymeric stresses in (6.26). We obtain a differential equation

(6.28) \begin{equation} H^{(2)}_{\textit{inertia}}-\frac {1}{\textit{BR}}\,\frac {\partial }{\partial R}\left (R\frac {\partial H^{(2)}_{\textit{inertia}}}{\partial R} \right ) = \frac {1}{{G}} \int _{R}^{\infty }\left (f(s,\tau )-\frac {(U^{(0)}(s,\tau ))^2}{El\, s}\right ) {\textrm d}s := F(R,\tau ). \end{equation}

The solution to this differential equation with the boundary conditions (2.15b ) and (2.15e ) is

(6.29) \begin{align} H^{(2)}_{\textit{inertia}}(R,\tau )&= BI_1\big (\sqrt {B}\big )K_1^{-1}\big (\sqrt {B}\big ) \left (\int _{1}^{\infty }s\,K_0\big (\sqrt {B}s\big )F(s,\tau ) \, {\textrm d}s\right )K_0\big (\sqrt {B}R\big ) \notag \\ &\quad {}+ BI_0\big (\sqrt {B}R\big )\int _{R}^{\infty }s\,K_0\big (\sqrt {B}s\big )F(s,\tau )\, {\textrm d}s \notag \\ &\quad {}+ BK_0\big (\sqrt {B}R\big )\int _{1}^{R}s\,I_0\big (\sqrt {B}s\big )F(s,\tau )\, {\textrm d}s. \end{align}

Before we present numerical results to (6.29), we consider the steady interface shape in the presence of inertial effects. We substitute $U^{(0)} \approx {1}/{R}$ and let $\tau \rightarrow \infty$ to find

(6.30) \begin{equation} F(R,\infty ) = \frac {1}{{G}}\left (-\frac {1}{2ElR^2} + \frac {2\beta _p(1-2\alpha )}{R^4}\right ). \end{equation}

The steady interface shape is readily obtained by substituting (6.30) into (6.29).

At the rod ( $R = 1$ ), the steady rise height can be simplified and is given by

(6.31) \begin{align} H_{\textit{inertia}}^{(2)}(1,\infty ) &= \frac {\sqrt {B}}{{G}}\,K_1^{-1}\big (\sqrt {B}\big )\left (-\frac {1}{2El}\int _{1}^{\infty }R^{-1}K_0\big (\sqrt {B}R\big ) {\textrm d}R \right . \notag \\ &\quad {}+ 2\beta _p(1-2\alpha )\left .\int _{1}^{\infty }R^{-3}\,K_0\big (\sqrt {B}R\big ) {\textrm d}R\right ). \end{align}

Rod climbing, at the steady state, occurs when this rise height is positive, which corresponds to the criterion

(6.32) \begin{equation} \varLambda := 4\beta _p(1-2\alpha )El \gt \frac {\int _{1}^{\infty }R^{-1}\,K_0\big (\sqrt {B}R\big ) {\textrm d}R}{\int _{1}^{\infty }R^{-3}K_0\big (\sqrt {B}R\big )\: {\textrm d}R}. \end{equation}

Figure 5. Phase diagram for the existence of the rod-climbing phenomenon around an infinitely long rotating rod in a Giesekus fluid, incorporating small but finite inertial effects. All calculations were performed in the range $B \in (0.2,25)$ .

Figure 5 illustrates a phase diagram as a function of the parameters $\varLambda$ and $B$ for the rod-climbing criterion (6.32). When $ B = O(1)$ , the right-hand side of (6.32) takes a value approximately between 1 and 2, approaching 1 in the limit as $B \rightarrow \infty$ . A necessary condition, but not sufficient, for the rod climbing at the steady state to occur, regardless of surface tension effects, is that $\varLambda \gt 1.$ For a given Bond number $B$ , the climbing criterion (6.32) requires viscoelastic effects to surpass inertial effects for climbing to occur (represented by the elasticity number $El$ ), though the effective viscoelasticity is also reduced by the polymer solution’s diluteness (represented by $\beta _p$ ) and the polymer’s chain mobility (represented by $\alpha$ ). This highlights how the balance between viscoelasticity, inertia, polymer concentration and polymer mobility governs the climbing behaviour.

Figure 6. Steady-state interface shapes around an infinitely long rotating rod in a Giesekus fluid, incorporating small but finite inertial effects, for ( $a$ ) $\varLambda =1$ , ( $b$ ) $\varLambda =1.5$ , ( $c$ ) $\varLambda =2$ , and Bond numbers $B = 1,5, 25$ .

Figure 6 shows the steady-state interface shape $ H_{\textit{inertia}}(R,\infty )$ scaled by $2\beta _p(1 - 2\alpha )\, \textit{Wi}/\mathcal{G}$ , for various values of $\varLambda$ and $B$ . We consider three representative cases: $\varLambda = 1, 1.5, 2$ , which illustrate typical steady interface behaviours. For $\varLambda = 1$ , the interface height is negative when $B \geqslant 1$ ; for $ \varLambda = 2$ , it is positive in the same range. For $\varLambda = 1.5$ , the interface height may be either positive or negative depending on the Bond number: rod climbing occurs for $B \gtrsim 5$ , but not for smaller values.

Figure 7. Time evolution of the transient interface shapes around an infinitely long rotating rod in a Giesekus fluid, incorporating small but finite inertial effects, shown at times $\tau = 0.1, 1, 10, 25$ , for ( $a$ ) $El=0.5$ , ( $b$ ) $El=1$ , and ( $c$ ) $El=2$ . All calculations were performed using $\beta _p = 0.5$ , $\alpha = 0.25$ and $B = 9$ .

To examine the interface dynamics before reaching the steady state, we present the transient interface (6.29) obtained via numerical computations. In figures 7 and 8, we show the time evolution of the free surface (scaled by $2\beta _p(1-2\alpha )\,\textit{Wi}/\mathcal{G}$ ) for different elasticity numbers and Bond numbers. We observe that the impulsive motion of the rod initially causes a sudden dip in the fluid interface near the rod, an effect that is more pronounced at lower elasticity numbers due to stronger inertial forces. As time progresses, inertia continues to drive the interface downwards, while viscoelastic stresses act to push it upwards. When inertial effects outweigh viscoelastic effects, the interface height decreases monotonically until reaching a steady state, as illustrated for $El = 0.5$ in figure 7. At higher elasticity numbers, where viscoelastic effects are more dominant, the interface initially climbs due to the stronger viscoelastic response. However, because viscoelastic stresses relax faster than inertial effects, the climbing height eventually peaks and then decreases, governed by the slower relaxation of inertia towards its steady-state behaviour.

Figure 8. Time evolution of the transient interface shapes around an infinitely long rotating rod in a Giesekus fluid, incorporating small but finite inertial effects, shown at times $\tau = 0.1, 1, 10, 25$ , for ( $a$ ) $B=1$ , ( $b$ ) $B=5$ , and ( $c$ ) $B=25$ . All calculations were performed using $\beta _p = 0.5$ , $\alpha = 0.25$ and $ El = 2$ .

Depending on the spatial location of the interface, the viscoelastic fluid rises at varying rates due to the interplay between inertial and viscoelastic contributions. For $El = 0.5, 1, 2$ , we observe that the interface shape exhibits non-monotonic behaviour primarily within the range $ 1 \leqslant R \leqslant 1.5$ . At larger values of $R$ , the interface deformation becomes smaller and more regular.