1. Introduction
Viscoelastic shear flows exhibit distinct features compared with Newtonian flows, primarily due to the additional viscoelastic stresses arising from the stretching and relaxation of polymers, even in dilute solutions (Sanchez et al. Reference Sanchez, Jovanovic, Kumar, Morozov, Shankar, Subramanian and Wilson2022). In addition to the Newtonian turbulence that is dominated by large inertial forces, two viscoelastic turbulent states emerge: elastic turbulence (ET), which occurs at vanishingly small inertia (Groisman & Steinberg Reference Groisman and Steinberg2000; Steinberg Reference Steinberg2021), and elasto-inertial turbulence (EIT), for which the effects of inertia and elasticity are comparable (Samanta et al. Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013; Dubief, Terrapon & Hof Reference Dubief, Terrapon and Hof2023). The origins of these two turbulent states are of particular interest, and studies aimed at identifying their transition routes are still ongoing. For a comprehensive review of this topic, one can refer to Couchman et al. (Reference Couchman, Beneitez, Page and Kerswell2024).
In viscoelastic flows, experimental techniques do not allow for direct measurements of the polymer conformation tensor. Numerical simulations can compensate for this limitation, making them an efficient tool for investigating the transition routes to ET or EIT. However, as highlighted by Alves, Oliveira & Pinho (Reference Alves, Oliveira and Pinho2021), the lack of diffusion in the conformation tensor equations poses significant numerical difficulties in simulating viscoelastic fluid flows, especially at high Weissenberg numbers. To achieve numerical stability, adding an artificial conformation diffusion (ACD) term to the conformation equation appears as an efficient approach. This term is expressed as
$({\textit{Re}}Sc)^{-1}{\nabla} ^2 \boldsymbol {c}$
, where
${\textit{Re}}$
denotes the Reynolds number,
$Sc$
is the Schmidt number (defined by the ratio of the solvent viscosity to diffusivity of polymers in solution) and
$\boldsymbol {c}$
represents the polymer conformation tensor. Although there is always a minimal degree of molecular diffusion of the polymer stress in reality, quantified by
$Sc^{-1}\sim 10^{-6}$
, this level is too small to stabilise the computations (Dubief et al. Reference Dubief, Terrapon and Hof2023). Using a spectral method, Sureshkumar & Beris (Reference Sureshkumar and Beris1995) demonstrated that numerical stability requires a diffusion coefficient of
$({\textit{Re}}Sc)^{-1}\sim 10^{-3}$
, or in their case,
$Sc\approx 0.2$
. This study provided justification for the subsequent use of
$Sc \lt 1$
, although other simulations have employed different
$Sc$
values (Dubief et al. Reference Dubief, Terrapon and Hof2023).
Although the ACD has faced criticism, it remains a widely employed approach in polymer-flow simulations due to its simplicity and robustness. However, the recent identification of a new type of instability induced by ACD, referred to as polymer diffusive instability (PDI), complicates this situation. In viscoelastic plane Couette flows, Beneitez, Page & Kerswell (Reference Beneitez, Page and Kerswell2023) first reported the emergence of an unstable wall-mode PDI that arises at vanishing inertia and low diffusivity of the polymer stress. This instability was found to occur even at extremely small, non-zero diffusion levels, and was claimed to be robust against various boundary conditions applied to the polymer conformation equation. Direct numerical simulations indicated that PDI could trigger a self-sustaining chaotic state, offering a potential mechanism for the origin of ET in numerical simulations with ACD. Subsequently, Couchman et al. (Reference Couchman, Beneitez, Page and Kerswell2024) conducted a systematic study of PDI across a wide parameter space, employing both Oldroyd-B and the finitely extensible nonlinear elastic Peterlin (FENE-P) constitutive models for plane Couette flows and plane Poiseuille flows. Their investigations also considered configurations with finite Reynolds numbers, suggesting that PDI may also lead to EIT in addition to ET in polymer fluids. Furthermore, through two-dimensional numerical simulations, Beneitez et al. (Reference Beneitez, Page, Dubief and Kerswell2024) established a transition route to ET and EIT based on the secondary instability of PDI. In recent studies, Lewy & Kerswell (Reference Lewy and Kerswell2024) particularly focused on the PDI in highly concentrated polymeric solutions, and Beneitez, Mrim & Kerswell (Reference Beneitez, Mrim and Kerswell2025) confirmed the emergence of PDI even when explicit ACD is not present. Through asymptotic analysis by taking both the diffusivity and the solvent-to-total-viscosity ratio to be small, the instability system was significantly simplified, reducing the tenth-order ordinary differential equation system to a system of two algebraic equations. This simplification enables an in-depth exploration of the PDI mechanism.
Although the nature of PDI – whether it is physical or artificial – remains a topic of debate, its discovery poses challenges to previous transition simulations that employed ACD. The PDI is fundamentally different from numerical instability, as it persists even when the artificial diffusivity is very small. One may expect that suppressing or avoiding the unstable PDI is beneficial in simulations of polymer fluids incorporating ACD. It is important to emphasise that the inclusion of polymer diffusion increases the order of the differential system, requiring additional boundary conditions for the polymer conformation at the wall. Considering three representative sets of boundary conditions (the detailed forms of which are presented in § 2.2), Beneitez et al. (Reference Beneitez, Page and Kerswell2023) demonstrated the presence of unstable zones of PDI for all examined boundary conditions. However, variations in the neutral curves were observed depending on the choice of the conformation boundary conditions. Lewy & Kerswell (Reference Lewy and Kerswell2024) also reported that replacing the conformation boundary conditions from the vanishing diffusion to the Neumann conditions leads to a suppression of PDI in the high-concentration limit. These observations raise a pertinent question: Is there a boundary condition under which the unstable zone of PDI could be eliminated? If such a boundary condition exists, it could potentially offer a valid approach for simulating the transition process of polymer fluids with the introduction of ACD. Addressing this question is inherently challenging, as it necessitates exploring the elimination of PDI across the entire parameter space. Nevertheless, since the generic instability mechanism can be captured through its asymptotic structure (as has been done in Lewy & Kerswell Reference Lewy and Kerswell2024), we can develop the required conformation boundary conditions inspired by the asymptotic analysis. In particular, the perturbation behaviours observed in the near-wall region can offer valuable insights for selecting the boundary conditions. This is the primary focus of this paper.
The rest of this paper is organised as follows. In § 2, we present the mathematical framework of the linear viscoelastic instability in parallel wall-shear flows based on the Oldroyd-B constitutive model, along with conformation diffusion. This is followed by an exploration of its asymptotic structure in § 3, where we identify four regimes with reducing complexity from regime I to regime IV. Based on the simplified equation system in regime IV, we propose a set of new boundary conditions to eliminate PDI in § 4.1, and assess its applicability in the other asymptotic regimes in § 4.2 and the original system in § 4.3. In § 4.4, we further validate the applicability of the new boundary conditions within the FENE-P constitutive model. Then, the new boundary conditions are integrated in direct numerical simulation (DNS) code, showing the effectiveness in eliminating PDI for a representative configuration in § 4.5. Finally, we conclude with remarks in § 5.
2. Physical model and mathematical descriptions
2.1. Formulation
Following Couchman et al. (Reference Couchman, Beneitez, Page and Kerswell2024), we consider two types of parallel wall-shear flows of polymer solutions: a plane Poiseuille flow (PPF) driven by a pressure gradient and a plane Couette flow (PCF) driven by two counter-moving parallel plates. By selecting the half-channel width
$H$
as the unit length and the maximum velocity
$U_{\textit{max}}$
(the centreline velocity for PPF and the wall velocity for PCF) as the reference velocity, we define two dimensionless parameters, the Reynolds number
${\textit{Re}}$
and the Weissenberg number
$\textit{Wi}$
\begin{equation} {\textit{Re}}=\frac {\rho U_{\textit{max}}H}{\mu },\quad \textit{Wi}=\frac {\bar \lambda U_{\textit{max}}}{H}, \end{equation}
where
$\rho$
denotes the density, and
$\mu$
and
$\bar \lambda$
denote the total viscosity and the relaxation time of the polymer molecules (to their equilibrium states), respectively. Here, the total viscosity is the sum of the solvent viscosity
$\mu _s$
and the polymer viscosity
$\mu _p$
(Bird et al. Reference Bird, Curtiss, Armstrong and Hassager1987), i.e.
$\mu =\mu _s+\mu _p$
, and a viscosity ratio
$\beta$
is defined as
\begin{equation} \beta =\mu _s/\mu \in [0,1]. \end{equation}
The dimensionless Cartesian coordinate system
$(x,y,z)$
is employed, with
$x$
,
$y$
and
$z$
denoting the streamwise, wall-normal and spanwise directions, respectively. The velocity field
$\boldsymbol u=(u,v,w)$
is normalised by
$U_{\textit{max}}$
, the time
$t$
is normalised by
$H/U_{\textit{max}}$
and the pressure
$p$
is normalised by
$\rho U_{\textit{max}}^2/{\textit{Re}}$
. The conformation tensor
$\boldsymbol c$
and stress
$\boldsymbol \tau$
are normalised by
$k_BT/H_0$
and
$\mu _p U_{\textit{max}}/H$
, respectively, where
$k_B$
is the Boltzmann constant,
$T$
is temperature and
$H_0$
is the spring constant in the elastic dumbbell model of the polymer.
The dimensionless equations governing the viscoelastic fluid with ACD read (Couchman et al. Reference Couchman, Beneitez, Page and Kerswell2024)
\begin{align}& {\textit{Re}}\left (\frac {\partial \boldsymbol u}{\partial t}+(\boldsymbol u\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol u\right )=-\boldsymbol{\nabla }p+\beta {\nabla} ^2\boldsymbol u+(1-\beta )\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol \tau ,\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol u=0, \end{align}
\begin{align}&\qquad\qquad \frac {\partial \boldsymbol c}{\partial t}+(\boldsymbol u\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol c+\tau -\boldsymbol c\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol u-(\boldsymbol{\nabla }\boldsymbol u)^T\boldsymbol{\cdot }\boldsymbol c=\epsilon {\nabla} ^2\boldsymbol c, \end{align}
where
$\epsilon$
denotes the artificial diffusivity introduced in the polymer stress equation. It is related to the Schmidt number through
$Sc=(\epsilon {\textit{Re}})^{-1}$
. Notably, when
${\textit{Re}}=0$
, one can avoid introducing the Schmidt number and directly use
$\epsilon$
to denote the conformation diffusivity. For representative demonstration, the following asymptotic analysis will be based on the Oldroyd-B constitutive model, for which
$\boldsymbol \tau$
and
$\boldsymbol c$
are related by
\begin{equation} \boldsymbol \tau =\frac {\boldsymbol c-\boldsymbol I}{\textit{Wi}}, \end{equation}
with
$\boldsymbol I$
denoting the unity matrix. We would like to emphasise that the boundary conditions to be developed in this paper will also be applied to the widely used FENE-P constitutive model, because the asymptotic structures for both constitutive models are the same. This will be confirmed in § 4.4.
2.2. Linear stability theory (LST)
To analyse the linear instability, we decompose the flow field
$\boldsymbol \phi =(\boldsymbol u,p,\boldsymbol c)$
into a steady base flow
$\boldsymbol \varPhi =(\boldsymbol U,P,\boldsymbol C)$
and a harmonic perturbation
$\hat {\boldsymbol \phi }=(\hat {\boldsymbol u},\hat p,\hat {\boldsymbol c})$
\begin{equation} \boldsymbol \phi =\boldsymbol \varPhi (y)+\hat \epsilon \hat {\boldsymbol \phi }(y){\textrm {e}}^{{\textrm {i}}(kx+\gamma z-\omega t)}+ \textrm{c.c.}, \end{equation}
where
$\hat \epsilon \ll 1$
denotes the perturbation amplitude,
$k$
and
$\gamma$
denote the streamwise and spanwise wavenumbers,
$\omega$
is the complex frequency with its imaginary part denoting the growth rate and
$\textrm{c.c.}$
denotes the complex conjugate. In this paper, we focus on two-dimensional perturbations, for which
$\gamma =0$
. For parallel wall-shear flows,
$\boldsymbol U=(U(y),0,0)$
with
$y\in [-1,1]$
, and the base-flow profile is
\begin{equation} U=\left \{\!\!\begin{array}{ll} 1-y^2& \mbox{for a PPF}, \\ y&\mbox{for a PCF}. \end{array} \right . \end{equation}
For convenience, we introduce the wall shear
$\lambda :=U_y(0)$
, with
$\lambda =2$
and 1 for PPF and PCF, respectively. The non-zero elements of
$\boldsymbol C$
for the Oldroyd-B constitutive model are
\begin{equation} C_{11}=1+2\textit{Wi}^2U'^2,\quad C_{12}=C_{21}=\textit{WiU}',\quad C_{22}=C_{33}=1, \end{equation}
where, throughout this paper, a prime denotes the derivative with respect to its argument.
Substituting (2.5) into the governing equation system (2.3) and (2.4), and retaining the
$O(\hat \epsilon )$
terms, we obtain the following linear system:
\begin{align}&\quad\,\, {\textrm {i}} k {\textit{Re}}(U-c)\hat u+{\textit{ReU}}'\hat v+{\textrm {i}} k\hat p={\beta }{}\big(\hat u^{\prime\prime}-k^2\hat u\big)+{(1-\beta)}\left({\textrm {i}} k\hat \tau _{11}+\hat \tau _{12}^{\prime}\right)\!, \end{align}
\begin{align}&\, {\textrm {i}} k {\textit{Re}}(U-c)\hat v+\hat p'={\beta }{}\big(\hat v^{\prime\prime}-k^2\hat v\big)+{(1-\beta)}{}\left({\textrm {i}} k\hat \tau _{12}+\hat \tau _{22}^{\prime}\right)\!,\quad {\textrm {i}} k \hat u+\hat v'=0, \end{align}
\begin{align}& \big[{\textrm {i}} k (U-c)+\epsilon k^2\big]\hat c_{11}+\hat \tau _{11}+C_{11}^{\prime}\hat v-2\left({\textrm {i}} k C_{11}\hat u+C_{12}\hat u'+U'\hat c_{12}\right)=\epsilon \hat c_{11}^{\prime\prime}, \end{align}
\begin{align}& \big[{\textrm {i}} k (U-c)+\epsilon k^2\big]\hat c_{12}+\hat \tau _{12}+C_{12}^{\prime}\hat v-\left({\textrm {i}} k C_{12}\hat v+C_{22}\hat u'+U'\hat c_{22}\right)=\epsilon \hat c_{12}^{\prime\prime}, \end{align}
\begin{align}&\qquad\quad\,\, \big[{\textrm {i}} k (U-c)+\epsilon k^2\big]\hat c_{22}+\hat \tau _{22}-2\left({\textrm {i}} k C_{12}\hat v+C_{22}\hat v'\right)=\epsilon \hat c_{22}^{\prime\prime}, \end{align}
\begin{align}&\qquad\qquad\qquad\qquad \big[{\textrm {i}} k (U-c)+\epsilon k^2\big]\hat c_{33}+\hat \tau _{33}=\epsilon \hat c_{33}^{\prime\prime}, \end{align}
where
$c=\omega /k=c_r+{\textrm {i}} c_i$
and the Oldroyd-B constitutive model determines
\begin{align} (\hat \tau _{11},\hat \tau _{12},\hat \tau _{22},\hat \tau _{33})=\textit{Wi}^{-1}(\hat c_{11},\hat c_{12},\hat c_{22},\hat c_{33}). \end{align}
Note that, because
$\hat c_{33}$
is decoupled from the Oldroyd-B instability system, we do not need to consider its impact in the following, and directly introduce the boundary condition
$\hat c_{33}(\pm 1)=0$
.
At both walls
$y=\pm 1$
, the no-slip, non-penetration conditions are applied
\begin{align} \hat u(\pm 1)=\hat v(\pm 1)=0.\quad \end{align}
Following Beneitez et al. (Reference Beneitez, Page and Kerswell2023), we consider three representative types of boundary conditions for the conformation tensor. For the type-(i) boundary conditions (referred to as BC1), the diffusive effect is neglected at the wall, leading to
\begin{align} \big[\big({-{\textrm {i}} k c+\textit{Wi}^{-1}}\big)\hat c_{11}-{2\left(C_{12}\hat u'+U'\hat c_{12}\right)}\big]_{y\,=\,\pm 1}=0, \end{align}
\begin{align} \big[\big({-{\textrm {i}} k c+\textit{Wi}^{-1}}\big)\hat c_{12}-{\hat u'}\big]_{y\,=\,\pm 1}=0,\quad \hat c_{22}(\pm 1)=0. \end{align}
For the type-(ii) boundary conditions (referred to as BC2), only the terms
$\epsilon \partial _{yy}{\boldsymbol C}_{ij}$
are set to be zero at the wall
\begin{align} \big[\big({-{\textrm {i}} k c+\textit{Wi}^{-1}}+k^2\big)\hat c_{11}-{2\left(C_{12}\hat u'+U'\hat c_{12}\right)}\big]_{y\,=\,\pm 1}=0, \end{align}
\begin{align} \big[\big({-{\textrm {i}} k c+\textit{Wi}^{-1}}+k^2\big)\hat c_{12}-{\hat u'}\big]_{y\,=\,\pm 1}=0,\quad \hat c_{22}(\pm 1)=0. \end{align}
For the type-(iii) boundary conditions (referred to as BC3), the conformation tensor is assumed to show zero gradient at the wall (
$\partial _y{\boldsymbol C}_{ij}|_{y\,=\,\pm 1}=0$
)
\begin{align} \hat c_{11}^{\prime}(\pm 1)=\hat c_{12}^{\prime}(\pm 1)=\hat c_{22}^{\prime}(\pm 1)=0. \end{align}
The equation system (2.8) and (2.9), together with the boundary conditions (2.10) and one of (2.11) to (2.13), constitutes an eigenvalue problem, where
$c$
serves as the eigenvalue. The controlling parameters for this system include
$\textit{Wi}$
,
${\textit{Re}}$
,
$k$
,
$\beta$
and
$\epsilon$
. In this paper, the eigenvalue system is solved using the following two numerical approaches:
-
(i) spectral collocation (SC) method: this method is based on Chebyshev polynomials and is a global method that excels in calculating the eigenvalues across a wide spectrum. A detailed description can be found in Zhang (Reference Zhang2021);
-
(ii) Newton iterative (NI) method: the equations are discretised using the fourth-order compact Malik scheme (Malik Reference Malik1990), and the eigenvalue solution is obtained by NI method to minimise the determinant of the coefficient matrix to a sufficiently small value (e.g.
$10^{-9}$
). This local method focuses on a single eigenvalue based on an initial guess and is well known for its high accuracy. Our previous eigenvalue calculations in viscoelastic pipe flows (Dong & Zhang Reference Dong and Zhang2022) and Newtonian boundary-layer flows (Ji, Dong & Zhao Reference Ji, Dong and Zhao2023; Song & Dong Reference Song and Dong2023) employed this method.
In this paper, unless otherwise specified, most results are obtained using the NI method, with initial values for iteration derived either from the SC method or from a continuation of existing eigenvalue solutions.
3. Asymptotic analyses of PDI
3.1. Guideline of the asymptotic analyses
As a guideline for the subsequent asymptotic analyses, this subsection outlines the scaling relations for the various asymptotic regimes along with their corresponding structures, while detailed mathematical derivations are provided in the following sections. Based on the linear system (2.8), Couchman et al. (Reference Couchman, Beneitez, Page and Kerswell2024) computed the neutral curves for PDI in both PPF and PCF, revealing the emergence of the thin diffusive layer in the
$O(\sqrt {\epsilon })$
vicinity of the wall. In this layer, the conformation diffusivity becomes the dominant influence. At the critical state where PDI emerges, a clear scaling relation of
$k\sim \epsilon ^{-1/2}$
is observed, with the critical Weissenberg number reaching several tens. Notably, the lower-branch neutral curve of PDI manifests even when
$\textit{Wi}$
approaches infinity, with the wavenumber decreases accordingly. In light of these findings, we conduct an asymptotic analysis under the assumption that
$\epsilon \ll 1$
while taking
$\textit{Wi}=O(1)$
, referred to as regime I. As will be shown in § 3.2.2, the most unstable PDI in regime I appears in the limit
$\textit{Wi}\to \infty$
and
$k\to 0$
, with the scaling
$\bar k=k\epsilon ^{1/2}\sim \textit{Wi}^{-3/2}$
. We then derive the instability properties in these limits, referred to as regime II. Furthermore, to simplify the asymptotic system, we consider the dilute limit – similar to the centre-mode analysis in Dong & Zhang (Reference Dong and Zhang2022). In this limit, the rescaled wavenumber,
$\tilde k=k\epsilon ^{1/2}\textit{Wi}^{3/2}$
, scales as
$(1-\beta )^{-3/2}$
, so that the resulting asymptotic system is characterised by a single controlling parameter,
$\breve \sigma =\tilde k^{2/3}(1-\beta )$
. This regime is referred to as regime III. Finally, we show that, if the limit
$\breve \sigma \to \infty$
is considered, unstable PDI is still maintained (termed as regime IV); however, this asymptotic state exhibits a parameter-free character, which facilitates the derivation of analytical solutions. In summary, this paper considers the following four regimes:
-
(i) regime I:
$\textit{Wi}=O(1)$
,
$1-\beta =O(1)$
,
$k\sim \epsilon ^{-1/2}$
and
$c-U(-1)\sim \epsilon ^{1/2}$
; -
(ii) regime II:
$\textit{Wi}\gg 1$
,
$1-\beta =O(1)$
,
$k\epsilon ^{1/2}\sim \textit{Wi}^{-3/2}\ll 1$
and
$c-U(-1)\sim \epsilon ^{1/2}\textit{Wi}^{1/2}$
; -
(iii) regime III:
$\textit{Wi}\gg 1$
,
$\beta \to 1$
,
$k\epsilon ^{1/2}\textit{Wi}^{3/2}\sim (1-\beta )^{-3/2}\gg 1$
and
$c-U(-1)\sim \epsilon ^{1/2}\textit{Wi}^{1/2}\tilde k^\alpha$
, where
$\alpha =-1$
for BC1/BC2 and
$=-1/3$
for BC3; -
(iv) regime IV:
$\textit{Wi}\gg 1$
,
$\beta \to 1$
,
$k^{2/3}\epsilon ^{1/3}Wi(1-\beta )\gg 1$
and
$c-U(-1)\sim \epsilon ^{1/2}\textit{Wi}^{1/2}\tilde k^{-1}\breve \sigma ^{1/2}$
for BC1/BC2 and
$\sim \epsilon ^{1/2}\textit{Wi}^{1/2}\tilde k^{-1/3}\breve \sigma ^{-1/2}$
for BC3.
The asymptotic structures of these regimes are depicted in figure 1. In regime I, the instability primarily manifests within the diffusive layer (layer (ii)) in the
$O(\sqrt {\epsilon })$
vicinity of the wall, while in regimes II and III, the diffusive layer splits into two distinct layers: an upper diffusive layer (UDL, layer (ii-1)) and a lower diffusive layer (LDL, layer (ii-2)) (with the mathematical details illustrated in § 3.3). Remarkably, the thicknesses of the UDL and LDL are of different magnitudes for regimes II and III. For the former, the thicknesses of the UDL and LDL are
$O(\textit{Wi}^{3/2}\epsilon ^{-1/2})$
and
$O(\textit{Wi}^{1/2}\epsilon ^{-1/2})$
, respectively. In contrast, for the low-concentration regime III, the rescaled wavenumber
$\tilde k=k\epsilon ^{1/2}\textit{Wi}^{3/2}$
must be large, resulting in the modified thicknesses of
$O(\tilde k^{-1}\textit{Wi}^{3/2}\epsilon ^{-1/2})$
and
$O(\tilde k^{-1/3}\textit{Wi}^{1/2}\epsilon ^{-1/2})$
, respectively. Specifically, for regime III, the magnitude of the phase speed correction,
$c-U(-1)$
, may vary depending on the specific conformation boundary conditions imposed. Regime IV is actually the asymptotic state of regime III by taking
$\breve \sigma$
to be asymptotically large. In this limit, the LDL further splits into two sublayers: a bulk sublayer (layer (iii)) and a wall sublayer (layer (iv)). It is worth noting that regime I converges to regime II as
$\textit{Wi}$
becomes large, regime II converges to regime III as
$\beta$
approaches unity and regime III converges to regime IV as
$\breve \sigma$
becomes large.
A summary of the reduction process from the instability system (2.8) to regime IV. Panels (a), (b) and (c) illustrate the asymptotic structures for regimes I, II/III and IV, respectively. Layer (i) denotes the main layer; layer (ii) denotes the diffusive layer; layers (ii-1) and (ii-2) denote the UDL and LDL, respectively; layers (iii) and (iv) denote the bulk sublayer and the wall sublayer, respectively. The elongated chains represent the polymer states in the high and low shear rates.
3.2. Regime I analysis
3.2.1. Formulation for regime I
This subsection focuses on
$\textit{Wi}=O(1)$
and
$k\sim \epsilon ^{-1/2}$
, and thus we introduce
\begin{equation} \bar k=\epsilon ^{1/2}k=O(1). \end{equation}
Without loss of generality, we only consider the perturbations near the lower wall
$y=-1$
. By balancing the unsteady terms with the convective terms in (2.8a
,
b
), we can estimate the magnitude of the phase speed, which can be expressed as
\begin{equation} c=U(-1)+\epsilon ^{1/2}c_1+\cdots , \end{equation}
where
$U(-1)=0$
for a PPF and
$-1$
for a PCF.
In this regime, the perturbations are mainly concentrated in the near-wall diffusive layer, as depicted in figure 1(a). In this layer, we introduce a local coordinate
\begin{equation} Y=\epsilon ^{-1/2}(y+1)=O(1), \end{equation}
and the velocity and conformation tensor of the base flow are approximated as
\begin{equation} U=\sqrt {\epsilon }\lambda Y+O(\epsilon ),\quad (C_{11},C_{12},C_{22},C_{33})=\big(2\textit{Wi}^2\lambda ^2+1, \textit{Wi}\lambda ,1,1\big)+O(\sqrt {\epsilon }). \end{equation}
From a scaling estimate from the equation system (2.8) and (2.9), we can expand the perturbation field as
\begin{align} \left (\hat u,\hat v,\hat p\right )&=\big(\hat u_0,\hat v_0,\epsilon ^{-1/2}{\beta }{}\hat p_0\big)+\cdots ,\nonumber\\ \left (\hat c_{11},\hat c_{12},\hat c_{22}\right )&=\epsilon ^{-1/2}\textit{Wi}^2\left (\hat c_{110},\hat c_{120},\hat c_{220}\right )+\cdots . \end{align}
Substituting the above expansion into the instability system (2.8) and (2.9), and collecting the leading-order terms, we obtain
\begin{align}&\qquad\quad \hat u_0^{\prime\prime}-\bar k^2\hat u_0+\sigma \textit{Wi}\left({\textrm {i}}\bar k\hat c_{110}+\hat c_{120}^{\prime}\right)-{\textrm {i}}\bar k\hat p_0=0, \end{align}
\begin{align}& \hat v_0^{\prime\prime}-\bar k^2\hat v_0+ \sigma \textit{Wi}\left({\textrm {i}}\bar k\hat c_{120}+\hat c_{220}^{\prime}\right)-\hat p_0^{\prime}=0,\quad {\textrm {i}}\bar k\hat u_0+\hat v_0^{\prime}=0, \end{align}
\begin{align}&\, \hat c_{110}^{\prime\prime}= S\hat c_{110}-2{\textrm {i}} \bar k\big(2\lambda ^2+\textit{Wi}^{-2}\big)\hat u_0-2\textit{Wi}^{-1}\lambda \hat u_0^{\prime}-2\lambda \hat c_{120}, \end{align}
\begin{align}&\quad\,\,\, \hat c_{120}^{\prime\prime}= S\hat c_{120}-{\textrm {i}} \bar k\big(2\lambda ^2+\textit{Wi}^{-2}\big)\hat v_0-\textit{Wi}^{-2}\hat u_0^{\prime}-\lambda \hat c_{220}, \end{align}
\begin{align}&\qquad\qquad \hat c_{220}^{\prime\prime}= S\hat c_{220}-2\big({\textrm {i}}\bar k \lambda \textit{Wi}^{-1}\hat v_0+\textit{Wi}^{-2}\hat v_0^{\prime}\big), \end{align}
where
$S={\textrm {i}}\bar k(\lambda Y-c_1)+\textit{Wi}^{-1}+\bar k^2$
, and
$\sigma =(1-\beta )/\beta$
. Here, the equation for
$\hat c_{33}$
is omitted, as it does not influence the other perturbation quantities. We assume that
$\epsilon {\textit{Re}} \ll 1$
, which allows us to neglect the inertia terms in the leading-order momentum equation. Considering that the perturbations should damp to match the main layer in the upper limit, we impose the attenuation conditions for the perturbation
\begin{equation} (\hat u_0,\hat v_0,\hat p_0,\hat c_{110},\hat c_{120},\hat c_{220})\to 0\quad \mbox{as }Y\to \infty . \end{equation}
Following the boundary conditions from equations (2.10) to (2.13), we derive the wall boundary conditions for the asymptotic system presented in (3.5): the BC1 conditions read
\begin{align}& \hat u_0(0)=\hat v_0(0)=0,\quad \big(-{\textrm {i}} \bar kc_1+\textit{Wi}^{-1}\big)\hat c_{110}(0)=2\lambda \left (\frac {\hat u_0^{\prime}(0)}{Wi }+{ \hat c_{120}(0)}\right )\!, \end{align}
\begin{align}&\qquad\qquad \big(-{\textrm {i}} \bar kc_1+\textit{Wi}^{-1}\big)\hat c_{120}(0)=\frac {\hat u_0^{\prime}(0)}{\textit{Wi}^2},\quad \hat c_{220}(0)=0; \end{align}
the BC2 conditions read
\begin{align}& \hat u_0(0)=\hat v_0(0)=0,\quad \big(-{\textrm {i}} \bar kc_1+\textit{Wi}^{-1}+\tilde k^2\big)\hat c_{110}(0)=2\lambda \left (\frac {\hat u_0^{\prime}(0)}{Wi }+{ \hat c_{120}(0)}\right )\!, \end{align}
\begin{align}&\qquad\qquad \big(-{\textrm {i}} \bar kc_1+\textit{Wi}^{-1}+\tilde k^2\big)\hat c_{120}(0)=\frac {\hat u_0^{\prime}(0)}{\textit{Wi}^2},\quad \hat c_{220}(0)=0; \end{align}
the BC3 conditions read
\begin{align} \hat u_0(0)=\hat v_0(0)=\hat c_{11}^{\prime}(0)=\hat c_{12}^{\prime}(0)=\hat c_{22}^{\prime}(0)=0. \end{align}
In this regime, the eigenvalue system discussed in § 2.2 is recast to (3.5) with boundary conditions (3.6) and one of (3.7) to (3.9), with
$c_1$
being the eigenvalue. In this system, only three independent controlling parameters (
$\textit{Wi}$
,
$\beta$
and
$\bar k$
) are present, indicating a simplification from the original system. Such an asymptotic-analysis-inspired simplification was also achieved in Lewy & Kerswell (Reference Lewy and Kerswell2024). This system can be solved using the local NI numerical approach, with the initial value derived from solving (2.8) and (2.9) using the global SC method. In the NI calculations, we select a computational domain of
$Y\in [0,1000]$
, and discretise the governing equations into 2001 grid points, with a clustering of points in the near-wall region. The results are confirmed to be of sufficient accuracy by comparing with those obtained using a doubled computational domain and twice the number of grid points.
3.2.2. Numerical results for regime I
By selecting a viscoelastic PPF with
$\lambda =2$
and
$\beta =0.9$
, we calculate the rescaled growth rate
$c_{1i}$
of PDI for the three wall boundary conditions represented in (3.7) to (3.9). The resulting contours in the
$\textit{Wi}$
-
$\bar k$
plane are depicted in figure 2. Although unstable zones emerge at moderate and high
$\textit{Wi}$
values for all considered boundary conditions, the overall shapes of these unstable zones vary significantly. The largest unstable zone is observed in panel (a), corresponding to BC1. For each
$\textit{Wi}$
value exceeding the critical threshold
$\textit{Wi}_c$
, two points with zero growth rates appear, referred to as the lower- and upper-branch neutral points. As
$\textit{Wi}$
increases, the most unstable state, indicated by the peak growth rate, occurs at a very low
$\bar k$
value, in close proximity to the lower-branch neutral point. This reveals a regime of particular interest in the limit of
$\textit{Wi}\gg 1$
and
$\bar k\ll 1$
.
Contours of the rescaled growth rate
$c_{1i}$
of the regime-I PDI for various wall boundary conditions for PPF with
$\beta =0.9$
, where the circles mark the critical point for the instability onset (
$\textit{Wi}_c$
,
$\bar k_c$
): (a) BC1; (b) BC2; (c) BC3.
In figure 3, we choose BC1 and plot the neutral curves for representative
$\beta$
values for both PPF and PCF. Following Couchman et al. (Reference Couchman, Beneitez, Page and Kerswell2024), the horizontal axis is rescaled by
$\lambda \textit{Wi}$
. For a given
$\beta$
value, the neutral curves for PPF and PCF are quite similar, with the
$\bar k$
band for PCF appearing in a slightly lower range. Reducing the polymer concentration (i.e. increasing
$\beta$
) leads to a shift of the unstable zone to higher
$\lambda \textit{Wi}$
values, while the wavenumber band of
$\bar k$
remains almost unchanged. The critical values of
$\lambda \textit{Wi}_c$
and
$\bar k_c$
are summarised in table 1. For each
$\beta$
value, the critical values of
$\lambda \textit{Wi}_c$
are almost identical for PCF and PPF cases, aligning with the observation in Couchman et al. (Reference Couchman, Beneitez, Page and Kerswell2024). However, the critical values of
$\bar k_c$
exhibit slight differences. Remarkably, for the lower-branch neutral curve, the rescaled wavenumber
$\bar k$
scales as
$\textit{Wi}^{-3/2}$
. By further fitting the coefficient of the scaling relation
${\mathcal C}=\bar k/\textit{Wi}^{-3/2}$
, we obtain the coefficients for each
$\beta$
value, as summarised in table 2. These coefficients increase remarkably as
$\beta$
approaches unity.
Neutral curves of regime-I PDI in the
$\bar k$
-
$\lambda \textit{Wi}$
plane for BC1 with representative
$\beta$
values, where lines with and without are for the PCF (
$\lambda =1$
) and PPF (
$\lambda =2$
) flows, respectively. The fitting curve
$\bar k\sim \textit{Wi}^{-3/2}$
indicates the scaling law of the lower-branch neutral curve.
Critical values for various
$\beta$
and
$\lambda$
values.
Fitting of the scaling coefficient
${\mathcal C}=\bar k/\textit{Wi}^{-3/2}$
for the lower-branch neutral curve of PDI instability in a viscoelastic PPF.
The properties of the most unstable regime-I PDI mode for different
$\beta$
values under BC1. Panels(a) and (b) show rescaled wavenumber
$\bar k=k\epsilon ^{1/2}$
for PPF (
$\lambda =2$
) and PCF (
$\lambda =1$
), respectively; (c) rescaled growth rate
$c_{1i}$
. The light dashed lines in (a, b) indicate the neutral curves extracted from figure 3.
Figures 4(a) and 4(b) illustrate the parameter pairs
$(\lambda \textit{Wi},\bar k)$
for the most unstable PDI mode in PPF and PCF, respectively. Analysing the results across different polymer concentrations, a noticeable decrease in
$\bar k$
with an increase in
$\textit{Wi}$
is evident, with only a limited discrepancy observed between PPF and PCF. Panel (c) shows the dependence of the rescaled growth rate
$c_{1i}$
on
$\lambda \textit{Wi}$
, where a scaling relation of
$c_{1i}\sim \textit{Wi}^{1/2}$
is clearly identified in the high-
$\textit{Wi}$
limit for each curve.
The observations above motivate us to investigate the instability mechanism associated with the lower-branch PDI mode in the high-
$\textit{Wi}$
limit, as this regime also encompasses the most unstable state. Given that the instability properties of PPFs and PCFs are quite close, we will focus on PPFs for demonstration in the following. In particular, we will focus on the scaling relations
$\bar k\sim \textit{Wi}^{-3/2}$
and
$c_{1i}\sim \textit{Wi}^{1/2}$
in the subsequent subsection.
3.3. Regime II analysis
In regime II, we consider a large Weissenberg number with
$\bar k\sim \textit{Wi}^{-3/2}$
. Thus, for convenience, we further rescale the wavenumber as
\begin{equation} \tilde k=\textit{Wi}^{3/2}\bar k=\textit{Wi}^{3/2}\epsilon ^{1/2}k=O(1). \end{equation}
Under this condition, the diffusive layer splits into two distinct regions: LDL and UDL, with the asymptotic structure illustrated in figure 1(b).
3.3.1. Lower diffusive layer for regime II
Balancing the conformation tensor and viscosity in the momentum equation (3.5a
), we obtain
$\hat u_0^{\prime\prime}\sim \sigma \textit{Wi}{\textrm {i}} \bar k\hat c_{110}$
, which leads to
$\hat c_{110}\sim Y^{-2}\textit{Wi}^{-1}\bar k^{-1}\hat u_0$
. Considering the conformation tensor equation (3.5d
), we need to balance the diffusive term (
$\hat c_{110}^{\prime\prime}$
) with the convective term (e.g.
$-4{\textrm {i}}\bar k\lambda ^2\hat u_0$
), which yields
$\hat c_{110}\sim Y^2\bar k\hat u_0$
. These two relations together determine the thickness of the diffusive layer, yielding
$Y\sim \textit{Wi}^{1/2}$
. For convenience, we introduce a local coordinate
\begin{equation} \tilde Y=\textit{Wi}^{-1/2}Y=\textit{Wi}^{-1/2}\epsilon ^{-1/2}(1+y)=O(1). \end{equation}
Simultaneously, we can derive the magnitude of the conformation tensor from the relations established above,
$\hat c_{110}\sim \textit{Wi}^{-1/2}\hat u_0$
. Similarly, we find that
$\hat v_0\sim \textit{Wi}^{-1}\hat u_0$
,
$\hat p\sim \textit{Wi}^{-3/2}\hat u_0$
,
$\hat c_{120}=\textit{Wi}^{-3/2}\hat u_0$
and
$\hat c_{220}=\textit{Wi}^{-5/2}\hat u_0$
. We can then expand the perturbation field as follows:
\begin{equation} (\hat u_0,\hat v_0,\hat p_0,\hat c_{110},\hat c_{120},\hat c_{220})=\left (\tilde u_0,\frac {\tilde v_0}{\textit{Wi}},\frac {\tilde p_0}{\textit{Wi}^{3/2}},\frac {\tilde c_{110}}{\textit{Wi}^{1/2}},\frac {\tilde c_{120}}{\textit{Wi}^{3/2}},\frac {\tilde c_{220}}{\textit{Wi}^{5/2}}\right )+\cdots . \end{equation}
By balancing the inertia term in the conformation tensor equation, we find that the expression for the phase speed expansion should be modified to
\begin{equation} c=U(-1)+\epsilon ^{1/2}\textit{Wi}^{1/2}\tilde c_1+\cdots . \end{equation}
It can be observed that the scaling of the phase speed for the regime-II analysis is consistent with that of the most unstable mode depicted in figure 4.
Substituting the above expansions into the instability equations (3.5), we obtain
\begin{align}&\quad\,\, \tilde u_0^{\prime\prime}+\sigma \big({\textrm {i}}\tilde k\tilde c_{110}+\tilde c_{120}^{\prime}\big)=0,\quad \tilde v_0^{\prime\prime}+ \sigma \big({\textrm {i}}\tilde k\tilde c_{120}+\tilde c_{220}^{\prime}\big)-\tilde p_0^{\prime}=0, \end{align}
\begin{align}&\qquad {\textrm {i}}\tilde k\tilde u_0+\tilde v_0^{\prime}=0,\quad \tilde c_{110}^{\prime\prime}=\tilde S\tilde c_{110}-2 \big(2{\textrm {i}} \tilde k\lambda ^2\tilde u_0+\lambda \tilde u_0^{\prime}+\lambda \tilde c_{120}\big), \end{align}
\begin{align}& \tilde c_{120}^{\prime\prime}= \tilde S\tilde c_{120}-2{\textrm {i}} \tilde k\lambda ^2\tilde v_0-\tilde u_0^{\prime}-\lambda \tilde c_{220},\quad \tilde c_{220}^{\prime\prime}=\tilde S\tilde c_{220}-2\big({\textrm {i}}\tilde k \lambda \tilde v_0+\tilde v_0^{\prime}\big), \end{align}
where
$\tilde S={\textrm {i}}\tilde k(\lambda \tilde Y-\tilde c_1)+1$
. It is seen that the streamwise pressure gradient is absent in the leading-order streamwise momentum equation in this regime. Under these scaling relations, the wall boundary conditions for both BC1 and BC2, as illustrated in (3.7) and (3.8), respectively, can be recast into a unified form
\begin{align} \tilde u_0(0)=\tilde v_0(0)=0,\quad { \tilde S_0}\tilde c_{110}(0)={2\lambda \left[\tilde u_0^{\prime}(0)+\tilde c_{120}(0)\right]},\quad { \tilde S_0}\tilde c_{120}(0)={\tilde u_0^{\prime}}(0),\quad \tilde c_{220}=0, \end{align}
where
$\tilde S_0=\tilde S(0)=-{\textrm {i}} \tilde k\tilde c_1+1$
. Therefore, we use BC12 to represent this type of boundary condition. For BC3, we can express the wall boundary conditions as
\begin{align} \tilde u_0(0)=\tilde v_0(0)=\tilde c_{110}^{\prime}(0)=\tilde c_{120}^{\prime}(0)=\tilde c_{220}^{\prime}(0)=0. \end{align}
This system does not guarantee the attenuation conditions in the limit of
$\tilde Y\to \infty$
. Therefore, we must consider a thicker layer (UDL) where the streamwise pressure gradient comes back to the leading order. By analysing the balance of the momentum equation, we find that the thickness of the UDL should be comparable to the rescaled streamwise length scale, specifically,
$(1+y)\sim \bar k\sim \textit{Wi}^{-3/2}$
.
3.3.2. Upper diffusive layer for
$1-\beta =O(1)$
Following the aforementioned scaling estimate, we introduce
\begin{equation} \bar Y=\textit{Wi}^{-1}\tilde Y=\textit{Wi}^{-3/2}Y=\textit{Wi}^{-3/2}\epsilon ^{-1/2}(1+y)=O(1). \end{equation}
From a scaling estimate, we expand the perturbation field in this layer as
\begin{equation} (\hat u_0,\hat v_0,\hat p_0,\hat c_{110},\hat c_{120},\hat c_{220})=\left (\bar u_0,\bar v_0,\frac {\bar p_0}{\textit{Wi}^{1/2}},\frac {\bar c_{110}}{\textit{Wi}^{3/2}},\frac {\bar c_{120}}{\textit{Wi}^{3/2}},\frac {\bar c_{220}}{\textit{Wi}^{5/2}}\right )+\cdots . \end{equation}
Substituting these equations to (3.5), we obtain the leading-order governing equations in this layer
\begin{align}&\quad \sigma \big({\textrm {i}}\tilde k\bar c_{110}+\bar c_{120}^{\prime}\big)-{\textrm {i}}\tilde k\bar p_0=0,\quad {\textrm {i}}\sigma \tilde k\bar c_{120}-\bar p_0^{\prime}=0, \end{align}
\begin{align}& {\textrm {i}}\tilde k\bar u_0+\bar v_0^{\prime}=0,\quad {\textrm {i}}\bar k\lambda \bar Y\bar c_{110}-2 \big(2{\textrm {i}} \tilde k\lambda ^2\bar u_0+\lambda \bar c_{120}\big)=0, \end{align}
\begin{align}&\,\,\,\, {\textrm {i}}\bar k\lambda \bar Y\bar c_{120}-2{\textrm {i}} \tilde k\lambda ^2\bar v_0=0,\quad {\textrm {i}}\bar k\lambda \bar Y\bar c_{220}-2{\textrm {i}}\tilde k \lambda \bar v_0=0. \end{align}
The system can be reduced to
\begin{equation} \bar Y^2\bar v_0^{\prime\prime}-2\bar Y\bar v'_0+(2-\tilde k^2\bar Y^2)\bar v_0=0, \end{equation}
whose solution, by excluding the exponentially growing part, is expressed as
\begin{equation} \bar v_0=d_0\bar Y{\textrm {e}}^{-\tilde k\bar Y}, \end{equation}
with
$d_0$
being a constant. We can therefore obtain the other perturbation quantities from (3.19)
\begin{equation} \left (\bar u_0,\bar c_{110},\bar c_{120},\bar c_{220}\right )=\big({{\textrm {i}}{\bar k}^{-1}(1-\bar k\bar Y)},-4{\textrm {i}}\lambda ,2\lambda ,2\big)d_0{\textrm {e}}^{-\bar k\bar Y}. \end{equation}
It is obvious that the perturbations attenuate exponentially in the large-
$\bar Y$
limit. We can also estimate the asymptotic behaviours of these perturbations in the lower limit
\begin{equation} \left (\bar u_0, \bar v_0, \bar p_0,\bar c_{110},\bar c_{120},\bar c_{220}\right )\to \big({{\textrm {i}} }{\bar k}^{-1},\bar Y, -2{\textrm {i}}\lambda \sigma , -4{\textrm {i}}\lambda , 2\lambda , 2\big)\quad \mbox{as }\bar Y\to 0. \end{equation}
Comparing the expansions between (3.18) with (3.12), we find that the magnitude of
$\hat c_{110}$
in this layer is smaller than that in LDL, indicating that
$\tilde c_{110}$
in LDL should approach zero in the large-
$\tilde Y$
limit. The other perturbations should match in the overlapping region between the two layers. Thus, we can obtain the upper boundary conditions for system (3.14)
\begin{equation} \left (\tilde u_0^{\prime}, \tilde c_{110},\tilde c_{120}^{\prime}, \tilde c_{220}^{\prime}\right )\to 0,\quad \tilde p_0\to -2\sigma \bar k\lambda \tilde u_0\quad \mbox{as }\tilde Y\to \infty . \end{equation}
3.3.3. Numerical results for regime II
In regime II, the eigenvalue system is reformulated as (3.14) with boundary conditions (3.15) (or (3.16)) and (3.24), where
$\tilde c_1$
is the eigenvalue. The number of the controlling parameters is further reduced to two:
$\tilde k$
and
$\beta$
. Choosing a few representative
$\beta$
values, we plot the dependence of the rescaled growth rate
$\tilde c_{1i}$
on the rescaled wavenumber
$\tilde k$
for a PPF (
$\lambda =2$
), as shown in figure 5. The results for both boundary conditions (3.15) and (3.16) are examined. For each curve, the PDI becomes unstable when the rescaled wavenumber
$\tilde k$
exceeds a critical value, known as the neutral point. The rescaled growth rate
$\tilde c_{1i}$
increases with
$\tilde k$
until reaching its most unstable state. However, for further higher
$\tilde k$
values, the rescaled growth rate decreases monotonically. For BC12, represented by (3.15), the neutral wavenumber
$\tilde k$
shifts to higher values as
$\beta$
approaches unity, while
$\tilde c_{1i}$
decreases notably. In contrast, for BC3, represented by (3.16), the variation of
$\tilde k$
on
$\beta$
is not monotonic. Notably, for
$\tilde k\gg 1$
, the rescaled growth rate
$\tilde c_{1i}$
for BC12 reduces as the dilute limit is approached, while the opposite is true for BC3. In the inset of each panel, we apply a suitable rescaling so that the curves tend to collapse in the dilute limit,
$\beta \to 1$
. Remarkably, for BC12, the growth rate
$\tilde c_{1i}$
scales as
$\tilde k^{-1}$
, while for BC3,
$\tilde c_{1i}$
scales as
$\tilde k^{-1/3}$
. On the other hand, both conditions share the same horizontal scaling parameter,
$\tilde k(1-\beta )^{3/2}$
, implying the fundamental scaling relation
$\tilde k\sim (1-\beta )^{-3/2}$
in the dilute limit. Motivated by Dong & Zhang (Reference Dong and Zhang2022), in which the centre-mode instability system is further reduced by taking the dilute limit,
$\beta \to 1$
, we explore PDI in the dilute limit in the following subsection, referring this specific scaling as regime III.
Dependence of the rescaled growth rate
$\tilde c_{1i}$
on the rescaled wavenumber
$\tilde k$
with
$\beta =0.7$
, 0.8, 0.9, 0.95 and 0.98 for regime-II PDI mode for a PPF. The inset in each panel shows the curves obtained after applying an appropriate rescaling: (a) BC12; (b) BC3.
3.4. Regime III analysis
3.4.1. Formulation for regime III
Now we take
$\beta \to 1$
, and introduce
$\sigma =(1-\beta )/\beta \ll 1$
. Then, from the scaling relation revealed by figure 5, we set
$\tilde k\sim \sigma ^{-3/2}$
and introduce
\begin{equation} \breve \sigma =\tilde k^{2/3}\sigma =O(1). \end{equation}
From a scaling estimate, we find that the thickness of LDL becomes
$\tilde Y\sim \tilde k^{-1/3}$
. Thus, we introduce a local coordinate
\begin{equation} \breve Y=\tilde k^{1/3}\tilde Y=O(1), \end{equation}
and expand the perturbation field as
\begin{align}&\qquad\qquad\quad (\hat u_0,\hat v_0,\hat p_0)=\big(\breve U_0,\textit{Wi}^{-1}\tilde k^{2/3}\breve V_0,\textit{Wi}^{-3/2}\tilde k\breve P_0\big)+\cdots , \end{align}
\begin{align}& (\hat c_{110},\hat c_{120},\hat c_{220})=\big(\textit{Wi}^{-1/2}\tilde k^{1/3}\breve C_{110},\textit{Wi}^{-3/2}\tilde k\breve C_{120},\textit{Wi}^{-5/2}\tilde k\breve C_{220}\big) +\cdots . \\[9pt] \nonumber \end{align}
Examining the boundary conditions (3.15c
) for BC12, we find that only when
$\tilde c_1\sim \tilde k^{-1}$
can these boundary conditions appear at leading-order impact. This scaling is selected as the magnitude of the phase speed correction for BC12. In contrast, for BC3, the boundary conditions (3.16) always appear in leading order, and to construct an eigenvalue system, the phase speed correction must appear in the governing equations with
$\tilde c_1\sim \tilde k^{-1/3}$
. These estimates also align with the numerical observations from figure 5 by applying
$\beta \to 1$
. Therefore, we express the phase speed expansion in a unified form
\begin{equation} c=U(-1)+\epsilon ^{1/2}\textit{Wi}^{1/2}\tilde k^{\alpha }\breve c_1 +\cdots , \end{equation}
where
$\alpha =-1$
for BC12 and = −1/3 for BC3.
Substituting the above expansions into (3.14) and collecting the leading-order terms, we obtain
\begin{align}& \breve U_0^{\prime\prime}+\breve \sigma \big({\textrm {i}}\breve C_{110}+\breve C_{120}^{\prime}\big)=0,\quad \breve V_0^{\prime\prime}+\breve \sigma {\textrm {i}}\breve C_{120}-\breve P_0^{\prime}=0, \end{align}
\begin{align}& {\textrm {i}}\breve U_0+\breve V_0^{\prime}=0,\quad \breve C_{110}^{\prime\prime}=\breve S\breve C_{110}-2\big(2{\textrm {i}} \lambda ^2\breve U_0+\lambda \breve C_{120}\big), \end{align}
\begin{align}&\,\, \breve C_{120}^{\prime\prime}= \breve S\breve C_{120}-2{\textrm {i}} \lambda ^2\breve V_0,\quad \breve C_{220}^{\prime\prime}=\breve S\breve C_{220}-2{\textrm {i}}\lambda \breve V_0. \end{align}
where
$\breve S={\textrm {i}}\lambda \breve Y$
for BC12 and
$\breve S={\textrm {i}}(\lambda \breve Y-\breve c_1)$
for BC3.
For BC12, the phase speed correction
$\breve c_1$
does not appear in the equation system (3.29); it only appears in the wall boundary conditions
\begin{align} \breve U_0(0)=\breve V_0(0)=\breve C_{120}(0)=\breve C_{220}(0)=0,\quad {\breve S_0} \breve C_{110}(0)={2\lambda \breve U_0^{\prime}(0)}. \end{align}
where
$\breve S_0=-{\textrm {i}} \breve c_1+1$
. Note that the boundary condition for
$\tilde c_{120}$
, (3.15d
), comes to the second order. In contrast,
$\breve c_1$
appears in the leading-order equation system for BC3, and the boundary conditions are
\begin{align} \breve U_0(0)=\breve V_0(0)=\breve C_{110}^{\prime}(0)=\breve C_{120}^{\prime}(0)=\breve C_{220}^{\prime}(0)=0. \end{align}
From a scaling analysis, we obtain that the UDL thickness in this regime becomes
$\bar Y\sim \tilde k^{-1/3}$
. Following the same mathematical procedure, we obtain the upper boundary conditions of the LDL perturbations from matching with the UDL,
\begin{equation} \big(\breve U_0^{\prime}, \breve C_{110},\breve C_{120}^{\prime}, \breve C_{220}^{\prime}\big)\to 0,\quad \breve P_0\to -2\breve \sigma \lambda \breve U_0\quad \mbox{as }\breve Y\to \infty . \end{equation}
Eliminating
$\hat u$
,
$\hat p$
and
$\hat {\boldsymbol c}$
governing equations (3.29), we obtain a seventh-order ordinary differential equation for
$\breve V_0$
\begin{equation} \big[{\mathcal D}^4-\breve S{\mathcal D}^2+{\textrm {i}}\lambda {\mathcal D}+2\lambda ^2\breve \sigma\big]\big[{\mathcal D}^3-\breve S{\mathcal D}+{\textrm {i}}\lambda\big]\breve V_0=0, \end{equation}
where
${\mathcal D}:=d/{\rm d}\breve Y$
. The no-slip, non-penetration boundary conditions are readily rewritten as
\begin{equation} \breve V_0(0)=\breve V_0^{\prime}(0)=0. \end{equation}
To convert the boundary conditions of the conformation tensor, we need to observe the near-wall behaviours of the perturbation field.
Let us expand the perturbation field in the near-wall region
\begin{align}&\quad \breve C_{110}=\bar a_0+\bar a_1\breve Y+\bar a_2\frac {\breve Y^2}2+\bar a_3\frac {\breve Y^3}6+\cdots , \end{align}
\begin{align}& \breve C_{120}=\bar b_0+\bar b_1\breve Y+\bar b_2\frac {\breve Y^2}2+\bar b_3\frac {\breve Y^3}6+\bar b_4\frac {\breve Y^4}{24}\cdots , \end{align}
where
$(\bar a_0,\bar a_1,\bar a_2,\bar a_3,\bar b_0,\bar b_1,\bar b_2,\bar b_3,\bar b_4)$
are constants.
For BC12, we have
$\breve S(0)=\breve C_{120}(0)=0$
, then, it is easy to derive that
$\bar a_2=\bar b_0=\bar b_2=\bar b_3=0$
. Substituting (3.35) into (3.29) and (3.30e
), we obtain a set of algebraic equations
\begin{align}& {\textrm {i}}\breve V_0^{(3)}(0)+\breve \sigma ({\textrm {i}}\bar a_0+\bar b_1)=0,\quad {\textrm {i}} \breve V_0^{(4)}(0)+{\textrm {i}}\breve \sigma \bar a_1=0,\quad {\textrm {i}}\breve V_0^{(5)}(0))=0, \end{align}
\begin{align}&\quad {\textrm {i}}\breve V_0^{(6)}(0)+\breve \sigma ({\textrm {i}} \bar a_3+\bar b_4)=0,\quad \bar a_3={\textrm {i}}\lambda \bar a_0+4\lambda ^2 {\textrm {i}}\breve V_0^{\prime\prime}(0)-2\lambda \bar b_1, \end{align}
\begin{align}&\qquad\quad \bar b_4=2{\textrm {i}}\lambda \bar b_1-2{\textrm {i}}\lambda ^2 \breve V_0^{\prime\prime}(0),\quad (1-{\textrm {i}}\breve c_1)\bar a_0=2{\textrm {i}}\lambda \breve V_0^{\prime\prime}(0). \end{align}
Eliminating the
$(\bar a_0,\bar a_1,\bar a_3,\bar b_1,\bar b_4)$
, we obtain the conformation boundary conditions to (3.33) for BC1
\begin{equation} \breve V_0^{(5)}(0)=0,\quad (1-{\textrm {i}}\breve c_1)\breve V_0^{(6)}(0)-2{\textrm {i}}\breve c_1 \lambda ^2\breve \sigma \breve V_0^{\prime\prime}(0)=0. \end{equation}
For BC3, following the same procedure, we derive
\begin{align} \breve V_0^{(5)}(0)+{\textrm {i}}\breve c_1\breve V_0^{(3)}(0)-\frac {\lambda \breve V_0^{(4)}(0)}{\breve c_1}&=0,\nonumber\\ \breve V_0^{(6)}(0)-{\textrm {i}}\lambda \breve V_0^{(3)}(0)+{\textrm {i}}\breve c_1\breve V_0^{(4)}(0)+2\lambda ^2\breve \sigma \breve V_0^{\prime\prime}(0)&=0. \end{align}
Additionally, the upper boundary conditions (3.32) can be converted to
\begin{equation} \big(\breve V_0^{\prime\prime},\breve V_0^{(3)},\breve V_0^{(4)}\big)\to 0\quad \mbox{as }\breve Y\to \infty . \end{equation}
3.4.2. Numerical results for regime III
In this regime, the eigenvalue system is characterised by (3.33), complemented by the boundary conditions (3.34), (3.37) (or (3.38)) and (3.39). The eigenvalue
$\breve c_1$
is solely governed by a single parameter,
$\breve \sigma$
, whose real and imaginary parts for PPF are plotted in figure 6. For BC12, the PDI mode becomes unstable when
$\breve \sigma$
exceeds 0.86, with its rescaled growth rate
$\breve c_{1i}$
increasing monotonically with
$\breve \sigma$
. The rescaled phase speed correction
$\breve c_{1r}$
initially declines with increasing
$\breve \sigma$
, although a reversal occurs when
$\breve \sigma \gt 3.3$
. In contrast, for BC3, the neutral point emerges at
$\breve \sigma =0.0735$
, with the rescaled growth rate peaking at
$\breve \sigma =0.416$
; the rescaled phase speed decreases monotonically with increasing
$\breve \sigma$
. It is crucial to emphasise that despite the magnitude of the rescaled
$\breve c_1$
for BC12 is greater than that for BC3, it does not mean that the PDI for BC12 is more unstable, due to the differing scaling laws for the two phase speed corrections, varying by a factor
$\tilde k^{-2/3}$
. The insets in figures 6(a) and 6(b) further present log–log plots, clearly revealing an excellent scaling relationship between
$\breve c_1$
and
$\breve \sigma$
in the limit of
$\breve \sigma \to \infty$
. Motivated by this observation, we proceed to simplify the instability system further by taking the asymptotic limit of removing the only remaining controlling parameter. This approach leads to a parameter-free system, allowing the analytical derivation of solutions. We will investigate this simplified system in the following subsection, referring to it as regime IV.
Dependence on
$\breve \sigma$
of the real and imaginary parts of
$\breve c_{1}$
for regime III with BC12 (a) and BC3 (b), where
$\lambda =2$
(PPF). The circle in each panel marks the neutral point. The inset shows the log–log plot to illustrate the large-
$\breve \sigma$
asymptote, with the squares in (a) and (b) denoting the asymptotic predictions by (3.56) and (3.67), respectively.
The curves in figure 7 illustrate the variation with
$\tilde k$
of the rescaled growth rate
$\tilde c_{1i}$
predicted by the regular-concentration regime II discussed in § 3.3. Specifically, we choose
$\breve \sigma =\sigma /\tilde k^{2/3}=3.5$
to match the dilute-limit scaling, where
$\tilde k\to \infty$
signifies the dilute limit
$\sigma \to 0$
. For both BC12 and BC3, a sharp increase in
$\tilde c_{1i}$
with
$\tilde k$
is evident near the lower-branch neutral point, followed by a gradual decrease with a scaling
$\tilde c_{1i}\sim \tilde k^{-1}$
, which represents the dilute-limit scaling. The neutral wavenumber
$\tilde k$
is found to be greater for BC3. In the high-
$\tilde k$
limit, the regime-II calculations of
$\tilde c_{1i}$
converge toward the predictions for regime III, as indicated by the circles, suggesting a natural connection between the two regimes. When comparing the results of the two boundary conditions, it is evident that, although the growth rate for BC12 at a moderate
$\tilde k$
is higher than that for BC3, it decreases at a much higher rate as
$\tilde k$
increases.
Comparison of the rescaled growth rate
$\tilde c_{1i}$
in the dilute and high-
$\tilde k$
limits obtained by regime II and III calculations, where
$\breve \sigma =\tilde k^{2/3}\sigma =3.5$
and
$\lambda =2$
(PPF): (a) BC12; (b) BC3.
Comparison of the eigenfunctions between the regime-II and regime-III calculations in the dilute limit, where
$\breve \sigma =3.5$
and
$\lambda =2$
(PPF). For regime II, we choose
$\beta =0.997$
(
$\sigma =0.003$
) and
$\tilde k=3.97\times 10^4$
. Solid lines: regime II; circles: regime III: (a) BC12; (b) BC3. We choose
$\breve U_0^{\prime}(0)=1$
for normalisation.
In figure 8, we compare the norms of the eigenfunctions obtained from calculations based on regimes II and III. The dilute limit with
$\beta =0.997$
(or
$\sigma =0.003$
) and
$\breve \sigma =3.5$
is considered. Both boundary conditions, BC12 and BC3, exhibit similar near-wall behaviours of the perturbation velocity field, specifically:
$\breve V_0\sim \breve Y^2$
,
$\breve U_0\sim \breve Y$
. However, the near-wall conformation perturbations show evident discrepancy:
$\breve C_{120}$
for BC12 is almost proportional to
$\breve Y$
, while for BC3, it exhibits a zero gradient. The alignment of the calculations from the two regimes further validates the accuracy of the predictions made by our asymptotic analyses.
3.5. Regime IV analysis
As shown in figure 6, the system (3.33) becomes unstable when
$\breve \sigma$
exceeds a certain threshold, marked by a dot, and a clear scaling relation between
$\breve c_{1}$
and
$\breve \sigma$
is evident as
$\breve \sigma \to \infty$
. To facilitate an analytical derivation of the dispersion relation by further simplifying the instability formulation, we perform an asymptotic analysis in this subsection by taking
$\breve \sigma$
to be asymptotically large. In this regime, no controlling parameter remains. A scaling estimate determines that the LDL further splits into two distinct layers, a bulk sublayer with
$\breve Y=O(1)$
and a wall sublayer with
$\breve Y\sim \breve \sigma ^{-1/4}$
.
3.5.1. Regime-IV solution for BC12
For BC12, figure 6(a) reveals that both the real and imaginary parts of
$\breve c_1$
increase with
$\breve \sigma$
, and fitting the data show that
$\breve c_{1}\sim \breve \sigma ^{1/2}$
as
$\breve \sigma \to \infty$
. Therefore, we expand the phase speed as
\begin{equation} \breve c_{1}=\breve \sigma ^{1/2}(\breve c_{10}+\breve \sigma ^{-1/4}\breve c_{11}+\cdots ). \end{equation}
In the bulk sublayer, the solution is expanded as
\begin{equation} \breve V_0=\sum _{l=0}^{\cdots }\breve \sigma ^{-l/4}\breve V_{0l}. \end{equation}
Substituting into (3.33), we obtain the governing equations to the first three orders
\begin{align} {\mathcal L}\breve V_{00}=0,\quad {\mathcal L}\breve V_{01}=0,\quad {\mathcal L}\breve V_{02}=0, \end{align}
where
${\mathcal L}:=[{\mathcal D}^3-{\textrm {i}}\lambda \breve Y{\mathcal D}+{\textrm {i}}\lambda ]$
. Excluding the exponentially growing component, we derive the solutions
\begin{equation} \breve V_{0l}= d_{1l}\tilde \xi +d_{2l}\left (\tilde \xi \int _\infty ^{\tilde \xi }\mbox{Ai}(\tilde \xi ){\rm d}\tilde \xi -\mbox{Ai}^{\prime}(\tilde \xi )\right )\! ,\quad \mbox{with }l=0,1,2, \end{equation}
where Ai denotes the Airy function of the first kind (Abramowitz & Stegun Reference Abramowitz and Stegun1964),
$\tilde \xi =({\textrm {i}}\lambda )^{1/3}\breve Y$
and
$(d_{1l},d_{2l})$
are constants to be determined. Here, the exponentially growing part is dropped to satisfy the matching condition with the UDL solution. As the wall is approached, the asymptotic behaviour of
$\breve V_{0l}$
is
\begin{equation} \breve V_{0l}\to {\mathcal V}_{l0}+{\mathcal V}_{l1}\breve Y+{\mathcal V}_{l2}\breve Y^2/2+\cdots \quad \mbox{as }\breve Y\to 0, \end{equation}
where
\begin{equation} {\mathcal V}_{l0}=-d_{2l}\mbox{Ai}^{\prime}(0),\quad {\mathcal V}_{l1}=(d_{1l}-d_{2l}/3)({\textrm {i}}\lambda )^{1/3},\quad {\mathcal V}_{l2}=d_{2l}\mbox{Ai}(0)({\textrm {i}}\lambda )^{2/3}. \end{equation}
For a non-trivial expansion, we assume that
$d_{10}$
and
$d_{20}$
cannot both be zero. Since the conditions in (3.34) are not satisfied for this case, we need to consider the underneath wall sublayer.
In the wall sublayer, we introduce
$\eta =\breve \sigma ^{1/4}\breve Y=O(1)$
, and expand the solution as
\begin{equation} \breve V_0=\sum _{l=0}^{\cdots }\breve \sigma ^{-l/4}\breve v_l. \end{equation}
The governing equations to the first three orders read
\begin{align} \big(\bar {\mathcal D}^4+2\lambda ^2\big)\bar {\mathcal D}^3\breve v_{l}=0\quad \mbox{with }l=0,1,2, \end{align}
where
$\bar {\mathcal D}:=d/{\rm d}\eta$
. The wall boundary conditions read
\begin{align}&\qquad\quad \breve v_{l}(0)=\breve v_{l}^{\prime}(0)=\breve v_{l}^{(5)}(0)=0\quad \mbox{with }l=0,1,2, \end{align}
\begin{align}& \breve v_{0,1}^{(6)}(0)+2\lambda ^2\breve v_{0,1}^{\prime\prime}(0)=0,\quad \breve v_{2}^{(6)}(0)+2\lambda ^2\breve v_{2}^{\prime\prime}(0)=\breve v_0^{(6)}(0)/({\textrm {i}}\breve c_{10}). \end{align}
The leading-order solution is
\begin{equation} \breve v_0={\mathcal V}_{00}+a_0\,{\textrm {e}}^{\bar \lambda _1\eta }+b_0\,{\textrm {e}}^{\bar \lambda _2\eta }, \end{equation}
where
$\bar \lambda _{1,2}=(-1\pm {\textrm {i}})(2\lambda ^2)^{1/4}/\sqrt {2}$
, and
$(a_0,b_0)$
are constants to be determined. To match the bulk sublayer, the exponentially growing complementary solution is excluded. Applying the boundary conditions at the wall, it is determined that
\begin{align} &{\mathcal V}_{00}+a_0+b_0=0,\quad \bar \lambda _1a_0+\bar \lambda _2b_0=0,\quad \bar \lambda _1^5 a_0+\bar \lambda _2^5b_0=0,\nonumber\\&\quad \bar \lambda _1^6 a_0+\bar \lambda _2^6b_0+2\lambda ^2(\bar \lambda _1^2 a_0+\bar \lambda _2^2b_0)=0. \end{align}
The equation (3.51d ) is satisfied automatically, and the conditions (3.51b ) and (3.51c ) are equivalent. Then, the solutions to (3.51a , b ) are
\begin{equation} a_0=\frac {-\bar \lambda _2{\mathcal V}_{00}}{\bar \lambda _2-\bar \lambda _1},\quad b_0=\frac {\bar \lambda _1{\mathcal V}_{00}}{\bar \lambda _2-\bar \lambda _1}. \end{equation}
Considering the matching condition, the second-order solution is
\begin{equation} \breve v_1={\mathcal V}_{10}+{\mathcal V}_{01}\eta +a_1{\textrm {e}}^{\bar \lambda _1\eta }+b_1{\textrm {e}}^{\bar \lambda _2\eta }. \end{equation}
Applying the wall boundary conditions, we obtain that
${\mathcal V}_{01}$
must be zero (the solutions for
$a_1$
and
$b_1$
are omitted for brevity because they do not affect the determination of the dispersion relation).
The third-order solution is
\begin{equation} \breve v_2={\mathcal V}_{20}+{\mathcal V}_{11}\eta +{\mathcal V}_{02}\eta ^2/2+a_2{\textrm {e}}^{\bar \lambda _1\eta }+b_2{\textrm {e}}^{\bar \lambda _2\eta }. \end{equation}
Applying the wall boundary condition (3.49b ), we obtain
\begin{align} {\mathcal V}_{20}+a_2+b_2=0,\quad {\mathcal V}_{11}+\bar \lambda _1a_2+\bar \lambda _2b_2=0,\quad \bar \lambda _1^5 a_2+\bar \lambda _2^5b_2=0, \end{align}
\begin{align} \bar \lambda _1^6 a_2+\bar \lambda _2^6b_2+2\lambda ^2\big(\bar \lambda _1^2 a_2+\bar \lambda _2^2b_2+{\mathcal V}_{02}\big)=\big(\bar \lambda _1^6a_0+\bar \lambda _2^6b_0\big)/({\textrm {i}}\breve c_{10}). \end{align}
Combining (3.52) and (3.55d
), and considering
${\mathcal V}_{02}/{\mathcal V}_{00}=-({\textrm {i}}\lambda )^{2/3}\mbox{Ai}(0)/\mbox{Ai}^{\prime}(0)$
, we derive
\begin{equation} \breve c_{10} =-\frac {a_0\bar \lambda _1^2+b_0\bar \lambda _2^2}{{\textrm {i}} {\mathcal V}_{02}}=\left \{\!\!\begin{array}{ll} 1.1249+0.6495{\textrm {i}}&\mbox{for }\lambda =2;\\ 0.8929 + 0.5155{\textrm {i}}&\mbox{for }\lambda =1. \end{array} \right . \end{equation}
Here, the positive imaginary part corresponds to the unstable PDI. The asymptotic prediction (3.56) is represented by the squares in figure 6(a), and excellent agreement with the numerical solutions is achieved in the limit of
$\breve \sigma \to \infty$
.
3.5.2. Regime-IV solution for BC3
For BC3, on the contrary, the complex phase speed
$\breve c_1$
decreases with
$\breve \sigma$
at a rate of
$\breve c_1\sim \breve \sigma ^{-1/2}$
. Consequently, the expansion in (3.40) is changed to
\begin{equation} \breve c_{1}=\breve \sigma ^{-1/2}\big(\breve c_{10}+\breve \sigma ^{-1/4}\breve c_{11}+\cdots\big). \end{equation}
In the bulk sublayer, the expansion of
$\breve V_0$
remains identical to that given by (3.41), and the first two orders of the governing equations are unchanged from (4.3a
,
b
). Therefore, the solutions for
$\breve V_{00}$
and
$\breve V_{01}$
are the same as those presented in (3.43), with the asymptotic behaviours in the limit of
$\breve Y\to 0$
following directly from (3.44) with
$l=0,1$
.
In the wall sublayer, the leading-order solution takes the same form as in (3.50). Applying boundary conditions (3.51a
,
b
),
$\breve v_0^{(4)}=0$
and
$\breve v_0^{(6)}+2\lambda ^2\breve v_0^{\prime\prime}=0$
, we find that
${\mathcal V}_{00}=0=a_0=b_0$
, implying
$\breve v_0=0$
. Consequently, using (3.45), we obtain
$d_{20}=0$
,
${\mathcal V}_{01}=d_{10}({\textrm {i}}\lambda )^{1/3}$
, and thus,
$\breve V_{00}=d_{10}({\textrm {i}}\lambda )^{1/3}\breve Y$
.
Following the same procedure, we can determine the coefficients to second-order solution (3.53):
${\mathcal V}_{10}=0$
and
\begin{equation} a_1=-b_1={\mathcal V}_{01}/(\bar \lambda _2-\bar \lambda _1). \end{equation}
This also implies
$d_{21}=0$
and
$\breve V_{01}=d_{11}({\textrm {i}}\lambda )^{1/3}\breve Y$
.
However, due to the presence of
$-{\textrm {i}}\breve c_1$
in
$\breve S$
for BC3, the third-order governing equations include the impact of the rescaled phase speed, expressed as
\begin{equation} {\mathcal L}\breve V_{02}=-{\textrm {i}}\breve c_{10}{\mathcal D}\breve V_{00}=-({\textrm {i}}\lambda )^{1/3}{\textrm {i}}\breve c_{10}d_{10}. \end{equation}
The solution can be expressed as
\begin{equation} \breve V_{02}=-\frac {{\textrm {i}}\breve c_{10}d_{10}}{({\textrm {i}}\lambda )^{2/3}}+ d_{12}\tilde \xi +d_{22}\left (\tilde \xi \int _\infty ^{\tilde \xi }\mbox{Ai} (\tilde \xi ){\rm d}\tilde \xi -\mbox{Ai}^{\prime}(\tilde \xi )\right )\!. \end{equation}
Thus, the coefficients in the asymptotic behaviour (3.44) for
$l=2$
read
\begin{equation} {\mathcal V}_{20}=-\frac {{\textrm {i}}\breve c_{10}d_{10}}{({\textrm {i}}\lambda )^{2/3}}-d_{22}\mbox{Ai}^{\prime}(0),\quad {\mathcal V}_{21}=\left (d_{12}-d_{22}/3\right )({\textrm {i}}\lambda )^{1/3},\quad {\mathcal V}_{22}=d_{22}\mbox{Ai}(0)({\textrm {i}}\lambda )^{2/3}. \end{equation}
In the wall sublayer, the third-order governing equation remains identical to (3.47), and so the solution can be expressed in the same form as (3.54). However, the coefficients
$a_2$
and
$b_2$
change, since the lower boundary conditions are now modified to
\begin{align} {\mathcal V}_{20}+a_2+b_2=0,\quad {\mathcal V}_{11}+\bar \lambda _1a_2+\bar \lambda _2 b_2=0,\quad \lambda {\mathcal V}_{20}=\breve c_{10}{\mathcal V}_{01}. \end{align}
The fourth-order solution
$\breve v_3$
does not provide valuable information for determination of the dispersion relation, and so we directly move to the fifth order. The governing equation is
\begin{equation} \big(\bar {\mathcal D}^4+2\lambda ^2\big)\bar {\mathcal D}^3\breve v_{4}=-2{\textrm {i}}\lambda ^3{\mathcal V}_{10}-2{\textrm {i}}\lambda ^3\big[a_1{\textrm {e}}^{\bar \lambda _1\eta }(3+\bar \lambda _1\eta ) +b_1{\textrm {e}}^{\bar \lambda _2\eta }(3+\bar \lambda _2\eta )\big], \end{equation}
which is subject to the wall boundary conditions
\begin{equation} \breve v_4(0)=\breve v_4'(0)=0,\quad \lambda \breve v_4^{(4)}(0)=\breve c_{10}\breve v_3^{(5)}(0)+\cdots ,\quad \breve v_4^{(6)}(0)+2\lambda ^2\breve v_4^{\prime\prime}(0)={\textrm {i}}\lambda \breve v_1^{(3)}, \end{equation}
where
$\cdots$
denotes terms associating with high-order expansions of
$\breve c_{1}$
. The solution can be expressed as
\begin{align} \displaystyle \breve v_4&=\displaystyle {\mathcal V}_{40}+{\mathcal V}_{31}\eta +{\mathcal V}_{22}\frac {\eta ^2}2+a_4{\textrm {e}}^{\bar \lambda _1\eta }+b_4{\textrm {e}}^{\bar \lambda _2\eta }\nonumber \\ &\quad\displaystyle +{\textrm {i}}\lambda \left [-{\mathcal V}_{10}\frac {\eta ^3}{6}+\frac {a_1{\textrm {e}}^{\bar \lambda _1\eta }}{8\bar \lambda _1^2} (\bar \lambda _1\eta ^2-3\eta ) +\frac {b_1{\textrm {e}}^{\bar \lambda _2\eta }}{8\bar \lambda _2^2}(\bar \lambda _2\eta ^2-3\eta )\right ]. \end{align}
Substituting into the fourth relation in (3.64), and considering the solutions (3.58), we obtain
\begin{equation} {\mathcal V}_{22}=\frac {{\textrm {i}}}{\sqrt {2}}{\mathcal V}_{01}. \end{equation}
Combining with (3.62c ) and (3.61), we can derive
\begin{equation} \breve c_{10}=-\frac {({\textrm {i}}\lambda )^{1/3}\mbox{Ai}^{\prime}(0)}{2\sqrt {2}\mbox{Ai}(0)}=\left \{\!\!\begin{array}{ll} 0.2812+0.1624{\textrm {i}}&\mbox{for }\lambda =2,\\ 0.2232 + 0.1289{\textrm {i}}&\mbox{for }\lambda =1. \end{array} \right . \end{equation}
Again the positive imaginary part of
$\breve c_{10}$
implies an unstable PDI. This asymptotic prediction is plotted in figure 6(b), which aligns with the numerical curves of regime III when
$\breve \sigma$
is asymptotically large.
4. Development of new conformation boundary conditions to eliminate unstable PDI
4.1. Derivation of the new conformation boundary conditions based on regime IV
Thanks to the previous asymptotic analyses, the original system (2.8) and (2.9), initially characterised by five controlling parameters, has been substantially simplified into regime IV, which features no remaining controlling parameters. In this subsection, we aim to construct a stable asymptotic solution based on regime IV, and derive a new set of conformation boundary conditions from the near-wall behaviours of this asymptotic solution.
4.1.1. Mathematical details
We first summarise the results obtained in § 3.5. For BC12, the phase speed
$\breve c_1$
does not appear in the leading-order equations within the bulk sublayer, and consequently, the dispersion relation is largely determined by the wall sublayer, subject to the relation between the constants
${\mathcal V}_{00}$
and
${\mathcal V}_{02}$
imposed by the bulk-sublayer solution. In principle, many choices are possible for the set of conformation boundary conditions at the wall; however, the present set promotes an unstable PDI and is therefore unfavourable. In contrast, for BC3, the rescaled phase speed contributes to both the bulk and wall sublayers, although it emerges only at third order in the bulk sublayer. Given that both sets of boundary conditions yield unstable solutions, we now seek an alternative set of boundary conditions that ensures a negative growth rate, thus guaranteeing stability.
In fact, setting
$\breve S={\textrm {i}}\lambda \breve Y-{\textrm {i}} \breve c_{1}$
represents a more generic form. For BC1, one may assume that
$\breve c_{1}$
appears with a small prefactor. Guided by this asymptotic structure, we aim to construct a set of conformation boundary conditions that ensures
$\breve c_1=O(1)$
and yields a negative imaginary part,
$\breve c_{1i}\lt 0$
. Considering that the wall sublayer is not expected to generate active perturbations, we expand the phase speed according to the bulk-sublayer property, where the second-order correction scales as
$O(\breve \sigma ^{-1})$
. Therefore, the phase speed expansion differs from the earlier expansions (3.40) and (3.57), expressed as
\begin{equation} \breve c_1=\breve c_{10}+\breve \sigma ^{-1}\breve c_{11}+\cdots . \end{equation}
Correspondingly, the perturbation in the bulk region should be expanded as
\begin{equation} \breve V_0=\breve V_{00}+\breve \sigma ^{-1}\breve V_{01}+\cdots . \end{equation}
Comparing with the expansions for the BC12 and BC3, we note that the terms at
$O(\breve \sigma ^{-1/4},\breve \sigma ^{-1/2},\breve \sigma ^{-3/4})$
become trivial in the solution to be constructed. The governing equations are
\begin{align} {\mathcal L}_1\breve V_{00}=0,\quad {\mathcal L}_1\breve V_{01}=-{\textrm {i}}\breve c_{11}{\mathcal D}\breve V_{00}, \end{align}
where
${\mathcal L}_1:=[{\mathcal D}^3-{\textrm {i}}(\lambda \breve Y-\breve c_{10}){\mathcal D}+{\textrm {i}}\lambda ]$
.
The leading-order solution takes the same form as (3.43), with the definition of
$\tilde \xi$
changed to
$\tilde \xi =({\textrm {i}}\lambda )^{1/3}\breve Y+\tilde \xi _0$
, where
$\xi _0=-({\textrm {i}}\lambda )^{1/3}\breve c_{10}/\lambda$
. The second-order solution can be expressed as
\begin{equation} \breve V_{01}= \frac {-{\textrm {i}}\breve c_{11}}{({\textrm {i}}\lambda )^{2/3}}\left (d_{10}+d_{20}\int _\infty ^{\tilde \xi }\mbox{Ai}(\tilde \xi ){\rm d}\tilde \xi \right ) +d_{11}\tilde \xi +d_{21}\left (\tilde \xi \int _\infty ^{\tilde \xi }\mbox{Ai}(\tilde \xi ){\rm d}\tilde \xi -\mbox{Ai}^{\prime}(\tilde \xi )\right )\!. \end{equation}
The lower-limit asymptotic behaviours defined in (3.44) are expressed as
\begin{align}&\qquad\qquad\qquad\quad {\mathcal V}_{00}=d_{10}\tilde \xi _0+d_{20}\left (\tilde \xi _0\int _\infty ^{\tilde \xi _0}\mbox{Ai}(\tilde \xi ){\rm d}\xi -\mbox{Ai}^{\prime}(\tilde \xi _0)\right )\!, \end{align}
\begin{align}&\qquad\quad {\mathcal V}_{01}=\left (d_{10}+d_{20}\int _\infty ^{\tilde \xi _0}\mbox{Ai}(\tilde \xi ){\rm d}\xi \right )({\textrm {i}}\lambda )^{1/3},\quad {\mathcal V}_{02}=d_{20}\mbox{Ai}(\xi _0)({\textrm {i}}\lambda )^{2/3}, \end{align}
\begin{align}& {\mathcal V}_{10}=\frac {-{\textrm {i}}\breve c_{11}}{({\textrm {i}}\lambda )^{2/3}}\left (d_{10}+d_{20}\int _\infty ^{\tilde \xi _0}\mbox{Ai}(\tilde \xi ){\rm d}\tilde \xi \right )+ d_{11}\tilde \xi _0+d_{21}\left (\tilde \xi _0\int _\infty ^{\tilde \xi _0}\mbox{Ai}(\tilde \xi ){\rm d}\xi -\mbox{Ai}^{\prime}(\tilde \xi _0)\right )\!, \end{align}
\begin{align}&\qquad\qquad\quad {\mathcal V}_{11}=\frac {-{\textrm {i}}\breve c_{11}d_{20}}{({\textrm {i}}\lambda )^{1/3}}\mbox{Ai}(\tilde \xi _0)+ ({\textrm {i}}\lambda )^{1/3}\left ( d_{11}+d_{21}\int _\infty ^{\tilde \xi _0}\mbox{Ai}(\tilde \xi ){\rm d}\xi \right )\!, \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\quad {\mathcal V}_{12}={-{\textrm {i}}\breve c_{11}d_{20}}\mbox{Ai}(\tilde \xi _0)+ d_{21}\mbox{Ai}(\tilde \xi _0)({\textrm {i}}\lambda )^{2/3}. \end{align}
Focusing now on the wall sublayer, where
$\eta =O(1)$
, we adopt the same expansion for the solution
$\breve V_0$
as in (3.46). The leading-order solution remains identical in form to (3.50), while the second-order solution takes the form
$\breve v_1={\mathcal V}_{01}\eta +a_1{\textrm {e}}^{\bar \lambda _1\eta }+b_1{\textrm {e}}^{\bar \lambda _2\eta }$
. Our strategy is to minimise the perturbation in the wall sublayer. Being similar to the case for BC3, if we take
\begin{equation} d_{10}\tilde \xi _0+d_{20}\left (\tilde \xi _0\int _\infty ^{\tilde \xi _0}\mbox{Ai}(\tilde \xi ){\rm d}\xi -\mbox{Ai}^{\prime}(\tilde \xi _0)\right )=0, \end{equation}
then the leading-order perturbation
$\breve v_0$
vanishes. Next, we further seek to reduce the magnitude of
$\breve v_{1}$
. Notably, if we take
${\mathcal V}_{01}=0$
and
${\mathcal V}_{10}=0$
, then the solution is simplified to
$\breve v_1=0$
. These conditions can indeed be satisfied, leading to the following required condition:
\begin{align} \left (d_{10}+d_{20}\int _\infty ^{\tilde \xi _0}\mbox{Ai}(\tilde \xi ){\rm d}\xi \right )({\textrm {i}}\lambda )^{1/3}=0. \end{align}
Combining conditions (4.6) and (4.7), we can directly obtain the dispersion relation
\begin{equation} \mbox{Ai}^{\prime}(\tilde \xi _0)=0. \end{equation}
The derivative of the Airy function Ai
$'$
has zero values only in the negative real axis, which corresponds to a negative value of
$\tilde c_{1i}$
, indicating a stable discrete mode. The first zero point of Ai
$'(\tilde \xi _0)$
appears at
$\tilde \xi _0=-1.0188$
. Thus, for
$\lambda =2$
(the PPF base flow), the numerical solution of the rescaled phase speed is
\begin{equation} \breve c_{10}=1.40-0.809{\textrm {i}}, \end{equation}
with the relation between
$d_{10}$
and
$d_{20}$
expressed as
$ d_{20}=d_{10}\int _{\tilde \xi _0}^\infty \mbox{Ai}(\tilde \xi ){\rm d}\tilde \xi$
.
Solving the third- and fourth-order equations in the wall sublayer, and applying the matching condition and the no-slip conditions, we derive
\begin{equation} \breve v_2={\mathcal V}_{02}\eta ^2/2,\quad \breve v_3=0. \end{equation}
The fifth- and sixth-order governing equations are
\begin{align} \big(\bar {\mathcal D}^4+2\lambda ^2\big)\bar {\mathcal D}^3\breve v_4=-{\textrm {i}}\breve c_{10}2\lambda ^2{\mathcal V}_{02}\eta ,\quad \big(\bar {\mathcal D}^4+2\lambda ^2\big)\bar {\mathcal D}^3\breve v_5={\textrm {i}}\lambda ^3{\mathcal V}_{02}\eta ^2. \end{align}
The solutions are
\begin{equation} \breve v_4={\mathcal V}_{10}-\frac {{\textrm {i}}\breve c_{10}{\mathcal V}_{02}\eta ^4}{24}+a_4{\textrm {e}}^{\bar \lambda _1\eta }+b_4{\textrm {e}}^{\bar \lambda _2\eta },\quad \breve v_5={\mathcal V}_{11}\eta +\frac {{\textrm {i}}\lambda {\mathcal V}_{02}\eta ^5}{120}+a_5{\textrm {e}}^{\bar \lambda _1\eta }+b_5{\textrm {e}}^{\bar \lambda _2\eta }. \end{equation}
Applying the no-slip and non-penetration conditions at the wall, we obtain
\begin{align} {\mathcal V}_{10}+a_4+b_4=0,\quad \bar \lambda _1a_4+\bar \lambda _2b_4=0, \end{align}
\begin{align} a_5+b_5=0,\quad \bar \lambda _1a_5+\bar \lambda _2b_5+{\mathcal V}_{11}=0. \end{align}
The seventh-order solution reads
\begin{equation} \breve v_6=\frac {{\mathcal V}_{12}\eta ^2}{2}-\frac {{\textrm {i}}\breve c_{10}{\mathcal V}_{04}\eta ^6}{6!}-\frac {{\textrm {i}}\breve c_{10}\eta }{4}\left (\frac {a_4{\textrm {e}}^{\bar \lambda _1\eta }}{\lambda _1}+\frac {b_4{\textrm {e}}^{\bar \lambda _2\eta }}{\lambda _2}\right )+a_6{\textrm {e}}^{\bar \lambda _1\eta }+b_6{\textrm {e}}^{\bar \lambda _2\eta }. \end{equation}
The no-slip, non-penetration conditions determine
\begin{align} a_6+b_6=0,\quad \bar \lambda _1a_6+\bar \lambda _2b_6-\frac {{\textrm {i}}\breve c_{10}}{4}\left (\frac {a_4}{\bar \lambda _1}+\frac {b_4}{\bar \lambda _2}\right )=0. \end{align}
The next task is to find boundary conditions from high-order derivatives of
$\breve v_{4,5,6}$
. A few representative quantities are listed as follows:
\begin{align}&\qquad\qquad\,\, \breve v_4^{(4)}(0) =2\lambda ^2{\mathcal V}_{10}-{\textrm {i}}\breve c_{10}{\mathcal V}_{02},\quad \breve v_5^{(4)}(0) =-2\lambda ^2\left (a_5+b_5\right )=0, \end{align}
\begin{align}&\qquad\qquad \breve v_4^{(5)}(0) =-2\lambda ^2(\bar \lambda _1a_4+\bar \lambda _2b_4)=0,\quad \breve v_5^{(5)}(0) =2\lambda ^2{\mathcal V}_{11} +{\textrm {i}}\lambda {\mathcal V}_{02}, \end{align}
\begin{align}& \breve v_4^{(3)}(0)=-2\lambda ^2\left (\frac {a_4}{\bar \lambda _1}+\frac {b_4}{\bar \lambda _2}\right )\!,\quad \breve v_6^{(5)}(0) =-2\lambda ^2(\bar \lambda _1a_6+\bar \lambda _2b_6) +2\lambda ^2\frac {5{\textrm {i}}\breve c_{10}}{4}\left (\frac {a_4}{\bar \lambda _1}+\frac {b_4}{\bar \lambda _2}\right )\!. \end{align}
By adjusting
${\mathcal V}_{10}$
and
${\mathcal V}_{11}$
, it is possible to obtain
$\breve v_{4}^{(4)}(0)=\breve v_{5}^{(5)}(0)=0$
. Remarkably, performing
${\textrm {i}}\breve c_{10}\times$
(4.16e
) + (4.16f
), we can derive
\begin{equation} \breve v_6^{(5)}(0)+{\textrm {i}}\breve c_{10}\breve v_4^{(3)}(0)=0. \end{equation}
Combining with
$\breve v_l^{(5)}(0)=\breve v_l^{(4)}(0)=0$
for
$l\in [0,5]$
, we may construct the boundary condition to the perturbation
$\breve V_0$
\begin{equation} \breve{V}_{0}^{(4)}(0)=0,\quad \breve{V}_{0}^{(5)}(0)+{\textrm{i}}\breve{c}_{10}\breve{V}_{0}^{\prime\prime\prime}(0)=0. \end{equation}
To convert (4.18) to the conformation boundary conditions for system (3.29), we expand
$\breve C_{110}$
and
$\breve C_{120}$
in the near-wall region. This expansion takes the same form as in (3.35), but with a change in the coordinate from
$\breve Y$
to
$\eta$
. Substituting into (3.29) leads to
\begin{align} {\textrm {i}}& \breve V_0^{(3)}+\breve \sigma ({\textrm {i}} \bar a_0+\bar b_1)=0,\quad \bar a_1+\bar b_2=0,\quad {\textrm {i}} \breve V_0^{(5)}+\breve \sigma ({\textrm {i}} \bar a_2+\bar b_3)=0, \end{align}
\begin{align}&\quad \bar a_2=-{\textrm {i}}\breve c_{10} \bar a_0-2\lambda \bar b_0,\quad \bar b_2=-{\textrm {i}}\breve c_{10} \bar b_0,\quad \bar b_3=-{\textrm {i}}\breve c_{10} \bar b_1+{\textrm {i}}\lambda \bar b_0. \end{align}
Using the relation (4.18), we can derive
\begin{equation} \bar b_0=\bar b_2=\bar a_1=0, \end{equation}
indicating
$\breve C_{110}^{\prime}(0)=\breve C_{120}(0)=\breve C_{120}^{\prime\prime}(0)=0$
. Therefore, the new conformation boundary conditions can be expressed as
\begin{equation} \breve C_{110}^{\prime}(0)=\breve C_{120}(0)=0, \end{equation}
which is precisely the focus of this paper.
4.1.2. Discussion and validation
The physical interpretation of boundary conditions (4.21) is that, at the wall, the streamwise stretching of the polymer
$\hat c_{11}$
remains with zero gradient in the wall-normal direction, while the shear deformation
$\hat c_{12}$
is completely suppressed. This is consistent with the analysis of polymer conformations in dilute solutions reported in Mavrantzas & Beris (Reference Mavrantzas and Beris1999), which shows that the polymer chains become flat near the wall, with negligible dimension in the perpendicular direction. On the contrary, the BC12 conditions impose an abrupt change of the conformation diffusion (from zero at the wall to a finite value in the internal flow region, as analysed in § 3.5.1), which drives additional instability mechanisms. Similarly, for BC3, even though the conformation tensor remains continuous in the near-wall region, the mismatch of the perturbations in the wall sublayer and those in the bulk sublayer also leads to instability, as analysed in § 3.5.2. In fact, from the physical point of view, we are unable to determine which set of conformation boundary conditions is more appropriate. However, as demonstrated by the asymptotic analysis in § 4.1.1, the boundary conditions (4.21) minimise the perturbation in the wall sublayer, avoiding additional instability induced from the wall, which is favourable in numerical simulations.
In § 2.7 of Lewy & Kerswell (Reference Lewy and Kerswell2024), the mechanism of PDI was discussed based on the leading-order terms derived from the asymptotic analysis in the limit of
$\beta \to 0$
. Through an energy estimate, the growth rate was expressed as (equation (28) in that paper)
\begin{equation} kc_i\overline {|\hat Q|^2}=\cdots +\frac {1}{2}\left [\epsilon \partial _y|\hat Q|^2+\frac {\beta k^2T_{11}}{1-\beta }\partial _y|\hat \omega |^2\right ]_{-1}^1, \end{equation}
where
$\hat Q=[-{\textrm {i}} k \partial _y \hat c_{11}-(\partial _y^2+k^2)\hat c_{12}]/Wi$
,
$\hat \omega ={\textrm {i}} k\hat v-\partial _y \hat u$
,
$T_{11}=(C_{11}-1)/Wi$
and the terms denoted by ’
$\cdots$
’ indicate the volume terms with negative signs. To achieve an unstable PDI with a positive
$c_i$
, the contributions from the boundary must also be positive. If
$\epsilon$
is set to be zero (indicating no ACD is introduced), the first term of the boundary contribution vanishes. While the second term may still remain positive in this situation, the negative volume terms could overwhelm its contribution, resulting in a negative growth rate.
Although the concentration levels assumed in Lewy & Kerswell (Reference Lewy and Kerswell2024) do not align with those in the regime-III analysis, the growth rate estimate in (4.22) still remains an insightful observation, because the dominant contributors in the instability equations remain the same. Based on the near-wall behaviour described in (3.35) with (4.20), we can estimate that
$\hat Q(0)\sim y$
in the near-wall region, indicating that
$\partial _y|\hat Q|^2=0$
at
$y=\pm 1$
. The implication is that the perturbation field governed by the boundary conditions (4.21) ensures that the contribution of the first term on the right-hand side of (4.22) is zero, thereby highlighting the effectiveness of this set of boundary conditions in eliminating unstable PDI.
4.1.3. Application of (4.21) to regime III
To confirm the effectiveness of the boundary conditions (4.21) for finite values of
$\breve \sigma$
, we perform numerical calculations of the regime-III system over the entire range of
$\breve \sigma$
, as shown in figure 9. The black dot-dashed lines demonstrate the rescaled phase speed and growth rate obtained by solving (3.29) with the new boundary conditions (4.21). Notably, these results do not vary with
$\breve \sigma$
, and the values agree precisely with the analytical prediction (4.9). For comparison, we also present the curves obtained based on BC12 and BC3. It is confirmed that our new boundary conditions effectively eliminate the unstable PDI.
Variation with
$\breve \sigma$
of the rescaled phase speed
$\breve c_{1r}$
(a) and growth rate
$\breve c_{1i}$
(b) for different boundary conditions, where
$\lambda =2$
(PPF). The circles in (b) mark the neutral points for BC12 and BC3.
4.2. Effectiveness of the new boundary conditions in regimes I and II
For regime II, the new boundary conditions can be readily translated to
\begin{equation} \tilde c_{110}^{\prime}(0)=\tilde c_{120}(0)=\tilde c_{220}(0)=0. \end{equation}
Therefore, the new eigenvalue system is represented by (3.14) with boundary conditions given by (3.24), (3.15a,b ) and (4.23). Given the existence of multiple eigenvalue solutions within this eigenvalue system, we identify the instability branch through a continuation method, which guarantees that all solutions originate from the same branch. To achieve this goal, we express the wall boundary conditions for the conformation tensor in terms of
\begin{equation} \tilde c_{110}^{\prime}(0)=0, \quad (1-\kappa )\left (\tilde c_{120}(0),\tilde c_{220}(0)\right )=\kappa \left (\tilde c_{120}^{\prime}(0),\tilde c_{220}^{\prime}(0)\right )\!, \end{equation}
where
$\kappa =1$
corresponds to the boundary conditions (3.31) (denoted by BC3) and
$\kappa =0$
indicates the newly developed boundary conditions (4.23). Gradually reducing
$\kappa$
from 1 to 0, we can trace the eigenvalue solutions along the same branch, as illustrated in figure 10(a). For all selected
$\beta$
values, the rescaled growth rate decreases with reducing
$\kappa$
, and for
$\kappa =0$
, all the PDI modes become stable. This figure shows results for one specific value of
$\tilde k$
, while figure 10(b) further explore the rescaled growth rate obtained by boundary conditions (4.23) (or (4.24) with
$\kappa =0$
) across a wide range of
$\tilde k$
. It is confirmed that the unstable PDI is effectively suppressed by the new boundary conditions throughout the entire parameter space, and the eigensolutions are more stable for small
$\tilde k$
values.
(a) Variation of the rescaled growth rate
$\tilde c_{1i}$
with
$\kappa$
for boundary conditions (4.24), where
$\tilde k=177$
and
$\lambda =2$
(PPF). The horizontal dashed line represent the zero growth line. (b) Variation of
$\breve c_{1i}$
with
$\tilde k$
for boundary conditions (4.23). The same
$\beta$
values as in figure 5 are selected.
For regime I, the eigenvalue system is described by (3.5), along with boundary conditions (3.6), (3.7a , b ) and the conformation boundary conditions
\begin{equation} \hat c_{110}^{\prime}(0)=\hat c_{120}(0)=\hat c_{220}(0)=0. \end{equation}
Applying the same continuation method as for regime II, we obtain the rescaled growth rate
$c_{1i}$
for various
$\textit{Wi}$
,
$\bar k$
and
$\beta$
values, as shown in figure 11. Once again, the eigenvalue calculations for the PPFs indicate that all modes become stable.
The rescaled growth rate
$c_{1i}$
for regime I with the new boundary conditions for various
$\beta$
values, where
$\lambda =2$
(PPF). (a) Dependence on
$\textit{Wi}$
with
$\bar k=0.00482$
; (b) dependence on
$\bar k$
with
$\textit{Wi}=40$
.
4.3. Application of the new boundary conditions for system (2.8) with Oldroyd-B model
The validations of our new conformation boundary conditions, as outlined above, are conducted based on the asymptotic systems that consider the asymptotic structures of PDI. To confirm the robustness of this set of boundary conditions, it is essential to also examine its performance within the original system (2.8), to rule out the emergence of new types of instabilities. At this stage, the eigenvalue system is represented by (2.8) with boundary conditions given by (2.10) and
\begin{align} \hat c_{11}^{\prime}(\pm 1)=0,\quad \hat c_{12}(\pm 1)=\hat c_{22}(\pm 1)=\hat c_{33}(\pm 1)=0. \end{align}
To obtain a broad range of the eigenspectra, the SC method is employed. In this subsection, the base flow is selected as the steady solution of the equation system (2.3) with
$\epsilon$
set to be zero; the impact of the ACD on the base flow will be discussed in § 4.5.
Given the five independent controlling parameters of the instability system, demonstrating the effectiveness of the boundary conditions (4.26) across the entire parameter space is not possible. Choosing two representative parameter sets, we compare the eigenspectra obtained using both boundary conditions of BC1 (2.11) and our boundary conditions (4.26), as illustrated by the circles and stars in figure 12, respectively. The parameters in panel (a) correspond to those in figure 1(b) of Couchman et al. (Reference Couchman, Beneitez, Page and Kerswell2024), showing good alignment in the results. Both calculations reveal a branch of discrete modes with small
$c_r$
values, marked by PDI in the figure. Notably, the least stable PDI under BC1 demonstrates a neutral nature, with
$c=c_r+{\textrm {i}} c_i=-0.00139+0{\textrm {i}}$
. When the wall conformation boundary conditions are switched from BC1 to (4.26), the PDI growth rates reduce, ensuring all the discrete modes becoming stable. For a higher
$\textit{Wi}$
value, the eigenvalue calculation based on BC1 reveals an unstable discrete PDI mode with
$c=0.00168+0.00744{\textrm {i}}$
, as illustrated in panel (b). Again, our new boundary conditions effectively eliminate this unstable PDI.
Comparison of the Oldroyd-B PPF eigenspectra obtained based on the boundary condition BC1 (2.11) and the present boundary condition (4.26), where 301 collection points are employed. Panels show (a)
${\textit{Re}}=1000$
,
$\textit{Wi}=22.1$
,
$k=52$
,
$\beta =0.9$
and
$\epsilon =10^{-5}$
; (b)
${\textit{Re}}=1$
,
$\textit{Wi}=2000$
,
$k=56$
,
$\beta =0.9$
and
$\epsilon =10^{-5}$
. The results labelled by ’BC1’ and ’Present’ are derived from boundary conditions (2.11) and (4.26), respectively.
Eigenspectra of system (2.8) for an Oldroyd-B PPF with various
$k$
, with
${\textit{Re}}=1$
,
$\textit{Wi}=100$
,
$\beta =0.9$
and
$\epsilon =10^{-5}$
. For all the calculations, 301 collocation points are employed.
Additionally, we have also selected
$\epsilon =10^{-5}$
,
$\beta =0.9$
,
$\textit{Wi}=100$
and
${\textit{Re}}=1$
for a PPF, and examined the eigenspectra for the Oldroyd-B instability system with our present boundary conditions by adjusting
$k$
across values of 0.1, 0.5, 1.0, 5.0, 20, 50, 100 and 200, as shown by the stars in figure 13. It has been confirmed that no unstable PDI emerges in all these calculations. Although the validation of the effectiveness of our present boundary conditions is not exhaustive across the entire parameter space, the selected cases provide a representative span across perturbations with all possible length scales. Hence, we are confident in the efficacy of our conformation boundary conditions in eliminating unstable PDI.
4.4. Application of the new boundary conditions for system (2.8) with FENE-P model
The FENE-P model is another widely used constitutive model for simulating viscoelastic fluids. In this model, the relation between the conformation stress and tensor (2.4) is modified by
\begin{equation} \boldsymbol \tau =\frac {f(tr\boldsymbol c)\boldsymbol c-\boldsymbol I}{\textit{Wi}},\quad f(s)=\left (1-\frac {s}{L^2}\right )^{-1}, \end{equation}
where
$L$
denotes the maximum extensibility of the polymer chains. In the limit of
$L\to \infty$
, the FENE-P model converges to the Oldroyd-B model. Since the asymptotic structures of the instabilities for both constitutive models are identical, the PDI mode also appears in the FENE-P fluids (Couchman et al. Reference Couchman, Beneitez, Page and Kerswell2024) when the ACD is introduced. The instability system for the FENE-P fluids is described by (2.8) with
\begin{equation} (\hat \tau _{11},\hat \tau _{12},\hat \tau _{22},\hat \tau _{33})=\frac {f(\textit{tr} \boldsymbol{C})}{\textit{Wi}}(\hat c_{11},\hat c_{12},\hat c_{22},\hat c_{33}). \end{equation}
At each wall, in addition to the no-slip, non-penetration conditions (2.10) for velocity perturbations, we also require the boundary conditions for the conformation tensor. For BC1, the boundary conditions are
\begin{align} \big[({-{\textrm {i}} k c+\textit{Wi}^{-1}f(\textit{tr}\boldsymbol{C})})\hat c_{11}-{2(C_{12}\hat u'+U'\hat c_{12})}\big]_{y\,=\,\pm 1}=0, \end{align}
\begin{align} \big[({-{\textrm {i}} k c+\textit{Wi}^{-1}f(\textit{tr}\boldsymbol{C})})\hat c_{12}-{\hat u'}\big]_{y\,=\,\pm 1}=0,\quad \hat c_{22}(\pm 1)=0. \\[9pt] \nonumber \end{align}
In contrast, to eliminate the PDI, it is also necessary to apply our boundary conditions (4.26).
In figure 14, we study the effect of our new boundary conditions on suppressing PDI in FENE-P fluids for two representative parameter sets with different
$\textit{Wi}$
values. The base flow is given by the analytical solution (2.14) of Zhang et al. (Reference Zhang, Lashgari, Zaki and Brandt2013). The blue circles denote the calculations based on the traditional BC1 (4.29), which show a branch of discrete PDI mode at small
$c_r$
values. In panel (a), the least stable mode is neutral with
$c=-0.00394+0{\textrm {i}}$
, while in panel (b), it is unstable with
$c=-0.00315+0.000386{\textrm {i}}$
. However, when the conformation boundary conditions are replaced by (4.26), as shown by the red stars, the least stable PDI mode becomes stable. We have also examined the eigenvalue spectra for various values of
$\textit{Wi}$
,
${\textit{Re}}$
and
$\beta$
under our present boundary conditions, and similar results were observed, although they are not shown here for brevity. Thus, we conclude that our present boundary conditions effectively eliminate the unstable PDI for both the Oldroyd-B and FENE-P models.
Comparison of the eigenspectra of FENE-P PPF instability system between BC1 (blue circles) and our boundary conditions (red stars), where
${\textit{Re}}=1000$
,
$k=65$
,
$\beta =0.7$
,
$\epsilon =10^{-5}$
and
$L=200$
. The number of collocation points is
$N=301$
. Panels show (a)
$\textit{Wi}=7.83$
; (b)
$\textit{Wi}=40$
.
4.5. Application in DNS
In the above calculations, the base flow was obtained by solving the equations without ACD, which can be expressed analytically. This choice is convenient for stability analysis and does not introduce significant errors, as the artificial diffusion is considered to be minimal. However, for DNSs, we must account for the impact of ACD on the base flow, since both the base flow and the perturbations are governed by the same system of equations and boundary conditions. Given that the perturbations
$(\hat c_{12},\hat c_{22},\hat c_{33})$
vanish at the wall, as indicated by (4.26b
,
c
,
d
), we can establish the boundary conditions for
$(C_{12},C_{22},C_{33})$
to be close to the base-flow state in (2.7). A straightforward choice is to set the boundary conditions for
$(C_{12},C_{22},C_{33})$
to be the same as those obtained without ACD. For instance, using the Oldroyd-B constitutive model, we can simply set
$(C_{12},C_{22},C_{33})=(Wi U',1,1)$
at
$y=\pm 1$
. On the other hand, because the boundary condition for
$\hat c_{11}$
, as presented in (4.26a
), shows a zero gradient, we need to establish the boundary condition of
$C_{11}$
in Neumann form. A reasonable choice for this boundary condition is
$C_{11}^{\prime}(\pm 1)=0$
. Although this does not satisfy the base flow (2.7a
) for a PPF, its impact on the base-flow velocity is minimal. In summary, we choose the wall boundary conditions for the Oldroyd-B mean conformation tensor as
\begin{align} C_{11}^{\prime}(\pm 1)=0,\quad C_{12}(\pm 1)=\textit{WiU}'(\pm 1),\quad C_{22}(\pm 1)=C_{33}(\pm 1)=1. \end{align}
For the Oldroyd-B model, we compare the base flow obtained from the ACD equations using our boundary conditions (4.30) with the analytical solutions given in (2.6) and (2.7). We find that a quantitative discrepancy arises only in the conformation component
$C_{11}$
, especially in the near-wall region, as displayed in figure 15(a). However, this discrepancy does not lead to noticeable differences in the eigenspectra calculations, as shown in panel (b). This indicates that the application of the conformation boundary conditions for the base flow (4.30) is indeed appropriate. The same test was also performed for the FENE-P model, yielding a similar conclusion.
(a) Comparison of
$C_{11}$
of a PPF obtained from the equations without ACD versus those with ACD subject to boundary conditions (4.30). (b) Comparison of the eigenspectra based on the base flows in (a), where the perturbation boundary conditions (4.26) are applied for both calculations. The controlling parameters are
${\textit{Re}}=1$
,
$\textit{Wi}=2000$
,
$k=56$
,
$\beta =0.9$
and
$\epsilon =10^{-5}$
.
In DNSs, we typically have two options for selecting the unknowns: (i) the instantaneous flow field
$\boldsymbol \phi$
, and (ii) the perturbation field
$\hat {\boldsymbol \phi }=\boldsymbol \phi -\boldsymbol \varPhi$
. For the latter, the base flow must be obtained analytically or numerically in advance, and the boundary conditions (4.26) are applied to the perturbation equations. For the former, we denote the base-flow conformation tensor without ACD by
$\bar {\boldsymbol C}$
, and the boundary conditions for the total conformation tensor are
\begin{align} \frac {\partial c_{11}}{\partial y}=0,\quad \big(c_{12},c_{22},c_{33}\big)=(\bar C_{12},\bar C_{22},\bar C_{33})\quad \mbox{at }y=\pm 1. \end{align}
For instance, in the Oldroyd-B model, this set of boundary conditions is expressed as
\begin{align} \frac {\partial c_{11}}{\partial y}=0,\quad (c_{12},c_{22},c_{33})=\left(Wi\frac {\partial U}{\partial y},1,1\right)\quad \mbox{at }y=\pm 1. \end{align}
For three-dimensional cases, we can simply set
$c_{13}=\bar C_{13}$
and
$c_{23}=\bar C_{23}$
at the wall.
We now perform DNSs to investigate the evolution of perturbations in Oldroyd-B polymeric fluids initialised with minimal random perturbations. In our DNS, the governing equations (2.3) are discretised using a spectral method in the streamwise and spanwise directions, and a finite-difference method based on collocated Chebyshev–Lobatto notes in the wall-normal direction. A second-order predictor–corrector scheme based on the Crank–Nicolson method is employed for time advancing. Periodic boundary conditions are imposed in the streamwise and spanwise directions, while the no-slip, non-penetration conditions, combined with either the conventional BC1 or our new conformation conditions are employed at the wall. Specifically, in the present two-dimensional simulations, we choose the domain length to be
$L_x=2\pi$
and
$L_y=2$
, with
$N_k=64$
Fourier modes in the streamwise direction and
$N_y=400$
Chebyshev–Lobatto nodes in the wall-normal direction. The time step is set as
$\varDelta _t=0.0001$
. The DNS solver was previously validated and used in Sun, Wan & Zhang (Reference Sun, Wan and Zhang2021) for simulating nonlinear perturbation evolution in viscoelastic pipe flows.
Figure 16(a) displays the temporal evolution of the perturbation kinetic energy
$E_k$
at a specific wavenumber,
$k=56$
, using boundary conditions BC1, which reveals a clear exponential growth within the time interval
$t\in [25,110]$
. The growth rate closely resembles the LST prediction,
$\omega _i=kc_i=0.4166$
, and the perturbation contours displayed in the inset highlight the PDI perturbation in the near-wall region, characterised by a small wavelength. In contrast, the perturbation evolution using our present boundary conditions (4.31) does not exhibit any evident instability at all over a considerably extended time interval, confirming the effectiveness of our boundary conditions in suppressing the unstable PDI.
Evolution of the perturbation kinetic energy of the Oldroyd-B PPF fluids for
$k=56$
obtained by DNS, where
${\textit{Re}}=1$
,
$\textit{Wi}=2000$
,
$k=56$
,
$\beta =0.9$
and
$\epsilon =10^{-5}$
. (a) Boundary conditions BC1 by setting
$\epsilon =0$
at the wall; (b) new boundary conditions (4.31). The inbox in (a) shows the contours of the perturbation wall-normal velocity of PDI at
$t=100$
.
(a) Eigenspectra for
${\textit{Re}}=190$
,
$\textit{Wi}=76$
,
$\beta =0.9$
,
$k=2$
and
$\epsilon =10^{-6}$
. (b) Eigenspectra for
${\textit{Re}}=8000$
,
$\textit{Wi}=20$
,
$\beta =0.9$
,
$k=1$
and
$\epsilon =10^{-6}$
. The insets in (a) and (b) provide zoomed-in views of the unstable modes. (c,d) Temporal evolution of perturbation kinetic energy
$E_k$
of the most unstable modes in (a,b). The insets in (c,d) present the contours of trace of the perturbation conformation tensor at
$t=5000$
.
We now evaluate whether the new set of conformation boundary conditions affects the calculations of other instability modes, including centre modes (Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018; Dong & Zhang Reference Dong and Zhang2022) and wall modes (Shenkar et al. Reference Shenkar, McMullen, McKeon and Graham2021). Figure 17(a) presents the eigenspectra in the
$c_r$
-
$c_i$
plane for
${\textit{Re}}=190$
,
$\textit{Wi}=76$
,
$\beta =0.9$
,
$k=2$
and
$\epsilon =10^{-6}$
, where an unstable centre mode, with
$c=0.999+0.00134{\textrm {i}}$
, emerges. The results obtained using both BC1 and the new boundary conditions demonstrate excellent agreement. Having chosen these parameters, we carry out DNS based on the new boundary conditions initialised from minimal random perturbations, with the evolution of the perturbation kinetic energy shown in panel (c). After a certain period, the centre mode becomes apparent, with the energy growth
$2\omega _i$
aligning perfectly with the LST prediction. The inset displays the contours of the trace of the perturbation conformation tensor, exhibiting the characteristic arrow shape of centre modes. Similarly, in panel (b), we plot the eigenspectra for the parameters featuring an unstable wall mode, indicated by
$c=0.249+0.00173{\textrm {i}}$
. Using the new boundary conditions, we perform DNS for this parameter set in panel (d), observing that the exponential growth of the wall mode aligns well with the LST prediction. The perturbation field is illustrated in the inset. These results confirm that the new set of conformation boundary conditions does not impact the computations of other instability modes, and thus is expected to be an effective means for simulating transition routes in viscoelastic flows.
5. Concluding remarks
Although ACD is regarded as an effective tool for suppressing numerical instability in simulations of transitional and turbulent polymer-solution flows, the recent discovery of PDI presents a significant challenge to its application in wall-shear flows. In this paper, we aim to tackle this issue by proposing a set of new boundary conditions for the conformation tensor to avoid the emergence of PDI.
To achieve this task, we first employ an asymptotic technique that reveals four distinct regimes across different parameter ranges. Regime I occurs at small conformation diffusivity,
$\epsilon \ll 1$
, with the streamwise wavenumber scaling as
$k\sim \epsilon ^{-1/2}$
. For both PPFs and PCFs, unstable PDI modes emerge when the Weissenberg number
$\textit{Wi}$
exceeds a threshold that is dependent on the concentration level. As
$\textit{Wi}$
increases, the peak of the PDI growth rate is observed at lower
$\bar k=k\epsilon ^{1/2}$
values. This observation prompts us to explore regime II, where
$k\epsilon ^{1/2}\sim \textit{Wi}^{-3/2}$
with
$\textit{Wi}\gg 1$
. In this regime, the near-wall diffusive layer splits into two sublayers, with the dispersion relation determined within the lower diffusive layer. Numerical solutions illustrate the emergence of unstable PDI for
$\tilde k=\textit{Wi}^{3/2}\epsilon ^{1/2}k$
ranging from a critical value to infinity. Furthermore, different choices of conformation boundary conditions result in an inverse dependence of the growth rate on polymer concentration. As a result, further investigation into PDI in the dilute limit, referred to as regime III, becomes of particular interest. This regime specifically focuses on the scaling relation
$k\epsilon ^{1/2}\textit{Wi}^{3/2}\sim (1-\beta )^{-3/2}$
, leading to a reduction in the number of controlling parameters from five (in the original instability system) to just one (
$\breve \sigma$
). To simplify the instability system even further, we subsequently consider the limit
$\breve \sigma \to \infty$
, leading to a parameter-free system (referred to as regime IV) that facilitates an analytical solution of the dispersion relation.
Then, we aim to derive the new conformation boundary conditions based on the simplified instability system inferred from regime IV. Based on the concept that the perturbations in the wall sublayer should be minimised, we construct an analytical solution with a negative growth rate. By examining the near-wall behaviour of this solution, we derive a set of conformation boundary conditions designed for eliminating unstable PDI. For the streamwise stretch component
$\hat c_{11}$
, we impose a Neumann condition
$\hat c_{11}^{\prime}=0$
, so that near the wall
$\hat c_{11}$
does not vary in the wall-normal direction. This prevents any artificial buildup or depletion of polymer stretching at the boundary and lets the bulk flow set its value. For all other conformation components (
$\hat c_{12},\hat c_{22}$
), we impose Dirichlet conditions, which enforces zero shear or normal stretching of the polymer in the wall-normal direction. There conditions are readily carried over to the governing systems in regimes I, II and III.
The effectiveness of our new conformation boundary conditions can also be interpreted through the analytical model for the onset of unstable PDI developed in Lewy & Kerswell (Reference Lewy and Kerswell2024). In this model, the PDI growth rate is primarily influenced by the shear of the conformation tensor and the shear of the vorticity at the boundaries (see (4.22)), in addition to negative volume contributions. Notably, the boundary conformation shear emerges as the dominant factor that ensures a positive PDI growth rate. By applying our new conformation boundary conditions, this term favourably becomes zero, indicating the elimination of any unstable PDI modes.
For representative configurations, this set of conformation boundary conditions has been validated as effective in eliminating unstable PDI within the ACD instability system, using both the Oldroyd-B and FENE-P constitutive models. Therefore, we can expect their effectiveness throughout the entire parameter space. Moreover, the new conformation boundary conditions can be conveniently integrated into the DNS code with ACD, as discussed in § 4.5. Comparing the DNS results based on the original boundary conditions in Beneitez et al. (Reference Beneitez, Page and Kerswell2023) and our present boundary conditions, we confirm that the latter does ensure stable simulations of polymer fluids without unstable PDI. Remarkably, the new set of conformation boundary conditions does not affect the calculations of other types of instabilities, including the centre and wall modes.
It needs to be emphasised that the new set of conformation boundary conditions proposed in this paper is not the unique way to eliminate unstable PDI, and we can expect more types of effective boundary conditions. With this set of boundary conditions ensuring the absence of unstable PDI throughout the parameter space, using ACD in simulations of transitional and turbulent polymer-solution flows can be considered as a reliable and efficient approach for future studies.
Funding
This research is supported by National Natural Science Foundation of China (92371104, 12588201), CAS project for Young Scientists in Basic Research (YSBR-087), and Strategic Priority Research Program of Chinese Academy of Sciences (XDB0620102).
Declaration of interests
The authors report no conflict of interest.
