1. Introduction
Natural convection of yield stress fluids arises in a wide range of industrial processes and technologies, including food canning and sterilisation (Datta & Teixeira Reference Datta and Teixeira1988; Ghani et al. Reference Ghani, Farid, Chen and Richards1999), food processing (Balmforth, Frigaard & Ovarlez Reference Balmforth, Frigaard and Ovarlez2014; Bonn et al. Reference Bonn, Denn, Berthier, Divoux and Manneville2017), pulp processing (Derakhshandeh et al. Reference Derakhshandeh, Kerekes, Hatzikiriakos and Bennington2011), electronics cooling (Bahiraei & Heshmatian Reference Bahiraei and Heshmatian2018), solar thermal collectors (Nagarajan et al. Reference Nagarajan, Subramani, Suyambazhahan and Sathyamurthy2014) and biotechnological processes such as the micro-polymerase chain reaction technique for DNA amplification (Zhang et al. Reference Zhang, Xu, Ma and Zheng2006a ). It also plays a role in geophysical contexts involving yield stress fluids, such as magma flows and convection in the Earth’s mantle (Orowan Reference Orowan1965; Petford Reference Petford2003; Balmforth et al. Reference Balmforth, Frigaard and Ovarlez2014). Given the widespread occurrence of naturally convecting flows, it has been the subject of sustained research interest, not only for Newtonian fluids (Ostrach Reference Ostrach1988), but increasingly for yield stress fluids.
Two canonical flow configurations to study internal natural convection comprise imposing a temperature difference across the horizontal or vertical walls of a square cavity. When the top and bottom boundaries are cooled and heated (respectively), one has the classical Rayleigh--Bénard configuration. When the thermal gradient is perpendicular to the direction of gravity and two thermally insulated boundaries confine the top and bottom of the fluid, we have the configuration of the differentially heated cavity; see, e.g. de Vahl Davis (Reference De Vahl Davis1983). This latter configuration is the focus of the present work.
Both configurations have garnered interest over the past several decades. Zhang, Vola & Frigaard (Reference Zhang, Vola and Frigaard2006b ) presented the first theoretical predictions of the stability of yield stress fluids in the Rayleigh–Bénard set-up. Their analysis showed that, unlike Newtonian fluids, viscoplastic fluids with a finite yield stress remain linearly stable at all Rayleigh numbers and that the Rayleigh–Bénard instability can only be triggered by a finite-amplitude perturbation. The weakly nonlinear stability analysis by Balmforth & Rust (Reference Balmforth and Rust2009) revealed a scaling law for the perturbation amplitude required to initiate Rayleigh–Bénard convection. Turan et al. (Reference Turan, Chakraborty and Poole2012a ) presented a systematic parametric study of the steady-state dynamics using numerical simulations of the regularised Bingham model. They demonstrated flow arrest at a critical yield stress and provided correlations for the average Nusselt number.
Contrary to the theoretical studies of Zhang et al. (Reference Zhang, Vola and Frigaard2006b ) and Balmforth & Rust (Reference Balmforth and Rust2009), subsequent experimental investigations by Darbouli et al. (Reference Darbouli, Metivier, Piau, Magnin and Abdelali2013) and Kebiche, Castelain & Burghelea (Reference Kebiche, Castelain and Burghelea2014) reported the onset of Rayleigh–Bénard convection from a motionless initial state at sufficiently large Rayleigh numbers, without the need to impose external disturbances. More recently, Metivier et al. (Reference Metivier, Brochard, Darbouli and Magnin2020) observed oscillatory sub-yield motion in Rayleigh–Bénard experiments with Carbopol gels, taking the form of periodic travelling waves hypothesised to arise from the material’s elasticity.
Numerical and experimental studies have attempted to reconcile the discrepancies between theoretical and numerical predictions based on viscoplastic models and experimental observations of yield stress fluids by, for example, comparing the necessary conditions for flow arrest (Li, Magnin & Métivier Reference Li, Magnin and Métivier2016), explaining the Rayleigh–Bénard convection of Carbopol gels as that of a Newtonian solvent saturating a porous medium (Metivier, Li & Magnin Reference Metivier, Li and Magnin2017), or examining the influence of boundary conditions (Ahmadi, Olleik & Karimfazli Reference Ahmadi, Olleik and Karimfazli2022).
Nevertheless, the contradiction regarding transient dynamics remains unresolved, underscoring the limitations of viscoplastic models in predicting the onset of motion and flow behaviour in Rayleigh–Bénard convection of yield stress fluids (such as the Carbopol gels) in experimental studies.
In the case of the differentially heated cavity, several numerical studies (Lyubimova Reference Lyubimova1977; Turan, Chakraborty & Poole Reference Turan, Chakraborty and Poole2010, Reference Turan, Poole and Chakraborty2012b ) have investigated the modulation of steady-state dynamics and heat transfer by a finite yield stress, as well as the critical criterion beyond which the flow becomes immobilised (Vikhansky Reference Vikhansky2009, Reference Vikhansky2010, Reference Vikhansky2011; Karimfazli & Frigaard Reference Karimfazli and Frigaard2013; Huilgol & Kefayati Reference Huilgol and Kefayati2015). Karimfazli et al. (Reference Karimfazli, Frigaard and Wachs2015) studied the transient dynamics of natural convection in smart viscoplastic fluids in which the yield stress is switched on and off. This enabled tuning the heat transfer rate between the convective and conductive limits for heat transfer control applications. Later, Karimfazli & Frigaard (Reference Karimfazli and Frigaard2016) derived the sufficient conditions for the onset of natural convection in viscoplastic fluids. Moreover, they showed that viscoplastic fluids may exhibit a finite time delay before flow onset in situations where Newtonian fluids begin to convect immediately.
Recently, Jadhav, Rossi & Karimfazli (Reference Jadhav, Rossi and Karimfazli2021) conducted an experimental study on natural convection of Carbopol gels in a differentially heated cavity. They measured time-resolved flow velocities and temperature through particle image velocimetry (PIV) and thermometry. While at steady state, their flow field measurements revealed classical viscoplastic flow patterns, such as the formation of the central convection cell and quiescent unyielded regions at the cavity corners, they also observed transient dynamics that are in conflict with theoretical predictions. Specifically, contrary to the delayed flow onset predicted by Karimfazli & Frigaard (Reference Karimfazli and Frigaard2016), immediate motion onset was observed at subcritical yield stress values. In these cases, the onset was characterised by a rapid increase in the system’s kinetic energy, which was not anticipated by viscoplastic models (Karimfazli & Frigaard Reference Karimfazli and Frigaard2016). More strikingly, experiments revealed immediate motion onset even at supercritical yield stress values, i.e. conditions under which the steady state is motionless. These apparent discrepancies further demonstrate that viscoplastic models cannot adequately capture the transient dynamics of natural convection in yield stress fluids.
Irrespective of the specific flow configuration, various hypotheses have been proposed to explain the discrepancies between theoretical and numerical predictions of motion onset in viscoplastic fluids and experimental observations of yield stress fluids. However, a coherent understanding has yet to emerge. The present paper aims to help reconcile some of these differences in the context of natural convection of a yield stress fluid in a cavity with differentially heated walls. Specifically, we investigate the influence of elasticity on the onset of material deformation and flow at both subcritical and supercritical yield stress values.
In our analysis, we characterise the onset dynamics under two initial temperature distributions. First, we consider a linear initial temperature field, i.e. the steady pure conduction profile established between two differentially heated walls. This configuration facilitates comparison with the numerical simulations of Karimfazli et al. (Reference Karimfazli, Frigaard and Wachs2015), where viscoplastic fluids are predicted to begin flowing immediately only if the steady state itself is flowing (Karimfazli & Frigaard Reference Karimfazli and Frigaard2016). Second, we explore a uniform initial temperature distribution, which is more representative of experimental set-ups in which the fluid is initially in thermal equilibrium with its surroundings. This enables comparison with the time-resolved flow field measurements of Jadhav et al. (Reference Jadhav, Rossi and Karimfazli2021). In this case, experimental results suggest immediate motion onset (Jadhav et al. Reference Jadhav, Rossi and Karimfazli2021), whereas numerical simulations of viscoplastic fluids predict delayed flow onset (Karimfazli & Frigaard Reference Karimfazli and Frigaard2016).
The remainder of this paper is organised as follows. In § 2 we discuss the problem set-up, its numerical simulation and present validation cases. In § 3 we present our results; this section is divided into two main subsections. In § 3.1 we present our results for natural convection starting with a linear initial temperature distribution, and discuss the comparison with viscoplastic models. A scaling analysis is presented, which shows how elasticity and yield stress affect motion onset. In § 3.2 we present our results for natural convection starting with a uniform initial condition for the temperature field and compare with existing experimental studies. Our analysis primarily focuses on the impact of elasticity on transient dynamics, as this is where we have observed the most significant differences from viscoplastic flow dynamics across the range of parameters investigated. Finally, we summarise and discuss the key findings of this study in § 4. Steady-state results are documented in Appendix A, while in the main text we focus on the effect of elasticity on motion onset and transient dynamics.
2. Problem set-up and numerical methodology
2.1. Set-up and governing equations
We consider a square cavity of side length
$\hat {L}$
, whose left side is maintained at a high temperature
$\hat {T}=\hat {T}_H$
and the right side maintained at a low temperature
$\hat {T}=\hat {T}_C$
, creating a horizontal thermal gradient. The top and bottom walls are thermally insulated, and all walls are no-slip. The direction of gravity,
$\hat {\boldsymbol{g}}$
, acts vertically, perpendicular to the thermal gradient direction. The origin of the coordinate system is placed at the centre of the cavity, with the
$x$
axis parallel to the direction of the thermal gradient (the horizontal direction) and the
$y$
axis parallel to the direction of gravity (the vertical direction). A schematic of the set-up is shown in figure 1.
Schematic of the heated cavity. The origin of the coordinate system is placed at the centre of the enclosure. The left wall is maintained at the high temperature,
$\hat T_H$
, the right wall at the low temperature,
$\hat T_C$
, and the top and bottom walls are thermally insulated. All walls are no-slip, and gravity acts downwards in the
$-\hat e_y$
direction.
The dynamics of the problem are governed by the following dimensionless groups:
\begin{align} {\boldsymbol{Ra} } = \frac {{\hat u}_{\textit{ref}}{\hat l}_{\textit{ref}}}{\hat {\kappa }} &= \frac {\hat {g}\hat {\alpha }_T\Delta \hat {T} \hat {L}^3}{\hat {\nu }\hat {\kappa }}, \end{align}
\begin{align} {\boldsymbol{Pr} } &= \frac {\hat {\nu }}{\hat {\kappa }}, \end{align}
\begin{align} {\textit{Wi}} = \hat {\lambda }\frac {{\hat u}_{\textit{ref}}}{{\hat l}_{\textit{ref}}} &= \frac {\lambda \hat {g}\hat {\alpha }_T\Delta \hat {T}\hat {L}}{\hat {\nu }}, \end{align}
\begin{align} \beta _s &= \frac {\hat {\nu }_s}{\hat {\nu }}, \end{align}
\begin{align} Bn = \frac {\hat {\tau }_y}{{\hat \tau }_{\textit{ref}}} &= \frac {\hat {\tau }_y}{\hat {\rho }_0 \hat {g}\hat {\alpha }_T\Delta \hat {T} \hat {L}}. \\[9pt] \nonumber \end{align}
Here
$\hat \rho _0$
denotes the density at the reference temperature
$\hat {T}_0=(\hat T_c + \hat T_H)/2$
;
$\hat {g}$
is the acceleration due to gravity,
$\hat {\alpha }_T$
the coefficient of thermal expansion,
$\hat \kappa$
the thermal diffusivity,
$\hat \nu _s=\hat \mu _s/\hat {\rho }_0$
the solvent kinematic viscosity,
$\hat \nu _p=\hat {\mu }_p/\hat {\rho }_0$
the plastic kinematic viscosity and
$\hat \nu = \hat \nu _s + \hat \nu _p=\hat {\mu }/\hat {\rho }_0$
the combined solvent and plastic kinematic viscosity;
$\hat \mu _s$
,
$\hat {\mu }_p$
and
$\hat {\mu }$
are the dynamic solvent viscosity, the dynamic plastic viscosity and the combined solvent and plastic viscosity;
$\beta _s$
is a viscosity ratio,
$\hat \lambda = \hat \mu _p /\hat G$
the relaxation time and
$\hat \tau _y$
the yield stress. The dimensional variables are denoted with a hat
$(\, \hat {}\, )$
and dimensionless variables without.
The Rayleigh number,
${\textit{Ra}}$
, represents the ratio of advective to conductive transport of thermal energy, the Prandtl number,
${\textit{Pr}}$
, represents the ratio of momentum to thermal diffusivity, the Weissenberg number,
${\textit{Wi}}$
, represents the ratio of elastic to viscous stresses and the Bingham number,
${\textit{Bn}}$
, represents the ratio of the yield stress to the characteristic buoyancy stress. The ranges explored for each dimensionless number are shown in table 1. We have kept
$\beta _s=0.5$
constant for all results shown in § 3 to ensure numerical stability. We note that the total effective viscosity is
$\hat {\tau }_y/\hat {\dot {\gamma }}+\hat {\mu }_p+\hat {\mu }_s$
. Therefore, for a finite yield stress, the contribution of the solvent viscosity approaches 0.5 only at very large strain rates. In all cases considered here,
$\dot {\gamma } = \mathcal{O}(10^{-3}) \ll 1$
. Therefore, the contribution of the solvent viscosity to the total viscosity is not dominant in the slow dynamics observed during the start-up flows considered here. The effect of smaller
$\beta _s$
is discussed in Appendix C.4.
Definitions of dimensionless numbers and dimensionless number space of
${\textit{Ra}}$
,
${\textit{Pr}}$
,
${\textit{Wi}}$
,
$\beta _s$
and
${\textit{Bn}}$
that is explored in the present study.
A mechanical analogue of the one-dimensional Saramito constitutive equation as a combination of a spring, two dashpots and a rigid element representing the yield stress.
Our problem set-up consists of a two-dimensional planar configuration governed by the two-dimensional mass and momentum conservation equations, coupled with an advection–diffusion equation for thermal energy transport and the Saramito constitutive equation to model the extra stresses in the elastoviscoplastic (EVP) fluid. In dimensionless form, the governing equations are as follows:
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad \boldsymbol{\nabla } \boldsymbol{\cdot } \boldsymbol{u} = 0, \end{align}
\begin{align}&\qquad\qquad\qquad \frac {\textit{Ra}} {\textit{Pr}} \left ( \frac {\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{u}\right ) = - \boldsymbol{\nabla } p + \beta _s \boldsymbol{\nabla }^2 \boldsymbol{u} + \boldsymbol{\nabla } \boldsymbol{\cdot } \boldsymbol{\tau }^e + T \hat {e}_y, \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad\quad {\boldsymbol{Ra} } \left (\frac {\partial T}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla } T \right ) = \boldsymbol{\nabla }^2 T, \end{align}
\begin{align}& \textit{Wi} \left ( \frac {\partial \boldsymbol{\tau }^e}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{\tau }^e - \boldsymbol{\nabla } \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\tau }^e - \boldsymbol{\tau }^e\boldsymbol{\cdot }\boldsymbol{\nabla } \boldsymbol{u}^{\mathrm{T}} \right ) + \max \left (0, \frac {\tau ^e_ {\textit{eff}} - Bn}{\tau ^e_ {\textit{eff}} }\right ) \boldsymbol{\tau }^e = 2(1-\beta _s) \unicode{x1D63F}. \\[7pt] \nonumber \end{align}
Here,
$\boldsymbol{u}$
denotes the velocity and
$T$
the temperature. Furthermore, we have used the Boussinesq approximation to account for buoyancy stresses arising from density variations due to temperature gradients. Water-based Carbopol gels used in natural convection experiments typically contain no more than 0.1 % Carbopol powder, so their density variation with temperature is almost identical to that of water. For temperature differences up to
$30\,^\circ$
C, the corresponding change in water density is about 1 %, thus satisfying the small density variation requirement of the Boussinesq framework and rendering this approximation suitable. Moreover, for a cavity width of approximately
$15\,\text{cm}$
and a Carbopol yield stress of the order of
$0.1\,\text{Pa}$
with a plastic viscosity around
$1\,\text{Pa}\,\text{s}$
, the resulting Rayleigh number is
${\textit{Ra}} \lesssim 10^6$
, which is comparable to the range considered in this study.
The dimensionless variables are defined as follows:
\begin{equation} \begin{aligned} \boldsymbol{x} = \frac {\hat {\boldsymbol{x}}}{\hat L}, \hspace{0.5cm} \boldsymbol{u} & = \frac {\hat {\boldsymbol{u}}\, \hat \nu }{\hat g \hat {\alpha }_T\Delta \hat T\hat L^2} , \hspace {0.5cm} t = \frac {\hat t\, \hat g\hat {\alpha }_T\Delta \hat T\hat L}{\hat \nu } , \\ p = \frac {\hat p}{\hat \rho _0\hat g\hat {\alpha }_T\Delta \hat T\hat L} , \hspace {0.5cm} T & = \frac {\hat T - \hat T_0}{\hat T_H - \hat T_C} , \hspace {0.5cm} \boldsymbol{\tau }^e = \frac {\hat {\boldsymbol{\tau }}^e}{\hat \rho _0\hat g\hat {\alpha }_T\Delta \hat T\hat L}. \end{aligned} \end{equation}
Equation (2.9) is the Saramito constitutive equation for EVP fluids (Saramito Reference Saramito2009), in which
$\boldsymbol{\tau }^e$
is the extra stress tensor,
$\tau ^e_ {\textit{eff}} = \sqrt {({1}/{2}) \boldsymbol{\tau }_d^e : \boldsymbol{\tau }_d^e} = \sqrt {({1}/{2})\tau ^e_{{d}_{ij}} \tau ^e_{{d}_{ij}}}$
is the second invariant of the deviatoric part of the extra stress tensor (
$\boldsymbol{\tau }_d^e = \boldsymbol{\tau }^e - ({1}/{2})\mathrm{tr}\,{\boldsymbol{\tau }^e}$
) and
$\unicode{x1D63F} = ({1}/{2}) (\boldsymbol{\nabla } \boldsymbol{u} + \boldsymbol{\nabla } \boldsymbol{u}^{\mathrm{T}})$
is the strain rate tensor. The Saramito constitutive equation models the EVP material as a Kelvin-Voigt viscoelastic solid prior to yielding, and as an EVP fluid in its yielded state. Though prototypical yield stress fluids like aqueous Carbopol gels show a Herschel--Bulkley-like shear thinning with a power-law index of approximately
$ 0.5$
, we use a power-law index equal to unity and focus on the effects of elasticity on the natural convection flow. A mechanical analogue of the one-dimensional formulation of the Saramito constitutive equation is illustrated in figure 2. In figure 2,
$\hat \mu _s$
denotes the dynamic solvent viscosity,
$\hat \mu _p$
the dynamic plastic viscosity and
$\hat {G}$
the elastic modulus. Here, two time scales capture the transient behaviour of the material: the relaxation time (
$\hat {\lambda }$
), which represents the time scale over which stresses relax following an imposed deformation, and the retardation time (
$\hat {t}_r = \hat {\mu }_s / \hat {G}$
), which characterises the delayed deformation response under an applied stress.
Boundary conditions imposed at the side walls are
\begin{align}&\,\,\,\, \boldsymbol{u}(x = \pm 1/2, y) = \mathrm{\textbf {0}}, \end{align}
\begin{align}&\qquad\quad \frac {\partial p}{\partial n} = 0, \end{align}
\begin{align}& T(x = \pm 1/2, y) = \mp 1/2, \end{align}
\begin{align}&\,\, \frac {\partial \boldsymbol{\tau }^e}{\partial n}(x = \pm 1/2, y) = \mathrm{\textbf {0}}, \end{align}
and the top and bottom walls
\begin{align}&\,\, \boldsymbol{u}(x, y=\pm 1/2) = \mathrm{\textbf {0}}, \end{align}
\begin{align}&\qquad\quad \frac {\partial p}{\partial n} = 0, \end{align}
\begin{align}& \frac {\partial T}{\partial n}(x, y = \pm 1/2) = 0, \end{align}
\begin{align}& \frac {\partial \boldsymbol{\tau }^e}{\partial n} (x, y = \pm 1/2) = \mathrm{\textbf {0}} . \end{align}
For initial conditions, we always start from a quiescent fluid,
$\boldsymbol{u} (t=0) =\boldsymbol{0}$
, and from zero extra stresses,
$\boldsymbol{\tau }^e (t=0) = \boldsymbol{0}$
. Regarding the temperature, we explore the effect of two different initial conditions on the natural convection. In § 3.1 we use
\begin{align} T(x,y,t=0) = -x, \end{align}
i.e. the initial temperature distribution linearly decreases from the left (hot) wall to the right (cold) wall, representative of pure conduction between differentially heated walls.
In § 3.2 we use
\begin{align} T(x,y,t=0) = 0, \end{align}
i.e. a uniform temperature distribution, representative of typical experimental conditions.
2.2. Numerical methodology and validation
The domain in figure 1 is discretised with
$N_x\times N_y = 384\times 384$
uniformly distributed grid points. We have performed grid convergence studies with a mesh of
$256\times 256$
and
$512\times 512$
grid points to show that this spatial discretisation sufficiently resolves the dynamics of our problem and our quantities of interest, reported in Appendix C. We have also evaluated our temporal resolution to ensure that all relevant time scales are adequately captured; details are provided in Appendix C.2.
The governing equations are solved using our in-house finite difference code for EVP fluids (Izbassarov et al. Reference Izbassarov, Rosti, Ardekani, Sarabian, Hormozi, Brandt and Tammisola2018). Spatial derivatives are discretised by a second-order central difference scheme, with the exception of the advection term in (2.9), for which we use the classical fifth-order weighted essentially non-oscillatory method (Shu Reference Shu2009). A low-storage, third-order Runge–Kutta scheme (Wesseling Reference Wesseling2009) is used to temporally integrate all evolution equations. The projection method (Chorin Reference Chorin1968) is used to impose the incompressibility constraint on
$\boldsymbol{u}$
. A fast Fourier transform-based fast Poisson solver is employed to solve the pressure-Poisson equation. For homogeneous Neumann boundary conditions, the pressure-Poisson equation is transformed into spectral space using the discrete cosine transform, resulting in a tri-diagonal system of linear equations. This system is efficiently solved via tri-diagonal inversion, and the solution is then transformed back to physical space (Schumann & Sweet Reference Schumann and Sweet1988). Further details of our code can be found in Izbassarov et al. (Reference Izbassarov, Rosti, Ardekani, Sarabian, Hormozi, Brandt and Tammisola2018).
To verify the accuracy of our solver for our set-up, we compare with available results in the literature for natural convection of a (i) Newtonian fluid, (ii) viscoplastic Bingham fluid, and (iii) viscoelastic FENE-P fluid. First, we compare with the classical work of de Vahl Davis (Reference De Vahl Davis1983) and also Karimfazli et al. (Reference Karimfazli, Frigaard and Wachs2015) for natural convection of Newtonian fluids at different
${\textit{Ra}}$
and
$Pr = 0.71$
. The comparison is presented in table 2, showing very good agreement with our solution for the average (
$\overline {\textit{Nu}}$
) and maximum (
$Nu_{\textit{max}}$
) Nusselt numbers, and the maxima of the horizontal (
$u$
) and vertical (
$v$
) velocity components. Second, in figure 3 we present the temporal evolution of the
$L^2$
norm of the velocity magnitude,
$\|{U}\|$
, normalised by its time history maximum for an EVP fluid. Simulations are performed at
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$Wi = 0.1$
and
$\beta _s = 0.15$
(parameters close to the viscoplastic limit) at different
${\textit{Bn}}$
. The results are compared with those of Karimfazli et al. (Reference Karimfazli, Frigaard and Wachs2015) for a viscoplastic Bingham fluid (shown as asterisk markers), showing very good agreement. Furthermore, the figure illustrates that at
$Wi = 0.1$
, the elasticity of an EVP fluid is sufficiently low to have a negligible effect on both the transient and steady results compared with viscoplastic behaviour. Accordingly,
$Wi = 0.1$
is included in various figures as a reference, serving as an approximation of the viscoplastic case. Finally, we simulate natural convection of a viscoelastic fluid using the FENE-P constitutive equation and compare with the results of Chauhan, Sahu & Sasmal (Reference Chauhan, Sahu and Sasmal2021) in figure 4. The simulation is performed at
${\textit{Ra}} = 10^6$
,
$Pr = 7$
,
$Wi = 3779.65$
,
$\beta _s = 0.5$
,
$\mathcal{L}_{\textit{FP}}^2 = 10$
(where
$\mathcal{L}_{\textit{FP}}$
is the FENE-P maximum extensibility). Using a mesh of
$N_x\times N_y = 384\times 384$
, we find excellent agreement with the mid-section velocity profiles and left wall Nusselt number profile, showing that our code and discretisation accurately resolve the viscoelastic stresses that modulate the flow and heat transfer in natural convection.
Comparison of the natural convection of a Newtonian fluid in a square cavity with differentially heated side walls with the results of de Vahl Davis (Reference De Vahl Davis1983) and Karimfazli et al. (Reference Karimfazli, Frigaard and Wachs2015) at different
${\textit{Ra}}$
and
$Pr = 0.71$
.
Evolution of
$\|{U}\|$
, normalised with its time history maximum, for natural convection of an EVP fluid in a square cavity with differentially heated side walls as a function of
${\textit{Bn}}$
. The solid line corresponds to our solution, and the asterisk markers (
$\ast$
) correspond to the results of Karimfazli et al. (Reference Karimfazli, Frigaard and Wachs2015) for a viscoplastic Bingham fluid. The dimensionless numbers are
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$Wi = 0.1$
,
$\beta _s = 0.15$
.
Profiles of
$(a)$
$u$
at
$x = 0$
,
$(b)$
$v$
at
$y = 0$
,
$(c)$
$Nu$
at
$x = -0.5$
for the natural convection of a FENE-P fluid in a square cavity with differentially heated side walls. The solid line shows our solution, and the asterisk markers (
$\ast$
) show the results of Chauhan et al. (Reference Chauhan, Sahu and Sasmal2021). The dimensionless numbers are
${\textit{Ra}} = 10^6$
,
$Pr = 7$
,
$Wi = 3779.65$
,
$\beta _s = 0.5$
,
$\mathcal{L}_{\textit{FP}}^2 = 10$
(
$\mathcal{L}_{\textit{FP}}$
is the FENE-P maximum extensibility). Note that the dimensionless numbers are defined and flow variables are normalised according to the definitions in the present study. According to the definitions in Chauhan et al. (Reference Chauhan, Sahu and Sasmal2021),
${\textit{Ra}} = 10^6$
,
$Pr=7$
,
$Wi = 10$
,
$\beta _s = 0.5$
,
$\mathcal{L}_{\textit{FP}}^2 = 10$
.
3. Results
3.1. Linear initial temperature,
$T(t=0) = -x$
In this section we examine the transient dynamics of natural convection of EVP fluids, starting from an initially quiescent state with a temperature distribution corresponding to steady pure conduction. Karimfazli et al. (Reference Karimfazli, Frigaard and Wachs2015) showed that in this configuration if the yield stress of a viscoplastic material exceeds a critical value (corresponding to a critical Bingham number,
${\textit{Bn}}_{\textit{cr}}=1/32$
), the fluid remains rigid, preserving the initial temperature distribution. Conversely, when the yield stress is below the critical threshold, the material yields immediately. The critical Bingham number,
${\textit{Bn}}_{\textit{cr}}$
, is defined such that
\begin{align} \begin{cases} \|{U}\| (\textit{Bn},\, t\to \infty ) \gt 0 & \text{if } \textit{Bn} \lt {\textit{Bn}}_{\textit{cr}},\\ \|{U}\| (\textit{Bn},\, t\to \infty ) = 0 & \text{if } Bn \geqslant {\textit{Bn}}_{\textit{cr}}; \end{cases} \end{align}
that is,
${\textit{Bn}}_{\textit{cr}}$
is the smallest Bingham number for which the steady state is motionless. Here, we investigate how the transient evolution from this initial temperature distribution is influenced by the presence of elasticity. Our analysis is presented in two parts: we first examine the system dynamics for subcritical
${\textit{Bn}}$
, followed by an exploration of the supercritical regime.
3.1.1.
$T(t=0) = -x$
and
${\textit{Bn}} \lt Bn_{\textit{cr}}$
Figure 5 displays the time evolution of the speed,
$U = \sqrt {u^2 + v^2}$
, at
${\textit{Ra}}=10^4$
and
${\textit{Bn}}=0.02$
for three values of
$Wi = \{0.1, 0.5, 5\}$
. The overlaid arrows represent quiver plots of the velocity vector. Each row corresponds to a different
${\textit{Wi}}$
. The first column shows the time evolution of the
$L^2$
norm of the velocity magnitude with the red markers indicating the time instance of the snapshots on each row -- a movie of this transient motion is provided in supplementary movie 1 available at https://doi.org/10.1017/jfm.2025.10732. Here, the
$L^2$
norm is denoted by
$\|{\boldsymbol{\cdot }}\|$
and defined as
\begin{align} \|{\xi }\| = \sqrt {\int _{-1/2}^{1/2} \int _{-1/2}^{1/2} \left |\xi \right |^2 \, \mathrm{d} x\, \mathrm{d} y}, \end{align}
where
$\xi$
is some observable. In the last three columns, the white dashed lines indicate estimated contours of the yield stress,
$\tau ^e_ {\textit{eff}} = 1.001 Bn$
.
Temporal evolution of
${\textit{Ra}}\,\|{U}\|$
with corresponding colour maps of
${\textit{Ra}} \,U$
at
$(a{,}b{,}c{,}d)$
$Wi = 0.1$
,
$(e{,}f{,}g{,}h)$
$Wi = 0.5$
,
$(i,j,k,l)$
$Wi = 5$
.
$(a)$
Cross markers (
$\times$
) show the viscoplastic solution (Karimfazli et al. Reference Karimfazli, Frigaard and Wachs2015). The time instance of each colour map is marked on
$(a{,}e{,}i)$
. White dashed lines on the colour maps show the yield surface. The
$Wi = 0.1$
case is overlaid for comparison in panels
$(e{,}i)$
, shown by the dashed line. The initial condition is
$T(x,y,t=0) = -x$
.
${\textit{Ra}} = 10^4$
$Pr = 10$
,
$\beta _s = 0.5$
,
${\textit{Bn}} = 0.02$
.
At
$Wi = 0.1$
, the time evolution of the velocity norm is nearly indistinguishable from that of a viscoplastic fluid (the black cross markers represent the velocity norm for a Bingham fluid at the same
${\textit{Bn}}$
,
${\textit{Ra}}$
and
${\textit{Pr}}$
). Similarly, the snapshots of speed closely resembles the predictions of Karimfazli et al. (Reference Karimfazli, Frigaard and Wachs2015) based on the Bingham model. At the onset of flow, the material remains primarily unyielded at the centre. This unyielded region is surrounded by a yielded layer, allowing it to undergo solid-body rotation. Smaller quiescent unyielded zones are also observed at the cavity corners as a result of low shear zones in these regions. As the flow evolves, the shear layer expands while the velocity norm decays to its steady-state value after an initial rapid increase (see the viscoplastic overshoot indicated in figure 5
a).
As
${\textit{Wi}}$
increases, the initial rapid increase in the velocity norm becomes more pronounced and occurs more quickly (see figure 5
e,i). The evolution of
$\|{U}\|$
for the
$Wi = 0.1$
case (serving as an approximation of the viscoplastic case), is shown by the dashed line on figure 5(e,i) for comparison. The snapshots in the second and last columns of figure 5 further confirm that the influence of
${\textit{Wi}}$
on the velocity magnitude, as well as the morphology of the yielded regions, is significantly more pronounced at early times compared with the steady state. Specifically, as
${\textit{Wi}}$
increases, both the size of the unyielded regions and the speed appear to increase. This is in apparent contrast to the mechanistic understanding of yield stress materials based on viscoplastic models. We attribute this to a shear-thinning effect in the Saramito model as
${\textit{Wi}}$
increases (Saramito Reference Saramito2007), consistent with similar trends in velocity magnitude reported by Chauhan et al. (Reference Chauhan, Sahu and Sasmal2021). The interplay between the increased deformability of the unyielded regions at higher
${\textit{Wi}}$
and the modulation of the shear layer due to shear-thinning effects alters the morphology of the central pseudoplug.
To better illustrate the influence of
${\textit{Wi}}$
on the dynamics discussed above, figure 6 presents the time evolution of the velocity norm (solid lines) alongside the normalised area of the yielded material,
$\mathcal{A}_{\!f}$
(dashed lines) for
${\textit{Bn}}=\{0.01, 0.015, 0.02\}$
at different values of
${\textit{Wi}}$
. Here,
$\mathcal{A}_{\!f}$
represents the estimated proportion of yielded material (where
$\tau ^e_ {\textit{eff}} \gt Bn$
). The results indicate that the initial overshoot in the velocity norm becomes more pronounced with increasing
${\textit{Wi}}$
and
${\textit{Bn}}$
.
Evolution of
${\textit{Ra}} \|{U}\|$
(solid line) and
$\mathcal{A}_{\!f}$
(dashed line) as a function of
${\textit{Wi}}$
at
$(a)$
${\textit{Bn}} = 0.01$
,
$(b)$
${\textit{Bn}} = 0.015$
,
$(c)$
${\textit{Bn}} = 0.02$
. The initial condition is
$T(x,y,t=0) = -x$
.
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$\beta _s = 0.5$
.
It also appears in figure 6
$(a)$
, that increasing
${\textit{Wi}}$
from
$Wi=0.1$
to
$Wi=0.5$
slightly slows down yielding. At the highest
${\textit{Wi}}$
considered (
$Wi=5$
), however, yielding appears to be delayed. The apparent delayed onset of yielding is consistently observed across different values of
${\textit{Bn}}$
.
A mechanistic understanding of the observed behaviour can be developed by considering the Saramito model’s mechanical analogy (figure 2). At
$t \to 0^+$
, the system is asymptotically equivalent to a Kelvin-Voigt model comprising a spring with elastic modulus
$\hat {G}$
and a viscous dashpot with viscosity
$\hat {\mu }_s$
. The imposed stress,
$\hat {\sigma }_B\sim \hat {\rho }_0\hat {g} \hat {\alpha }_T\Delta \hat {T} \hat {L}$
, arises from buoyancy and remains steady until motion starts. Yielding begins when the stress carried by the spring exceeds a critical threshold as, initially, the spring and the friction element are in series and carry the same stress.
At
$t \to 0^+$
, the stresses carried by the spring (
$\hat {\sigma }_S$
) and the dashpot (
$\hat {\sigma }_{Ds}$
) are given by
\begin{align} & \hat {\sigma }_S \sim \hat {\sigma }_B \big [ 1- \exp \big ({-}\hat {t}\hat {G}/\hat {\mu }_s\big ) \big ], \ & \hat {\sigma }_{Ds} \sim \hat {\sigma }_B \exp \big ({-}\hat {t}\hat {G}/\hat {\mu }_s\big ). \end{align}
This scaling analysis highlights that initially, the buoyancy stress is carried primarily by the solvent viscosity. The system deviates from the Kelvin-Voigt analogy when the stress in the spring, and consequently in the friction element, exceeds the yield stress, satisfying
\begin{align} \hat {\tau }_y \lesssim \hat {\sigma }_B \big [ 1- \exp \big ({-}\hat {t}\hat {G}/\hat {\mu }_s\big ) \big ]. \end{align}
The dimensionless time scale for yielding onset,
$t_{\textit{yield}}$
, is thus estimated as
\begin{align} t_{\textit{yield}} \sim -\textit{Wi} \ln \left (1-\frac {Bn}{Bn_{\textit{cr}}}\right ). \end{align}
This analysis highlights the interplay between yield stress and elasticity in governing both the deformation prior to yielding and the yielding onset time. A finite yield stress is required for flow onset to be delayed, with the delay time increasing logarithmically with
${\textit{Bn}}$
. Additionally, at finite
${\textit{Bn}}$
, the yielding delay scales linearly with
${\textit{Wi}}$
. Physically, this occurs because a higher
${\textit{Wi}}$
corresponds to a lower elastic modulus, allowing the spring to undergo greater deformation under a given stress. As a result, the dashpot (solvent) deforms more extensively and rapidly, sustaining the buoyancy stress for a longer duration before sufficient stress is imposed on the friction element to trigger yielding. Effectively, the retardation time
$\hat \mu _s/\hat G$
increases, retarding the rate of stress growth in the spring element, thus delaying yielding.
This mechanism also explains the pronounced initial increase in
$\|{U}\|$
. First, at higher
${\textit{Wi}}$
, the spring undergoes more extensive deformation for a given imposed stress. Second, the stress imposed on the spring before yielding increases with
${\textit{Bn}}$
. The combined effect of this pre-yield deformation and the subsequent yielding leads to the accentuated overshoot observed in
$\|{U}\|$
.
Furthermore, as
${\textit{Bn}}$
increases, the evolution of the yielded region, characterised by
$\mathcal{A}_{\!f}$
, exhibits oscillations at the onset of yielding (figure 6
b,c), with the oscillation frequency decreasing as
${\textit{Wi}}$
increases. This behaviour is expected, as a higher
${\textit{Wi}}$
corresponds to a lower elastic modulus, which in turn reduces the natural frequency. Conversely, the oscillation magnitude increases with
${\textit{Bn}}$
, which we attribute to the enhanced damping effect of the friction element, amplifying the elasto-inertial oscillations.
We compare the scaling relation for
$t_{\textit{yield}}$
with our numerical results in figure 7, which shows
$t_{\textit{yield}}$
as a function of (a)
${\textit{Wi}}$
and (b)
${\textit{Bn}}$
. Here,
$t_{\textit{yield}}$
is defined as the time at which
$\mathcal{A}_{\!f}(t = t_{\textit{yield}}) = 5 \times 10^{-3}$
. Bearing in mind the simplifying assumptions of the Saramito model’s one-dimensional mechanical analogue, the scaling analysis is found to show good quantitative agreement with the simulation results. In figure 7(a),
$t_{\textit{yield}}$
exhibits a linear dependence on
${\textit{Wi}}$
, consistent with the scaling prediction. Figure 7(b) shows that the logarithmic dependence on
${\textit{Bn}}$
captures the salient features of the numerical results, namely, a gradual increase in
$t_{\textit{yield}}$
at low
${\textit{Bn}}$
and a sharp rise as
${\textit{Bn}} \to Bn_{\textit{cr}}$
.
Comparison of scaling analysis predictions (
$t_{\textit{yield}} \sim -\textit{Wi} \ln {(1 - Bn/Bn_{\textit{cr}})}$
) with our numerical results for
$(a)$
${\textit{Wi}}$
,
$(b)$
${\textit{Bn}}$
. Here
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$\beta _s = 0.5$
,
$(a)$
${\textit{Bn}} = 0.02$
,
$(b)$
$Wi = 5$
.
Finally, figure 8 shows the time evolution of
$\|{U}\|$
and
$\|{\theta }\|$
(where
$\theta = T + x$
) to steady state for
${\textit{Bn}} = \{0.01, 0.015, 0.02\}$
and
$0.1 \leqslant {\textit{Wi}} \leqslant 5$
. Consistent with the preceding discussion, the influence of elasticity on the steady-state features appears negligible; a notable difference is only observed at the largest
${\textit{Wi}}$
. However, the time evolution of
$\|{U}\|$
shows increased material deformation during the early stages of the system response as pre-yield motion intensifies.
Evolution of
$(a{,}b{,}c)$
${\textit{Ra}} \|{U}\|$
and
$(d{,}e{,}f)$
$\|{\theta }\|$
, as a function of
${\textit{Wi}}$
at
$(a{,}d)$
${\textit{Bn}} = 0.01$
,
$(b{,}e)$
${\textit{Bn}} = 0.015$
,
$(c{,}f)$
${\textit{Bn}} = 0.02$
. The initial condition is
$T(x,y,t=0) = -x$
.
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$\beta _s = 0.5$
.
3.1.2.
$T(t=0) = -x$
and
${\textit{Bn}} \gt Bn_{\textit{cr}}$
When
${\textit{Bn}} \gt Bn_{\textit{cr}}$
, asymptotic analysis based on viscoplastic models predicts that the material remains rigid within the cavity (Karimfazli et al. Reference Karimfazli, Frigaard and Wachs2015; Karimfazli & Frigaard Reference Karimfazli and Frigaard2016). In this section we examine the dynamics in EVP fluids using a representative case where
${\textit{Bn}} = 0.04 \gt Bn_{\textit{cr}} =1/32$
.
As discussed in the previous section, at
$t\rightarrow 0^+$
, the imposed stress is initially supported by elasticity and solvent viscosity (represented by the spring and the dashpot with viscosity
$\hat {\mu }_s$
in the mechanical model shown in figure 2). Consequently, the stress applied to the viscoplastic module (comprising the friction element and the dashpot with viscosity
$\hat {\mu }_p$
) cannot exceed the total imposed buoyancy (see (3.3)). Regardless of the transient dynamics, the friction element therefore remains rigid.
Temporal evolution of
${\textit{Ra}}\,\|{U}\|$
with corresponding colour maps of
${\textit{Ra}} \,U$
at
$(a{,}b{,}c{,}d)$
$Wi = 0.1$
,
$(e{,}f{,}g{,}h)$
$Wi = 0.5$
,
$(i{,}j{,}k{,}l)$
$Wi = 5$
at supercritical
${\textit{Bn}} = 0.04$
. The time instance of each colour map is marked on
$(a{,}e{,}i)$
. The initial condition is
$T(x,y,t=0) = -x$
.
${\textit{Ra}} = 10^4$
$Pr = 10$
,
$\beta _s = 0.5$
.
Figure 9 displays colour maps of speed, overlaid with quiver plots of the velocity vector at
${\textit{Ra}} = 10^4$
,
${\textit{Bn}} = 0.04$
for three values of
$Wi =\{0.1, 0.5, 5\}$
. Each row corresponds to a different
${\textit{Wi}}$
. The first column shows the time evolution of the
$L^2$
norm of velocity, with the red markers indicating the time instances of the snapshots on each row.
All three illustrative cases reveal the immediate onset of systematic oscillatory motion at
$t = 0$
; that is,
$\|{U}\| \gt 0$
. The amplitude of
$\|{U}\|$
oscillations increases with
${\textit{Wi}}$
, while the frequency decreases. The temporal evolution of
$\|{U}\|$
is reminiscent of damped oscillations. The snapshots shown in the last three columns reveal material circulation with a periodically alternating direction. The transient, damped oscillatory motion is captured in supplementary movie 2; additional supplementary videos are provided in Iqbal (Reference Iqbal2025).
We associate the observed oscillation with the elasto-inertial response of the system, damped by the solvent viscosity. In this view, it is expected that the oscillation frequency decreases with decreasing elastic modulus, which corresponds to an increase in
${\textit{Wi}}$
. On the other hand, the oscillation amplitude is expected to increase with
${\textit{Wi}}$
, since the material with a lower elastic modulus deforms more extensively under a given imposed stress (buoyancy here). Similarly, the damping of these oscillations is well explained by dissipation due to solvent viscosity.
To further support this hypothesis, figure 10
$(a{-}c)$
display the time evolution of
$\|{U}\|$
over longer time periods on log-linear axes. The exponential decay of the oscillation amplitudes is clear in all three cases. It can also be seen that the decay time scale is almost identical for the three cases, consistent with the hypothesis that the decay is driven by solvent viscosity, which is constant and identical in the three cases. Figure 10
$(d{-}f)$
confirm that the temperature distribution remains dominantly steady and purely conductive, with small changes due to the pre-yield material deformation.
Evolution of
$(a{,}b{,}c)$
${\textit{Ra}}\|{U}\|$
and
$(d{,}e{,}f)$
$\|{\theta }\|$
at supercritical
${\textit{Bn}} = 0.04$
. The initial condition is
$T(x,y,t=0) = -x$
.
$(a{,}d)$
$Wi = 0.1$
,
$(b{,}e)$
$Wi = 0.5$
,
$(c{,}f)$
$Wi = 5$
. Here,
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$\beta _s = 0.5$
.
Lastly, while the motion patterns discussed so far resemble those observed in the natural convection of viscoplastic fluids in a square cavity (Karimfazli et al. Reference Karimfazli, Frigaard and Wachs2015), alternative modes of pre-yield deformation can be excited, for example, by increasing
${\textit{Ra}}$
when
${\textit{Bn}} \gt Bn_{\textit{cr}}$
. This occurs because the pre-yield deformation is governed by the interplay of solvent viscous stresses
$\hat {\sigma }_{Ds}$
, elastic stresses
$\hat {\sigma }_S$
and buoyancy
$\hat {\sigma }_B$
, whereas whether the material ultimately yields depends on the ratio of buoyancy to yield stress, represented by the Bingham number. As an illustrative case, figure 11 presents the motion development at
${\textit{Ra}}=10^6$
,
${\textit{Bn}}=0.04$
and
$Wi=0.5$
. Figure 11
$(a)$
shows the time evolution of
$\|{U}\|$
, with red circles marking the time instances corresponding to the speed colour maps on the second and third rows.
Although material deformation remains oscillatory, distinct deformation patterns emerge, which differ from unicellular convection. For example, in figure 11
$(c)$
the deformation appears perfectly confined within a circular region concentric with the cavity, where eight smaller, distinct cells can be identified. Similarly, figure 11
$(g)$
illustrates a multi-layered convective structure, with different layers rotating in opposite directions. As the deformation amplitude diminishes due to damping, the primary unicellular pattern progressively dominates (figure 11
h,i), while the other modes decay more rapidly. A movie of this case is provided in supplementary movie 3, highlighting the rich dynamics of the damped, elasto-inertial, multi-modal oscillations.
$(a)$
Temporal evolution of
${\textit{Ra}} \, \|{U}\|$
$(b{-}i)$
colour maps showing the temporal evolution of
${\textit{Ra}} \,U$
at
${\textit{Ra}} = 10^6$
at supercritical
${\textit{Bn}} = 0.04$
. The colour maps are arranged in chronological order from
$(b{-}i)$
, and their time instances are marked on
$(a)$
with circular markers. The initial condition is
$T(x,y,t=0) = -x$
. Here,
$Pr = 10$
,
$Wi = 0.5$
,
$\beta _s = 0.5$
.
3.2. Uniform initial temperature,
$T(t=0) = 0$
In this section we examine the transient dynamics of natural convection, starting from an initially quiescent fluid with a uniform temperature distribution corresponding to a common initial condition in experiments. Similar to the linear initial temperature distribution, the steady state in viscoplastic fluids is convective only if the yield stress is less than a critical value (
${\textit{Bn}} \lt Bn_{\textit{cr}}$
) (Karimfazli & Frigaard Reference Karimfazli and Frigaard2016). Moreover, if the yield stress is below the critical threshold, flow onset in viscoplastic materials is predicted to be delayed by a finite time (Karimfazli & Frigaard Reference Karimfazli and Frigaard2016). Here, we explore how the inclusion of the material elasticity influences the transient dynamics following this initial condition. Our analysis is structured in two parts: first, we investigate the system dynamics at subcritical
${\textit{Bn}}$
, followed by an examination of the supercritical case.
3.2.1.
$T(t=0) = 0$
and
${\textit{Bn}} \lt Bn_{\textit{cr}}$
Temporal evolution of
${\textit{Ra}}\,\|{U}\|$
with corresponding colour maps of
${\textit{Ra}} \,U$
at
$(a{,}b{,}c{,}d)$
$Wi = 0.1$
,
$(e{,}f{,}g{,}h)$
$Wi = 0.5$
,
$(i{,}j{,}k{,}l)$
$Wi = 5$
.
$(a)$
The time instance of each colour map is marked on
$(a{,}e{,}i)$
. The
$Wi = 0.1$
case is overlaid for comparison in panels
$(e{,}i)$
, shown by the dashed line. White dashed lines on the colour maps show the yield surface. The initial condition is
$T(x,y,t=0) = 0$
.
${\textit{Ra}} = 10^4$
$Pr = 10$
,
$\beta _s = 0.5$
,
${\textit{Bn}} = 0.01$
.
Figure 12 presents the time evolution of speed at
${\textit{Ra}}=10^4$
,
${\textit{Bn}}=0.01$
for three values of
$Wi = \{0.1, 0.5, 5\}$
. Each row corresponds to a different
${\textit{Wi}}$
. The first column shows the evolution of
$\|{U}\|$
, with red markers indicating the time instances of the snapshots in each row. The insets in these subfigures display
$\|{U}\|$
over a short period immediately after the temperature difference is applied across the cavity. The
$Wi = 0.1$
case is also overlaid on figure 12
$(e,i)$
for comparison with the approximate viscoplastic behaviour. The last three columns show colour maps of speed with overlaid quiver plots of the velocity vector. The dashed lines show the yield surface, representing isocontours of
$\tau ^e_ {\textit{eff}} =1.001Bn$
.
In all cases, motion begins immediately at
$t\rightarrow 0^+$
. A small, initial local maximum appears in
$\|{U}\|$
, becoming more pronounced as
${\textit{Wi}}$
increases. This is followed by a global maximum in
$\|{U}\|$
(the viscoplastic overshoot), after which the speed gradually decays to its steady-state value. Unlike the first local maximum, the amplitude of the global maximum in
$\|{U}\|$
does not change significantly with
${\textit{Wi}}$
. However, comparing similar time instances, the unyielded regions expand as
${\textit{Wi}}$
increases.
For the initial condition considered here, at
$t=0$
, the strongest buoyancy forces are at the vertical boundaries, where the temperature gradient is largest. In contrast to § 3.1, buoyancy evolves over time, diffusing across the cavity by conduction before yielding onset. This difference explains the distinct early time deformation patterns observed under the two initial conditions. When the initial temperature distribution is linear (see figure 5
b,f,j), buoyancy is uniform, leading to approximately circular motion with little azimuthal speed variation. Conversely, when the initial condition is uniform, buoyancy is highest near the vertical walls, resulting in larger speeds in these regions (figure 12
b,f,j). The latter figures also indicate that at this time, shear stress remains below the yield stress almost everywhere in the domain as the material remains mostly unyielded. The dynamics of the transient natural convection flow under a uniform temperature distribution and subcritical yield stress are illustrated in supplementary movie 4.
The immediate onset of motion contrasts with theoretical and numerical predictions based on viscoplastic models (Karimfazli & Frigaard Reference Karimfazli and Frigaard2016). However, both the immediate onset and the evolution of
$\|{U}\|$
are in qualitative agreement with experimental observations using Carbopol as the working fluid (Jadhav et al. Reference Jadhav, Rossi and Karimfazli2021).
To further illustrate the influence of elasticity on flow dynamics, figure 13 shows the time evolution of
$\|{U}\|$
(solid lines) and
$\mathcal{A}_{\!f}$
(dashed lines) for
${\textit{Bn}}=\{0.01, 0.015, 0.02\}$
at
$Wi=\{0.5,5\}$
. These figures confirm that both the initial deformation and the local maximum in
$\|{U}\|$
occur while the material is primarily unyielded (
$\mathcal{A}_{\!f}\ll 1$
), indicating that viscoelastic solid behaviour dominates at this stage. The higher values of
$\|{U}\|$
and speed at larger
${\textit{Wi}}$
are due to the increased material deformation as the elastic modulus decreases. We note, however, that the magnitude of the local maximum in
$\|{U}\|$
does not change notably with
${\textit{Bn}}$
.
Drawing analogies with the mechanical model shown in figure 2, the friction element remains primarily rigid at early times. During this period, the material behaviour is governed predominantly by the Kelvin-Voigt model. The stresses carried by the spring and the dashpot corresponding to the solvent viscosity may be described by (3.3), noting that here
$\hat {\sigma }_B$
slowly changes with time. This is because the time scale of thermal diffusion (
$\hat {t}_{\textit{TD}} \sim \hat {L}^2/\hat {\kappa }$
) is significantly larger than the retardation time scale (
$\hat {t}_{\mathrm{r}} \sim \hat {\mu }_s/\hat {G}$
); that is,
$\hat {t}_{\textit{TD}}/ \hat {t}_{\mathrm{r}} \sim \textit{Ra}/\textit{Wi} \gg 1$
.
It follows that
\begin{align} -\textit{Wi} \ln \left ( 1- \frac {\textit{Bn}}{{\textit{Bn}}_{\textit{cr}}}\right ) \lt t_{\textit{yield}}\, . \end{align}
As before, the delay in yielding onset increases with
${\textit{Bn}}$
. This is consistent with the
$\mathcal{A}_{\!f}$
trends in figure 13. Overall, pre-yield dynamics are primarily governed by the material’s elastic modulus and solvent viscosity while the yielding time is regulated by the yield stress and elastic modulus. Additionally, the frequency of the oscillations observed in the time evolution of
$\mathcal{A}_{\!f}$
at early times decreases with
${\textit{Wi}}$
while their amplitude increases with
${\textit{Bn}}$
.
Figure 14 shows the time evolution of
$\|{U}\|$
towards the steady state for
${\textit{Bn}} = \{0.01, 0.015, 0.02\}$
and
$0.1 \leqslant {\textit{Wi}} \leqslant 5$
. As before, the influence of elasticity is most pronounced during the transient dynamics, particularly near the time of motion onset. The steady state, however, appears largely unaffected by variations in
${\textit{Wi}}$
, with a notable difference observed only at
$Wi = 5$
, and this difference becomes more pronounced with increasing
${\textit{Bn}}$
.
Evolution of
$\|{U}\|$
, normalised with its time history maximum (solid line) and evolution of normalised yielded volume,
$\mathcal{A}_{\!f}$
(dashed line) as a function of
${\textit{Wi}}$
at
$(a)$
${\textit{Bn}} = 0.01$
,
$(b)$
${\textit{Bn}} = 0.015$
,
$(c)$
${\textit{Bn}} = 0.02$
. The initial condition is
$T(x,y,t=0) = 0$
. Here
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$\beta _s = 0.5$
.
The influence of
${\textit{Wi}}$
on the evolution of
$\|{U}\|$
at
$(a)$
${\textit{Bn}} = 0.01$
,
$(b)$
${\textit{Bn}} = 0.015$
,
${\textit{Bn}} = 0.02$
. The initial condition is
$T(x,y,t=0) = 0$
. Here
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$\beta _s = 0.5$
.
3.2.2.
$T(t=0) = 0$
and
${\textit{Bn}} \gt Bn_{\textit{cr}}$
When
${\textit{Bn}} \gt Bn_{\textit{cr}}$
and the fluid is initially at a uniform temperature equal to the average of the sidewall temperatures, viscoplastic fluids are predicted to remain rigid within the cavity (Karimfazli & Frigaard Reference Karimfazli and Frigaard2016). In this section we examine the dynamics of EVP fluids for a representative case with
${\textit{Bn}} = 0.04 \gt Bn_{\textit{cr}} = 1/32$
.
As discussed in previous sections, at
$t \rightarrow 0^+$
, the imposed stress is initially supported by elasticity and solvent viscosity. Consequently, the stress transmitted to the viscoplastic module cannot exceed the total imposed buoyancy (see (3.3)). Regardless of the transient dynamics, the friction element thus remains rigid.
Figure 15 presents colour maps of speed overlaid with quiver plots of the velocity vector at
${\textit{Ra}} = 10^4$
,
${\textit{Bn}} = 0.04$
for three values of
$Wi = \{0.1, 0.5, 5\}$
. Each row corresponds to a different
${\textit{Wi}}$
. The first column shows the time evolution of
$\|{U}\|$
, with red markers indicating the time instances corresponding to the snapshots in each row.
Temporal evolution of
${\textit{Ra}}\,\|{U}\|$
with corresponding colour maps of
${\textit{Ra}} \,U$
at
$(a{,}b{,}c{,}d)$
$Wi = 0.1$
,
$(e{,}f{,}g{,}h)$
$Wi = 0.5$
,
$(i{,}j{,}k{,}l)$
$Wi = 5$
at supercritical
${\textit{Bn}} = 0.04$
. The time instance of each colour map is marked on
$(a{,}e{,}i)$
. The initial condition is
$T(t=0) = 0$
. Here
${\textit{Ra}} = 10^4$
$Pr = 10$
,
$\beta _s = 0.5$
.
Many features are qualitatively similar to the case where
${\textit{Bn}} \gt Bn_{\textit{cr}}$
and
$T(t = 0) = -x$
(see § 3.1.2). In all three cases, oscillatory motion begins immediately. The amplitude of the oscillations increases with
${\textit{Wi}}$
, while their frequency decreases. The overall behaviour resembles damped elasto-inertial oscillations – the transient dynamics of this can be viewed in supplementary movie 5. The system eventually approaches a quiescent steady state with a purely conductive, linear temperature distribution as the oscillations are damped by solvent viscosity. All evidence supports the hypothesis that the observed motion is pre-yield and driven by viscoelastic deformation. Nonetheless, there are key differences compared with the case with an initial linear temperature distribution, which we highlight below.
First, the snapshots in the second column of figure 15 reveal a multicellular motion pattern, distinct from the unicellular structure observed when
$T(t = 0) = -x$
(see figure 9). We attribute this to the non-uniform initial distribution of the applied body force: buoyancy is largest near the vertical walls and approximately zero at the centre of the domain. As the system approaches the steady state and the buoyancy distribution becomes more uniform, the unicellular motion pattern re-emerges (see the last two columns in figure 15).
Second, figure 15 shows that the direction of circulation remains unchanged as the system evolves toward steady state, unlike the case where
$T(t = 0) = -x$
; that is, once the unicellular motion emerges, it remains clockwise. Comparing the amplitude of oscillations in
$\|{U}\|$
for the two initial conditions (figures 9 and 15), we see that the amplitude is an order of magnitude larger when
$T(t = 0) = -x$
, indicating significantly stronger inertia. Using the mechanical analogy of the Saramito model, we hypothesise that the system is partially over-damped: the EVP fluid is under-damped near the walls where buoyancy is high and over-damped elsewhere. The net result is some oscillation in
$\|{U}\|$
, but not strong enough to reverse the motion direction. The long-time evolution of
$\|{U}\|$
, shown in figure 16, confirms that oscillations decay exponentially, followed by a monotonic exponential decay of
$\|{U}\|$
, indicating dissipation of kinetic energy by solvent viscosity.
Evolution of
$\|{U}\|$
at supercritical
${\textit{Bn}} = 0.04$
for
$(a)$
$Wi = 0.1$
,
$(b)$
$Wi = 0.5$
,
$(c)$
$Wi = 5$
. The initial condition is
$T(t=0) = 0$
. Here
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$\beta _s = 0.5$
.
Oscillation frequency at supercritical
${\textit{Bn}} = 0.04$
as a function of
$(a)$
${\textit{Wi}}$
,
$(b)$
$Wi^{-1}$
for
$T(x,y,t=0)=-x$
(square markers) and
$T(x,y,t=0) = 0$
(circular markers). Here
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$\beta _s = 0.5$
.
Finally, figure 17 displays the dominant frequency of the oscillations’ Fourier spectrum for both initial conditions at
${\textit{Bn}} = 0.04$
(the supercritical regime). The dominant frequency decreases with
${\textit{Wi}}$
, and is lower when
$T(t = 0) = 0$
, as expected in a more strongly damped system. Moreover, the frequency decrease closely follows a power-law fit (
$f\propto Wi^{-1/2}$
), confirming the dominant role of elasticity in the observed oscillatory dynamics. We note this power-law dependence is reminiscent of the frequency scaling with spring stiffness in a damped mass-spring system.
4. Summary and discussion
We have studied the natural convection of EVP fluids in a square cavity with differentially heated side walls using two-dimensional direct numerical simulations based on the Saramito constitutive equation. Our objective was to investigate whether fluid elasticity can help bridge the gap between theoretical/numerical predictions and experimental observations of natural convection in yield stress fluids. The viscoplastic behaviour is recovered asymptotically as the elastic modulus approaches infinity. Alternatively, while the material remains unyielded, the Saramito model behaves as a Kelvin-Voigt solid.
We examined the transient dynamics starting from two initial temperature fields: (i) a linear temperature profile (
$T_{\textit{lin}}$
), which facilitates comparison with theoretical and numerical predictions based on viscoplastic models (Karimfazli et al. Reference Karimfazli, Frigaard and Wachs2015); and (ii) a uniform temperature profile (
$T_{\textit{uni}}$
), representative of experimental configurations (Jadhav et al. Reference Jadhav, Rossi and Karimfazli2021).
Overall, the steady-state behaviour is not significantly affected by elasticity over the range of parameters considered (see table 1). However, elasticity introduces notable differences in the transient dynamics. For both initial conditions, we characterised the transient behaviour in regimes where
${\textit{Bn}} \lt Bn_{\textit{cr}}$
and
${\textit{Bn}} \gt Bn_{\textit{cr}}$
, corresponding respectively to convective and quiescent steady states. The observed phenomena are interpreted using analogies with mechanical models that correspond to the Saramito constitutive equation.
For
$T_{\textit{lin}}$
and
${\textit{Bn}} \lt Bn_{\textit{cr}}$
, the velocity overshoot predicted by viscoplastic models is amplified by the inclusion of elastic effects. The overshoot amplitude increases with increasing
${\textit{Wi}}$
. Analysis of the time evolution of the yielded material area (
$\mathcal{A}_{\!f}$
) reveals that motion begins immediately, but yielding is delayed. This amplified overshoot is thus associated with viscoelastic deformation: increasing
${\textit{Wi}}$
corresponds to a decreasing elastic modulus, allowing greater deformation and, consequently, a higher velocity peak.
The delayed onset of yielding for
$T(t=0) = T_{\textit{lin}}$
and
${\textit{Bn}} \lt Bn_{\textit{cr}}$
contrasts with the behaviour predicted by viscoplastic models, which show immediate yielding (Karimfazli et al. Reference Karimfazli, Frigaard and Wachs2015). As
$t \rightarrow 0^+$
, the material behaves according to the Kelvin–Voigt model: the initially unyielded fluid corresponds to a rigid friction element in the mechanical analogue. At early times, the imposed body force (buoyancy) is balanced by the viscous stresses of the solvent. As deformation increases, more of the buoyancy force is balanced by the elastic component (spring). The spring and friction elements are effectively in series before yielding occurs. Therefore, yielding only begins once the stress transmitted through the spring exceeds the yield stress. A scaling analysis based on this analogy reveals that the yielding delay time scales linearly with
${\textit{Wi}}$
and logarithmically with
${\textit{Bn}}$
.
When
$T(t=0) = T_{\textit{lin}}$
and
${\textit{Bn}} \gt Bn_{\textit{cr}}$
, immediate motion is observed, and the material remains unyielded as it approaches the steady state. The deformation morphology resembles unicellular convection with a periodically alternating direction of circulation. These oscillations decay exponentially, with amplitudes increasing and frequencies decreasing with
${\textit{Wi}}$
. All evidence indicates a pre-yield regime governed by the Kelvin–Voigt model.
The immediate onset of motion is also observed when
$T(t=0) = T_{\textit{uni}}$
and
${\textit{Bn}} \lt Bn_{\textit{cr}}$
. A small local maximum in
$\|{U}\|$
appears prior to the global maximum and the approach to steady state. This immediate onset contrasts with predictions from viscoplastic models, where motion is delayed by a finite time (Karimfazli & Frigaard Reference Karimfazli and Frigaard2016). Nonetheless, the observed onset and evolution of
$\|{U}\|$
are qualitatively consistent with experimental data for Carbopol gels (Jadhav et al. Reference Jadhav, Rossi and Karimfazli2021). By tracking the time evolution of
$\mathcal{A}_{\!f}$
, we distinguish between the onset of motion and the onset of yielding – a distinction that is not possible in viscoplastic models, where the unyielded material behaves as a rigid solid.
We show that the immediate onset of motion corresponds to pre-yield viscoelastic solid deformation, described by the Kelvin–Voigt model. As expected, the amplitude of the initial (pre-yield) local maximum in
$\|{U}\|$
increases with
${\textit{Wi}}$
due to increased deformation at a lower elastic modulus. As with the case of
$T_{\textit{lin}}$
and
${\textit{Bn}} \lt Bn_{\textit{cr}}$
, the onset of yielding is delayed, and the delay increases with both
${\textit{Bn}}$
and
${\textit{Wi}}$
. Pre-yield oscillations may also occur, with amplitudes increasing and frequencies decreasing with
${\textit{Wi}}$
. We note that our predictions for pre-yield motion and yielding delay are based on the Saramito model. This model is well suited to our objective of clarifying the role of elasticity in natural convection of yield stress fluids, owing to its predictive capability for EVP flows (Fraggedakis, Dimakopoulos & Tsamopoulos Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016), robustness and numerical tractability. As a direction for future work, it would be valuable to examine how newer EVP models incorporating mechanisms such as kinematic hardening (Dimitriou & McKinley Reference Dimitriou and McKinley2019; Varchanis et al. Reference Varchanis, Makrigiorgos, Moschopoulos, Dimakopoulos and Tsamopoulos2019), pre-yield plastic deformation (Dimitriou & McKinley Reference Dimitriou and McKinley2019; Varchanis et al. Reference Varchanis, Makrigiorgos, Moschopoulos, Dimakopoulos and Tsamopoulos2019; Kamani, Donley & Rogers Reference Kamani, Donley and Rogers2021) and rate-dependent relaxation time or thixotropy (de Souza Mendes Reference de Souza Mendes2011; Dimitriou & McKinley Reference Dimitriou and McKinley2019; Varchanis et al. Reference Varchanis, Makrigiorgos, Moschopoulos, Dimakopoulos and Tsamopoulos2019; Kamani et al. Reference Kamani, Donley and Rogers2021), influence these pre-yield dynamics.
When
$T(t=0) = T_{\textit{uni}}$
and
${\textit{Bn}} \gt Bn_{\textit{cr}}$
, the response is qualitatively similar to the case with
$T_{\textit{lin}}$
and
${\textit{Bn}} \gt Bn_{\textit{cr}}$
: immediate motion of the unyielded material is observed, with an oscillatory evolution of
$\|{U}\|$
. However, the initial deformation morphology differs from the unicellular pattern. Early time motion exhibits multicellular structures due to the non-uniform buoyancy field. As the temperature field evolves toward steady state, the classical unicellular motion pattern is recovered. Furthermore, the oscillation amplitude is an order of magnitude smaller than in the
$T_{\textit{lin}}$
case, and the circulation direction remains clockwise throughout. Overall, in contrast with
$T_{\textit{lin}}$
, where the system appears under-damped, the dynamics for
$T_{\textit{uni}}$
suggest a partially over-damped system.
Finally, we have shown that over a wide range of
${\textit{Wi}}$
and
${\textit{Bn}}$
, the steady-state solution is not highly sensitive to the inclusion of elasticity. However, elastic effects substantially influence the transient dynamics. Critically, a quiescent Saramito fluid shows immediate motion onset when subject to a change in applied stress, regardless of the magnitude of the yield stress. In our model problem, this corresponds to the introduction of buoyancy at
$t = 0$
.
Unlike viscoplastic models, EVP fluids exhibit a one-way coupling between motion and yielding: deformation can occur prior to yielding. This distinction reconciles the apparent contradiction between predictions by viscoplastic models and experimental observations of natural convection in yield stress fluids in a differentially heated cavity. In EVP fluids, the onset of motion is immediate (as seen in experiments), while the yielding is delayed by a finite time (as predicted by viscoplastic models).
The insights from our study have implications for the modulation of highly elastic yield stress fluids in naturally convecting flows (Abdelgawad et al. Reference Abdelgawad, Haward, Shen and Rosti2024), enhanced understanding of start-up flow of yield stress fluids coupled with heat transfer (Divoux, Barentin & Manneville Reference Divoux, Barentin and Manneville2011; Coussot Reference Coussot2014) and offer insights for transient EVP flows in fast time scale control applications (Karimfazli et al. Reference Karimfazli, Frigaard and Wachs2015). More broadly, our findings contribute to the extensive body of literature on natural convection, highlighting the importance of elasticity in accurate modelling of EVP materials.
Building on our results, several observations emerge as promising directions for future work. For the parameter ranges considered here, the pre-yield decay is typically oscillatory, and the system often behaves as at least partially under-damped. These oscillations are not evident in the available experimental data (Jadhav et al. Reference Jadhav, Rossi and Karimfazli2021), likely for two reasons: first, the small amplitude and rate of deformation may be below current detection thresholds; second, it is known that the Saramito model does not fully replicate the creep response of commonly tested simple yield stress fluids (Fraggedakis et al. Reference Fraggedakis, Dimakopoulos and Tsamopoulos2016), which may be reproduced by more recently proposed EVP models (Dimitriou & McKinley Reference Dimitriou and McKinley2019; Varchanis et al. Reference Varchanis, Makrigiorgos, Moschopoulos, Dimakopoulos and Tsamopoulos2019; Kamani et al. Reference Kamani, Donley and Rogers2021). Encouragingly, while several pieces of the puzzle are coming together, others remain to be discovered and placed.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.10732.
Acknowledgements
The simulations were performed with computation time provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National Infrastructure for Computing (SNIC).
Funding
This work is supported by the funding from the European Research Council (ERC) through the Starting Grant MUCUS (grant no. 2019-StG-852529).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Steady-state dynamics
Here, we analyse the steady-state dynamics of the natural convection. We examine the effect of varying
${\textit{Ra}}$
,
${\textit{Bn}}$
and
${\textit{Wi}}$
on the flow dynamics and heat transfer. The heat transfer at the side walls of the cavity is quantified using the local Nusselt number,
$Nu$
:
\begin{equation} Nu (x=\pm 1/2,y) = -\frac {\partial T}{\partial x} (x=\pm 1/2,y)\, . \end{equation}
We calculate the average Nusselt number,
$\overline {\textit{Nu}}$
, over the side wall:
\begin{equation} \overline {\textit{Nu}} = \int _{-1/2}^{1/2} Nu (x=-1/2,y)\,\mathrm{d}y\, . \end{equation}
This quantifies the magnitude of the average heat transfer over the left side wall, with
$\overline {\textit{Nu}} = 1$
corresponding to pure conduction, and
$\overline {\textit{Nu}} \gt 1$
indicating enhanced heat transfer by convective motion. Furthermore, we define
$\theta = T + x$
to quantify the deviation of the temperature field from that of pure conduction.
Figure 18 shows the variation of
$\overline {\textit{Nu}}$
with
${\textit{Ra}}$
and
${\textit{Bn}}$
at
$Wi=0.5$
. These trends are consistent with previous numerical studies (Turan et al. Reference Turan, Chakraborty and Poole2010; Karimfazli et al. Reference Karimfazli, Frigaard and Wachs2015) on natural convection of viscoplastic fluids, but are briefly summarised here for completeness. Expectedly, heat transfer is enhanced by convection, and therefore,
$\overline {\textit{Nu}}$
increases with
${\textit{Ra}}$
at all Bingham numbers. Conversely, increasing
${\textit{Bn}}$
decreases heat transfer, as unyielded zones form, and the apparent viscosity increases with the yield stress, resulting in less effective convection. At
${\textit{Bn}} = 0.04$
, we find that
$\overline {\textit{Nu}} = 1$
for all
${\textit{Ra}}$
explored here, indicating that heat transfer occurs purely via conduction. We present colour maps of the velocity magnitude (
$U$
) and
$\theta$
fields as a function of
${\textit{Ra}}$
in figure 19, as well as profiles of
$u$
,
$v$
and
$Nu$
along several cross-sections. As convective transport increases with
${\textit{Ra}}$
, the velocity magnitude increases within the convection cell whose shape becomes more skewed, and the thickness of the momentum and thermal boundary layers decreases. Larger convective transport of the temperature can be seen in
$\theta$
(see figure 19
d−f) as temperature is transported by the convection cell from one heated wall to the other, particularly at the lower left and upper right corners, which is also reflected in the
$Nu$
profile in figure 19
$(i)$
. The colour maps are quantitatively similar to figures 15(d,e,f) and 16(d,e,f) in Karimfazli et al. (Reference Karimfazli, Frigaard and Wachs2015), indicating that the steady-state convection of EVP fluids shows good agreement with the viscoplastic solution at this level of elasticity.
We have found
$Wi = 0.5$
a representative value for Carbopol gels in a 15 cm cavity heated with
$\Delta \hat {T} = 30\,\mathrm{K}$
(a similar set-up as in Jadhav et al. Reference Jadhav, Rossi and Karimfazli2021). Using material properties for a 0.09 % w/w Carbopol gel (Lopez et al. Reference Lopez, Naccache, de Souza and Paulo2018):
$\hat {K} = 1.8$
Pa
$\mathrm{s}^n$
(consistency),
$\hat {\tau }_y = 3.0$
Pa,
$n = 0.4$
(power-law index),
$\hat {G}$
= 20.6 Pa (elastic modulus), and thermo-physical properties of water:
$\hat {\rho }_0 = 997$
kg m
$^{-3}$
,
$\hat {\kappa } = 1.4\times 10^{-7}$
$\mathrm{m}^2\,\mathrm{s}^{-1}$
,
$\hat {\alpha }_T = 2.37\times 10^{-4}$
$\mathrm{K}^{-1}$
we can estimate
${\textit{Wi}}$
in a
$\hat {L} = 15$
cm cavity heated with
$\Delta \hat {T} = 30\,\mathrm{K}$
. Using the Herschel-Bulkley constitutive equation,
$\hat \tau = \hat \tau _y + \hat K\dot {\hat {\gamma }}^n \implies \hat \mu _p = \hat K\dot {\hat {\gamma }}^{n-1}$
with
$\dot {\hat {\gamma }} = {\hat u}_{\textit{ref}}/{\hat l}_{\textit{ref}} = \hat \rho _0\hat g\hat {\alpha }_T\Delta \hat T\hat L/\hat \mu _p$
gives
$\hat \mu _p = [\hat K (\hat \rho _0\hat g\hat {\alpha }_T\Delta \hat T\hat L)^{n-1} ]^{1/n} \approx 0.13$
Pa s and
$\dot {\hat {\gamma }}\approx 80.8\, \mathrm{s}^{-1}$
. Therefore,
$Wi = \hat \mu _p \dot {\hat {\gamma }}/\hat G\approx 0.5$
. This indicates that, at steady state, the flow dynamics and heat transfer of Carbopol gels are not significantly modulated by this level of fluid elasticity. Note that this is in qualitative agreement with the PIV flow field measurements of Jadhav et al. (Reference Jadhav, Rossi and Karimfazli2021), which showed classical viscoplastic flow patterns at steady state.
Plot of
$\overline {\textit{Nu}}$
as a function of
${\textit{Ra}}$
and
${\textit{Bn}}$
. The circular markers correspond to
${\textit{Ra}} = 10^2$
, the triangular markers to
${\textit{Ra}} = 10^4$
and the square markers to
${\textit{Ra}} = 10^6$
. The remaining dimensionless numbers are held constant at
$Pr = 10$
,
$Wi = 0.5$
,
$\beta _s = 0.5$
.
$(a{,}b{,}c)$
Colour maps of
${\textit{Ra}}\, U$
at
$(a)$
${\textit{Ra}} = 10^2$
,
$(b)$
${\textit{Ra}} = 10^4$
,
$(c)$
${\textit{Ra}} = 10^6$
;
$(d{,}e{,}f)$
colour maps of
$\theta$
at
$(d)$
${\textit{Ra}} = 10^2$
,
$(e)$
${\textit{Ra}} = 10^4$
,
$(f)$
${\textit{Ra}} = 10^6$
; profiles of
$(g)$
${\textit{Ra}}\,\,u$
at
$x=0$
,
$(h)$
${\textit{Ra}}\,\,v$
at
$y = 0$
,
$(i)$
$Nu$
at
$x=-0.5$
as a function of
${\textit{Ra}}$
. The asterisk markers (
$\ast$
) correspond to the viscoplastic solution. Here
$Pr = 10$
,
$Wi = 0.5$
,
$\beta _s = 0.5$
,
${\textit{Bn}} = 0.01$
.
To study the effect of more elastic EVP fluids on the natural convection dynamics, such as the fluids newly created by Abdelgawad et al. (Reference Abdelgawad, Haward, Shen and Rosti2024), we investigate the effect of
${\textit{Wi}}$
. Figure 20
$(a,b,c)$
shows the variation of the velocity magnitude in the range
$0.1\leqslant Wi\leqslant 5$
. The most significant effect is on the morphology of the unyielded regions. As in the viscoplastic case (Karimfazli et al. Reference Karimfazli, Frigaard and Wachs2015), we find static unyielded zones in the corners and mobile unyielded zones (pseudoplugs) in the convection cell. With increasing
${\textit{Wi}}$
, we find that the six pseudoplugs merge into two elongated unyielded regions. This is owing to two effects. First, increasing
${\textit{Wi}}$
increases the thickness of the momentum and thermal boundary layers (Chauhan et al. Reference Chauhan, Sahu and Sasmal2021), which leads to low shear regions near the enclosure centre in which sufficient stresses do not develop to yield the material. Second, the unyielded regions become less stiff as
${\textit{Wi}}$
increases (decreasing elastic modulus), and hence, deform more under shear without yielding, as previously observed in porous media (Chaparian et al. Reference Chaparian, Izbassarov, Vita, Brandt and Tammisola2020b
; De Vita et al. Reference De Vita, Rosti, Izbassarov, Duffo, Tammisola, Hormozi and Brandt2018) and multiphase flows (Izbassarov & Tammisola Reference Izbassarov and Tammisola2020). In figure 20
$(g,h,i)$
we show the cross-sectional profiles of the velocities and Nusselt number. The profiles at
$Wi=0.1$
and
$Wi=0.5$
are near identical, further confirming that these cases are close to the viscoplastic limit, while the profiles at
$Wi = 5$
show small but discernible changes.
In figure 21 (
${\textit{Bn}} = 0.02$
) we find similar trends as in figure 20 (
${\textit{Bn}} = 0.01$
): there is negligible difference between
$Wi = 0.1$
and
$Wi = 0.5$
in the colour maps and profiles, indicating that we are close to the viscoplastic limit. The unyielded zones have of course increased in size at
${\textit{Bn}}=0.02$
, as seen in figure 21
$(a,b,c)$
, as fewer regions experience sufficient stress to overcome the larger yield stress. As
${\textit{Wi}}$
increases, the unyielded zones in the convection cell expand and merge into one large, mobile plug. We note that
$Wi = 5$
shows a larger difference compared with low
${\textit{Wi}}$
in the velocity and
$Nu$
profiles for
${\textit{Bn}} = 0.02$
(figure 21
g,h,i) compared with
${\textit{Bn}} = 0.01$
(figure 20
g,h,i). The colour maps of
$U$
and
$\theta$
also reflect this: in figures 20 and 21 we see larger speeds in
$U$
, and larger
$\theta$
at
$Wi = 5$
, indicating increased convective transport with
${\textit{Wi}}$
. This is more apparent at larger
${\textit{Bn}}$
, indicating the effect of
${\textit{Wi}}$
becomes stronger with
${\textit{Bn}}$
.
$(a{,}b{,}c)$
Colour maps of
${\textit{Ra}} \|{\boldsymbol{u}}\|$
at
$(a)$
$Wi = 0.1$
,
$(b)$
$Wi = 0.5$
,
$(c)$
$Wi = 5$
;
$(d{,}e{,}f)$
colour maps of
$\theta$
at
$(d)$
$Wi = 0.1$
,
$(e)$
$Wi = 0.5$
,
$(f)$
$Wi = 5$
; profiles of
$(g)$
${\textit{Ra}}\,\,u$
at
$x=0$
,
$(h)$
${\textit{Ra}}\,\,v$
at
$y = 0$
,
$(i)$
$Nu$
at
$x=-0.5$
as a function of
${\textit{Wi}}$
. The remaining dimensionless numbers are held constant at
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$\beta _s = 0.5$
,
${\textit{Bn}} = 0.01$
. The white isocontours represent the yield surface.
$(a{,}b{,}c)$
Colour maps of
${\textit{Ra}}\, U$
at
$(a)$
$Wi = 0.1$
,
$(b)$
$Wi = 0.5$
,
$(c)$
$Wi = 5$
;
$(d{,}e{,}f)$
colour maps of
$\theta$
at
$(d)$
$Wi = 0.1$
,
$(e)$
$Wi = 0.5$
,
$(f)$
$Wi = 5$
; profiles of
$(g)$
${\textit{Ra}}\,\,u$
at
$x=0$
,
$(h)$
${\textit{Ra}}\,\,v$
at
$y = 0$
,
$(i)$
$Nu$
at
$x=-0.5$
as a function of
${\textit{Wi}}$
. Here
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$\beta _s = 0.5$
,
${\textit{Bn}} = 0.02$
. The white lines in (a--c) show the yield surface.
The overall effect of
${\textit{Wi}}$
on the flow and heat transfer at steady state is illustrated in figure 22, showing maximum values
$u$
and
$v$
in the cavity and
$\overline {\textit{Nu}}$
over the side walls, as functions of
${\textit{Wi}}$
for
$(a,b,c)$
${\textit{Bn}} = 0.01$
and
$(d,e,f)$
${\textit{Bn}} = 0.02$
. We observe that
$u_{\textit{max}}$
,
$v_{\textit{max}}$
and
$\overline {\textit{Nu}}$
all increase with
${\textit{Wi}}$
(except for
$v_{\textit{max}}$
, which shows a non-monotonic variation). Similar observations were also found in the transient analysis (see, e.g. figures 5 and 8), where higher velocity magnitudes and system kinetic energy were seen at the largest
${\textit{Wi}}$
. The results, overall, indicate increased convection with
${\textit{Wi}}$
and are largely consistent with the findings by Chauhan et al. (Reference Chauhan, Sahu and Sasmal2021) on natural convection in viscoelastic, FENE-P fluids. Chauhan et al. (Reference Chauhan, Sahu and Sasmal2021) attributed the increased convection at low to moderate
${\textit{Wi}}$
to the associated shear-thinning effect in the FENE-P constitutive equation. The Saramito model shows a similar coupling between
${\textit{Wi}}$
and shear-thinning effects (Saramito Reference Saramito2007). To isolate the effect of elasticity on natural convection, we have performed additional simulations for natural convection in Oldroyd-B fluids, reported in Appendix B. Therein, we find a non-monotonic relationship for
$u_{\textit{max}}$
and
$\overline {\textit{Nu}}$
(an initial increase, followed by a decrease at large
${\textit{Wi}}$
), whereas
$v_{\textit{max}}$
shows monotonic increase. Moreover, figure 24 shows that one must go to significantly large
${\textit{Wi}}$
(
$Wi \gt 100$
) to observe a decrease in
$u_{\textit{max}}$
and
$\overline {\textit{Nu}}$
. Therefore, in the range of
${\textit{Wi}}$
we have investigated (
$Wi\leqslant 5$
) we conclude that elasticity and shear-thinning act in the same direction and work to increase
$u_{\textit{max}}$
,
$v_{\textit{max}}$
and
$\overline {\textit{Nu}}$
.
Steady state
$(a{,}d)$
${\textit{Ra}}\,u_{\textit{max}}$
,
$(b{,}e)$
${\textit{Ra}}\, v_{\textit{max}}$
,
$(c{,}f)$
$\overline {\textit{Nu}}$
as a function of
${\textit{Wi}}$
at
$(a{,}b{,}c)$
${\textit{Bn}} = 0.01$
and
$(d{,}e{,}f)$
${\textit{Bn}} = 0.02$
for natural convection of an EVP fluid. Here
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$\beta _s = 0.5$
.
We observe some differences in the trends with
${\textit{Wi}}$
for EVP fluids that are not seen for viscoelastic Oldroyd-B/FENE-P fluids. Proportionally larger changes in
$u_{\textit{max}}$
and
$\overline {\textit{Nu}}$
with
${\textit{Wi}}$
are observed for the Saramito fluid, as well as a non-monotonic variation of
$v_{\textit{max}}$
. This is due to several reasons. First, the effect of shear thinning is present in the Saramito constitutive equation (Saramito Reference Saramito2007), which accounts for the greater increase in
$u_{\textit{max}}$
and
$\overline {\textit{Nu}}$
compared with Oldroyd-B fluids (compare figures 22 and 24). Second, the Saramito constitutive equation can show enhanced elastic effects at finite
${\textit{Bn}}$
compared with viscoelastic fluids (Chaparian et al. Reference Chaparian, Ardekani, Brandt and Tammisola2020a
; Izbassarov et al. Reference Izbassarov, Rosti, Brandt and Tammisola2021; Habibi et al. Reference Habibi, Iqbal, Ardekani, Chaparian, Brandt and Tammisola2025). This is clearly seen in figures 20, 21 and 22, in which proportionally larger variations in
$u$
,
$v$
and
$Nu$
with
${\textit{Wi}}$
are observed at
${\textit{Bn}} = 0.02$
compared with
${\textit{Bn}} = 0.01$
. Finally, EVP fluids show the formation of yielded and unyielded regions: in particular, in figures 20
$(a,b,c)$
and 21
$(a,b,c)$
we see that the yield surface is close to the edge of the boundary layers. The interactions with the yield surface, we hypothesise, is the reason for the non-monotonic variation of
$v_{\textit{max}}$
in this range of
${\textit{Wi}}$
, which is not seen for Oldroyd-B fluids (figure 24) or FENE-P fluids (Chauhan et al. Reference Chauhan, Sahu and Sasmal2021).
To conclude our analysis, we assess the effect of the different initial temperature distributions (
$T(t=0) = -x$
and
$T(t=0) = 0$
) used in our study on the steady-state flow and heat transfer to identify possible hysteresis effects. We compare the mid-sectional velocity profiles and the Nusselt number profile on the left wall in figure 23 for the linear temperature distribution (solid line) and the uniform distribution (asterisk markers). Furthermore, we compare the velocity norm (
$\|{U}\|$
), the average Nusselt number on the left wall (
$\overline {\textit{Nu}}$
) and the normalised yielded area fraction (
$\mathcal{A}_{\!f}$
) at steady state in table 3. It can be seen that the velocity and
$Nu$
profiles are apparently indistinguishable, and
$\|{U}\|$
,
$\overline {\textit{Nu}}$
and
$\mathcal{A}_{\!f}$
at steady state agree exactly to two decimal places. This indicates that the initial temperature distribution has no significant effect on the steady-state flow and heat transfer and shows negligible hysteresis effects.
Comparison of steady state
$\|{U}\|$
,
$\overline {\textit{Nu}}$
and
$\mathcal{A}_{\!f}$
between
$T(t=0) = -x$
and
$T(t=0) = 0$
initial temperature distributions. Here
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$Wi = 0.5$
,
$\beta _s = 0.5$
,
${\textit{Bn}} = 0.01$
.
Comparison of profiles of
$(a)$
${\textit{Ra}}\, u$
at
$x=0$
,
$(b)$
${\textit{Ra}}\, v$
at
$y=0$
,
$(c)$
$Nu$
at
$x = -1/2$
for
$T(t=0) = -x$
(solid line) and
$T(t=0) = 0$
(asterisk markers
$\ast$
). Here
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$Wi = 0.5$
,
$\beta _s = 0.5$
.
Plots of
$(a)$
${\textit{Ra}}\,u_{\textit{max}}$
,
$(b)$
${\textit{Ra}}\,v_{\textit{max}}$
,
$(c)$
$\overline {\textit{Nu}}$
as a function of
${\textit{Wi}}$
for natural convection of an Oldroyd-B fluid. Here
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$\beta _s = 0.5$
.
Appendix B. The effect of
${\textit{Wi}}$
on natural convection of an Oldroyd-B fluid
Previous numerical studies have studied natural convection of viscoelastic, FENE-P fluids (Chauhan et al. Reference Chauhan, Sahu and Sasmal2021). These authors performed parametric studies of
${\textit{Wi}}$
to investigate the effect of elasticity on steady-state natural convection. However, as the authors discussed, this does not isolate the effect of elasticity due to the coupling of shear thinning with
${\textit{Wi}}$
in the FENE-P constitutive equation. The Saramito constitutive equation used in the present study experiences a similar coupling between shear thinning and
${\textit{Wi}}$
(Saramito Reference Saramito2007). To this end, we perform parametric studies of
${\textit{Wi}}$
for natural convection of an Oldroyd-B fluid to isolate the effect of fluid elasticity on the flow and heat transfer.
We perform the simulations at
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$\beta _s = 0.5$
. In figure 24 we plot the maximum horizontal velocity in the domain,
$u_{\textit{max}}$
, the maximum vertical velocity,
$v_{\textit{max}}$
, and
$\overline {\textit{Nu}}$
over the left side wall as a function of
${\textit{Wi}}$
. We find non-monotonic trends in
$u_{\textit{max}}$
and
$\overline {\textit{Nu}}$
: they increase for
$Wi \leqslant 100$
and decrease thereafter. Here
$v_{\textit{max}}$
shows monotonic increase. The non-monotonic relationship of
$u_{\textit{max}}$
and
$\overline {\textit{Nu}}$
with
${\textit{Wi}}$
corroborates the non-monotonic behaviour in
$\boldsymbol{u}$
and plateau in
$Nu$
observed by Chauhan et al. (Reference Chauhan, Sahu and Sasmal2021) at high
${\textit{Wi}}$
. However, for the range of
${\textit{Wi}}$
explored in the present study (
$Wi \leqslant 5$
), we conclude that
${\textit{Wi}}$
increases all three flow variables and acts in the same direction as shear thinning.
Appendix C. Convergence study and numerical tests
C.1. Convergence study
To verify the convergence of our spatial and temporal discretisation, we have performed convergence studies comparing meshes M1 =
$N_x\times N_y = 256\times 256$
and M3 =
$512\times 512$
. The mesh
$\mathrm{M}2 = 384\times 384$
was used in our study. We also compare our results with a smaller time step corresponding to a Courant--Friedrichs--Lewy (CFL) number of
$0.5$
, while a time-step corresponding to
$\mathrm{CFL}=0.9$
was used in the present work.
In figure 25 we plot the profiles of
$u$
and
$v$
at the vertical and horizontal mid-sections, respectively, and
$Nu$
at the left wall for meshes M2 (solid black line) and M3 (asterisk markers). We find excellent agreement between two meshes, showing that our spatial discretisation accurately resolves the flow variables for a representative case in our problem set-up.
In § 3.1.2, oscillations in the temporal evolution of
$\|{U}\|$
and
$\|{\theta }\|$
were observed at supercritical
${\textit{Bn}}$
i.e. when convection is suppressed by the yield stress. To verify that these oscillations are physical, we perform convergence studies for our spatial and temporal discretisation. In figure 26 we plot a close-up of the
$\|{U}\|$
and
$\|{\theta }\|$
evolution at
${\textit{Ra}} = 10^4$
and
${\textit{Ra}} = 10^6$
in figures 26
$(a,b)$
and 26
$(c,d)$
, respectively. We plot the base mesh, M2, which was used in the present study for reference, as well as the coarser mesh M1, the finer mesh M3 and a case using a lower time-step corresponding to
$\mathrm{CFL}=0.5$
. We find excellent spatial and temporal convergence for the solution of
$\|{U}\|$
across the different meshes and time steps, and only a small deviation in
$\|{\theta }\|$
between meshes M2 and M3.
Based on our convergence study, we find mesh
$\text{M2} =384\times 384$
a suitable mesh for our study.
Grid convergence for profiles of
$(a)$
${\textit{Ra}}\,\,u$
at
$x = 0$
,
$(b)$
${\textit{Ra}}\,\, v$
at
$y=0$
,
$(c)$
$Nu$
at
$x=-0.5$
. The dimensionless numbers are
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$Wi = 0.5$
,
$\beta _s = 0.5$
,
${\textit{Bn}} = 0.01$
. The solid line corresponds to the mesh
$\mathrm{M2} = 384\times 384$
used in the present study and the asterisk
$(\ast )$
markers correspond to mesh
$\mathrm{M3} = 512\times 512$
.
Evolution of
$(a{,}b{,}c)$
$\|{U}\|$
and
$(d{,}e{,}f)$
$\|{\theta }\|$
with time at
$(a{,}b)$
${\textit{Ra}} = 10^4$
,
$(c{,}d)$
${\textit{Ra}} = 10^6$
at three meshes M1, M2, M3 with different spatial resolution at CFL = 0.9, and mesh M2 at CFL = 0.5. Here
$Pr = 10$
,
$Wi = 0.5$
,
$\beta _s = 0.5$
,
${\textit{Bn}} = 0.04$
.
C.2. Effect of temporal resolution
To verify that our temporal resolution is sufficient to capture the fast dynamics of the start-up flow, we performed additional simulations for representative cases using (i) a half-smaller time step,
$\Delta t$
, and (ii) a halved time step combined with a twice larger sampling rate. The results, presented in figure 27, show that the higher temporally resolved simulations are visually indistinguishable from those in the present study. This confirms that the temporal resolution used in our simulations is adequate to resolve the fast transient dynamics of the start-up flow.
Effect of temporal resolution on the evolution of
${\textit{Ra}}\, \|{U}\|$
using a half-smaller time step,
$\Delta t$
, (square markers
$\square$
) and a half-smaller
$\Delta t$
and twice larger sampling rate (circular markers
$\circ$
) compared with the present work (solid line):
$(a)$
$T(t=0) = -x$
,
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$Wi = 5$
,
$\beta _s = 0.5$
,
${\textit{Bn}} = 0.02$
;
$(b)$
$T(t=0) = 0$
,
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$Wi = 5$
,
$\beta _s = 0.5$
,
${\textit{Bn}} = 0.01$
.
C.3. Effect of the extra stress boundary condition
In the present study, homogeneous Neumann boundary conditions were imposed for the extra stress tensor on all walls (i.e.
$\partial \boldsymbol{\tau }^e/\partial n = 0$
). To assess the influence of this, we performed additional simulations for selected representative cases using a linear extrapolation boundary condition for
$\boldsymbol{\tau }^e$
. The results are shown in figure 28, where the homogeneous Neumann condition is labelled ‘N’ and the linear extrapolation condition ‘LE’. The two cases appear visually indistinguishable across both initial temperature distributions and for different Weissenberg numbers. This confirms that the choice of boundary condition for
$\boldsymbol{\tau }^e$
does not introduce numerical artefacts or influence the simulation results.
Effect of the wall boundary condition imposed on the extra stress tensor for
$(a)$
$T(t=0) = -x$
,
$(b)$
$T(t=0) = 0$
. Here ‘N’ denotes the homogeneous Neumann condition used in the present work, shown by the cross markers (
$\times$
), and ‘LE’ denotes the linear extrapolation boundary condition, shown by the solid lines. Here
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$\beta _s = 0.5$
,
$(a)$
${\textit{Bn}} = 0.02$
,
$(b)$
${\textit{Bn}} = 0.01$
.
C.4. Effect of varying
$\beta _s$
Here we examine the effect of different
$\beta _s$
on the transient and steady-state flow. For prototypical yield stress fluids, the plastic viscosity,
$\hat \mu _p$
, is typically much larger than the solvent viscosity,
$\hat \mu _s$
, i.e.
$\beta _s \ll 1$
. We have used a value of
$\beta _s = 0.5$
that is larger than that for typical yield stress fluids to ensure numerical stability (as
$\beta _s\to 0$
, the diffusion term from solvent viscosity in (2.7) disappears and requires particular numerical treatment). Hence, this exploration serves as a sensitivity analysis on the effect of smaller
$\beta _s$
on our results.
Figure 29 shows the temporal evolution of
$\|{U}\|$
as a function of different
$\beta _s$
at
$Wi = 0.5$
,
${\textit{Bn}} = 0.01$
. As
$\beta _s$
decreases, quantitative differences emerge, namely an attenuated initial overshoot that has a lower, earlier peak and followed by a faster decay. This is due to the decreased retardation time,
$\hat \mu _s/\hat G$
, with decreasing
$\beta _s$
: lesser retardation of elastic stresses leads to earlier yielding onset, hence, the lower peak and earlier relaxation in the initial overshoot. Conversely, larger
$\beta _s$
corresponds to larger retardation time that retards elastic stress growth and delays yielding. Nevertheless, as figure 29 demonstrates, these differences are quantitative, not qualitative; immediate motion onset due to pre-yield elastic deformation is observed in all cases, underscoring the role of elasticity in driving the initial transient motion.
Effect of
$\beta _s$
on the transient evolution of
${\textit{Ra}}\, \|{U}\|$
for
$(a)$
$T(t=0) = -x$
,
$(b)$
$T(t=0) = 0$
. The remaining dimensionless numbers are held constant at
${\textit{Ra}} = 10^4$
,
$Pr = 10$
,
$Wi = 0.5$
,
${\textit{Bn}} = 0.01$
.
Open access
