1 Introduction
The interface of two approaching fluids in a porous medium is unstable if a low viscosity fluid displaces a large viscosity one. These instabilities usually result in the development of finger-shaped patterns and are known as viscous fingering (VF) or Saffman–Taylor instabilities (Saffman & Taylor Reference Saffman and Taylor1958). This phenomenon is favourable when, for example, strong mixing of the fluids is desired in low Reynolds number flows, e.g. micro-mixers (Jha, Cueto-Felgueroso & Juanes Reference Jha, Cueto-Felgueroso and Juanes2011), while it is unfavourable when the aim is to sweep the highest amount of the displaced fluid, e.g. oil industry (Doorwar & Mohanty Reference Doorwar and Mohanty2011), chromatographic columns (Rousseaux, Martin & De Wit Reference Rousseaux, Martin and De Wit2011) and contaminant transport (De Wit, Bertho & Martin Reference De Wit, Bertho and Martin2005).
There are two classes of VF: miscible and immiscible. In miscible VF, changes in the properties which lead to the growth of instabilities result from the injection of a solvent fully miscible with the displaced fluid. Consideration of molecular diffusion and mechanical dispersion in miscible displacements distinguishes the system from immiscible displacements where surface tension must be considered instead. The phenomenon has been widely studied and the effects of different parameters were analysed. For instance, in the case of miscible displacements which are the focus of the present study, the effects of diffusivity, heat transfer, chemical reaction, melting, heterogeneity, inlet flow conditions, flow configuration, etc., have been thoroughly studied (Tan & Homsy Reference Tan and Homsy1992; De Wit & Homsy Reference De Wit and Homsy1999; Chen, Chen & Miranda Reference Chen, Chen and Miranda2005; Pritchard Reference Pritchard2009; Mishra et al. Reference Mishra, Trevelyan, Almarcha and De Wit2010; Yuan & Azaiez Reference Yuan and Azaiez2014; Sajjadi & Azaiez Reference Sajjadi and Azaiez2016). The wide usage of nanofluids in different applications, including porous-medium related applications, raises the question of the potential effects of the presence of nanoparticles (NPs) on the instability. Given the fact that the presence of NPs alters the properties of the base fluid, its effect on VF instabilities should not be ignored. However, despite the importance of this question, there are only very few studies on the topic, dealing mainly with miscible displacements. In particular, it has been reported that the addition of NPs to the displacing fluid attenuates the VF instability of initially unstable binary isothermal systems (Ghesmat et al. Reference Ghesmat, Hassanzadeh, Abedi and Chen2011; Dastvareh & Azaiez Reference Dastvareh and Azaiez2017). NP deposition was however reported to enhance the instability. Regarding the effect of NP diffusivity, a linear stability analysis (LSA) conducted by Ghesmat et al. (Reference Ghesmat, Hassanzadeh, Abedi and Chen2011) shows a turning point in the instability of NP-laden unstable systems, where Brownian diffusivity attenuates the instability if it is smaller than the diffusion of the fluids, and makes the system more unstable otherwise. A subsequent study based on both LSA and nonlinear simulations (NLS) revealed that Brownian diffusivity only enhances the instability of initially unstable NP-laden systems and it cannot destabilize initially stable ones (Dastvareh & Azaiez Reference Dastvareh and Azaiez2017). In a more recent study on the effects of NPs at the mesoscopic scale for non-isothermal porous media consisting of regularly arranged impermeable disks, it was found that both the NP size and surface potential contribute to increasing the instability (Zargartalebi & Azaiez Reference Zargartalebi and Azaiez2018). It was also reported that for any given surface potential, NPs are no longer able to affect the instability when their size exceeds a critical diameter.
The objective of this study is to investigate the role of parameters affecting the development and growth of VF instabilities of non-isothermal nanoflows in homogeneous porous media. More specifically, the study will analyse the combined effects of Brownian diffusion and thermophoresis, the average motion of NPs resulting from the temperature gradient, on the flow dynamics, particle distribution and the instability. The use of thermophoretic effects and the potential for the transport of particles by temperature gradients has been applied in a number of processes that include thermal field-flow fractionation, where micro and nano-scale particles can be separated from the solvent by imposing temperature gradient (Giddings Reference Giddings1993), ultrafine particle collection (Wen & Wexler Reference Wen and Wexler2007) and the design of thermophoretic swimmers, where particles are self-propelled as a result of anisotropic heating of the surrounding fluid (Golestanian, Liverpool & Ajdari Reference Golestanian, Liverpool and Ajdari2007). Accordingly, thermophoresis as an additional transport mechanism of NPs may affect the NP-laden VF instability observed in applications such as nano-drug delivery (Ensign, Cone & Hanes Reference Ensign, Cone and Hanes2012), particle mixing in microfluidic devices (Stone, Stroock & Ajdari Reference Stone, Stroock and Ajdari2004) and heavy oil upgrading and recovery in the presence of nano-catalysts (Hashemi, Nassar & Pereira Almao Reference Hashemi, Nassar and Pereira Almao2014). However, this aspect of the flow has not been examined in previous studies dealing with non-isothermal nanoflow displacements in porous media.
In this study, first the effects of Brownian diffusion and thermophoresis are analysed in two representative unstable systems: HDC (hot fluid displaces cold fluid) and CDH (cold fluid displaces hot fluid). The analysis is conducted using the concentration distribution of the displaced fluid, resulting from direct nonlinear simulations of the flow. The observed phenomena are explained in terms of the viscosity and NP concentration distributions and through quantitative analyses. Physical interpretations of the observed trends are presented. Then the study is extended to examine the effects of thermophoresis in the case of initially stable systems.
2 Physical problem
A two-dimensional homogeneous porous medium with constant intrinsic permeability
$K$
or equivalently a horizontal Hele-Shaw cell is considered, where fluid B with initial concentration
$C_{b0}$
, viscosity
$\unicode[STIX]{x1D707}_{b0}=\unicode[STIX]{x1D707}(0,C_{b0},0)$
and temperature
$T_{b}$
is at rest. A miscible fluid A with initial concentration
$C_{a0}$
, viscosity
$\unicode[STIX]{x1D707}_{a0}=\unicode[STIX]{x1D707}(C_{a0},0,0)$
and temperature
$T_{a}$
containing NPs with initial concentration
$C_{n0}$
, is injected at a constant velocity
$U$
to displace fluid B. The viscosity of the resulting initial nanofluid (fluid A containing NPs) is
$\unicode[STIX]{x1D707}_{n0}=\unicode[STIX]{x1D707}(C_{a0},0,C_{n0})$
and the NPs are assumed to be in thermal equilibrium with the host fluid. Furthermore, the fluids are assumed to be incompressible and the nanofluid is dilute. The length and width of the medium are
$L$
and
$W$
as shown in figure 1.
Schematic view of the flow geometry.
3 Problem formulation
The flow is modelled with the following equations:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{V}_{D}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}P=-\frac{\unicode[STIX]{x1D707}}{K}\boldsymbol{V}_{D}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}C_{b}}{\unicode[STIX]{x2202}t}+\frac{1}{\unicode[STIX]{x1D719}}(\boldsymbol{V}_{D}\boldsymbol{\cdot }\unicode[STIX]{x1D735}C_{b})=D_{BA}\unicode[STIX]{x1D6FB}^{2}C_{b}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}C_{n}}{\unicode[STIX]{x2202}t}+\frac{1}{\unicode[STIX]{x1D719}}(\boldsymbol{V}_{D}\boldsymbol{\cdot }\unicode[STIX]{x1D735}C_{n})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(D_{Br}\unicode[STIX]{x1D735}C_{n}+D_{T}\frac{\unicode[STIX]{x1D735}T}{T}\right)-k_{dep}C_{n}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D719}}(\boldsymbol{V}_{D}\boldsymbol{\cdot }\unicode[STIX]{x1D735}T)=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D735}T)+\unicode[STIX]{x1D6FD}\left(D_{Br}\unicode[STIX]{x1D735}C_{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}T+D_{T}\frac{\unicode[STIX]{x1D735}T\boldsymbol{\cdot }\unicode[STIX]{x1D735}T}{T}\right). & \displaystyle\end{eqnarray}$$
The continuity equation and Darcy’s law (3.1), (3.2) are used for the conservation of mass and momentum, where
$\boldsymbol{V}_{D}(u,v)$
is Darcy’s velocity,
$P$
the local pressure,
$\unicode[STIX]{x1D707}$
the viscosity and
$K$
the intrinsic permeability. The flow is also governed by the advection diffusion equations for the transport of fluids and NPs. Since the total mass density of the system is conserved, the transport equations of only two components, e.g. B and the NPs, need to be considered ((3.3) and (3.4)). One expects that in multi-component systems like the one considered here, the flux of one component is affected not only by its own concentration gradient, but also by other driving forces including the concentration gradients of other components and the temperature gradient. However, since the mass fraction of the NPs is very low, the mass-averaged velocity of the system in the absence of bulk motion is essentially the same as that of the solvent (fluid) velocity (Deen Reference Deen1998). This implies that the diffusion of the NPs does not affect the molecular flux of the miscible base fluids. Accordingly, it is legitimate to neglect the NP flux effect on the transport of the fluids (3.3). Furthermore, the fluids molecular flux effect on the transport of the NPs and the temperature gradient effect on the transport of the fluids (Soret effect) are ignored.
It is well known that NPs increase the viscosity, thermal conductivity and the convective heat transfer coefficient of the base fluid (Murshed & Estelle Reference Murshed and Estelle2017; Ahmadi et al.
Reference Ahmadi, Mirlohi, Alhuyi Nazari and Ghasempour2018). Buongiorno (Reference Buongiorno2006) attributed the abnormal increase of the convective heat transfer coefficient to the NP/base fluid relative (slip) velocity. As a result of his analysis of the effect of seven different slip mechanisms, it was concluded that Brownian diffusion and thermophoresis (particle migrations under the effect of temperature gradient) are the most important slip mechanisms. Accordingly, he proposed a two component (NP–fluid) model for the mass and heat transport in nanofluids. Equation (3.4) is based on his model for the transport of NPs and is adapted for flows in porous media (Nield & Kuznetsov Reference Nield and Kuznetsov2009). The first and second terms in the right-hand side of this equation represent the net rate of change of NP concentration as a result of Brownian diffusion and thermophoresis, respectively. In equations (3.3) and (3.4),
$\unicode[STIX]{x1D719}$
is the medium porosity,
$C_{i}$
the mass fraction or volume fraction (in ideal solution),
$T$
the absolute temperature,
$D_{BA}$
the mutual diffusion coefficient of the fluids and
$D_{Br}$
and
$D_{T}$
are the diffusion coefficients of NPs resulting from Brownian motion and thermophoresis. Finally, the energy equation (3.5) is also used to model this non-isothermal system, where the last two terms represent the contributions of the Brownian motion and thermophoresis. In the energy equation,
$\unicode[STIX]{x1D706}$
is the thermal lag coefficient,
$\unicode[STIX]{x1D6FC}$
the thermal diffusivity and
$\unicode[STIX]{x1D6FD}$
the relative volumetric heat capacity of the NPs to that of the medium. Thus
$$\begin{eqnarray}\unicode[STIX]{x1D706}=\frac{\unicode[STIX]{x1D719}(\unicode[STIX]{x1D70C}c_{p})_{nf}}{(\unicode[STIX]{x1D70C}c_{p})_{m}},\quad \unicode[STIX]{x1D6FC}=\frac{k_{m}}{(\unicode[STIX]{x1D70C}c_{p})_{m}},\quad \unicode[STIX]{x1D6FD}=\frac{\unicode[STIX]{x1D719}(\unicode[STIX]{x1D70C}c_{p})_{p}}{(\unicode[STIX]{x1D70C}c_{p})_{m}},\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D70C}c_{p})_{m}=\unicode[STIX]{x1D719}(\unicode[STIX]{x1D70C}c_{p})_{nf}+(1-\unicode[STIX]{x1D719})(\unicode[STIX]{x1D70C}c_{p})_{s}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle k_{m}=\unicode[STIX]{x1D719}k_{nf}+(1-\unicode[STIX]{x1D719})k_{s}. & \displaystyle\end{eqnarray}$$
In the above equations,
$k$
is the thermal conductivity and
$\unicode[STIX]{x1D70C}c_{p}$
the volumetric heat capacity. The subscripts
$m$
,
$p$
,
$s$
and
$nf$
stand for the medium, NPs, solid and nanofluid respectively. In the present model
$\unicode[STIX]{x1D6FC}$
is assumed to be constant since the thermal conductivity of the porous medium is mostly dominated by the solid phase, rather than the NPs. Even in the case of a Hele-Shaw cell, it is legitimate to take a constant value for the thermal conductivity of nanofluid
$(k_{nf})$
as one of its properties. The value of
$(\unicode[STIX]{x1D70C}c_{p})_{nf}$
is also assumed constant by a similar argument. In (3.4),
$k_{dep}$
is the deposition rate described with the widely applied colloidal filtration theory (CFT) (Yao, Habibian & O’Melia Reference Yao, Habibian and O’Melia1971; Tufenkji & Elimelech Reference Tufenkji and Elimelech2004). A list of deposition models along with their pros and cons according to experimental data can be found in Zhang (Reference Zhang2012) and Zhang et al. (Reference Zhang, Murphy, Yu, Huh and Bryant2016).
To complete the formulation, a brief discussion of thermophoresis is presented. Thermophoresis is an additional transport mechanism of particles dispersed in a continuum aside from Brownian diffusion. It is associated with temperature gradients and is the equivalent of the Soret effect in gaseous and liquid mixtures. Physically, it represents the time-averaged chaotic impulse of the surrounding fluid on the particles. It must be noted that due to this random force, particles move randomly in all directions, but their average displacement is slightly biased to the opposite direction of the temperature gradient. From (3.4), the total mass flux of particles aside from convection is
$$\begin{eqnarray}\frac{\boldsymbol{J}_{p}}{\unicode[STIX]{x1D70C}_{p}}=-D_{Br}\unicode[STIX]{x1D735}C_{n}-D_{T}\frac{\unicode[STIX]{x1D735}T}{T}.\end{eqnarray}$$
The Brownian diffusion coefficient can be expressed with the Stokes–Einstein equation,
$D_{Br}=k_{B}T/3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}_{bf}d_{p}$
, where
$k_{B}$
is the Boltzmann constant,
$\unicode[STIX]{x1D707}_{bf}$
the viscosity of the base fluid and
$d_{p}$
the diameter of the particles. On the other hand, although different correlations for thermophoretic velocity and accordingly thermophoretic diffusivity have been introduced in gaseous mixtures, there are rather limited studies when it comes to liquids. A thermophoretic diffusion coefficient based on experimental data reported by McNab & Meisen (Reference McNab and Meisen1973) for relatively large
$(d_{p}\sim 1~\unicode[STIX]{x03BC}\text{m})$
latex spheres particles in water and n-hexane is expressed as follows (Buongiorno Reference Buongiorno2006):
$$\begin{eqnarray}D_{T}=-\unicode[STIX]{x1D6FD}\frac{\unicode[STIX]{x1D707}_{bf}}{\unicode[STIX]{x1D70C}_{bf}}C_{n},\end{eqnarray}$$
where
$\unicode[STIX]{x1D70C}_{bf}$
is the density of the base fluid and
$\unicode[STIX]{x1D6FD}=0.26(k_{bf}/k_{bf}+k_{p})$
. Since (3.9) is derived for very large particles, it would result in a very small thermophoretic diffusion coefficient for NPs. Furthermore, it suggests almost zero
$D_{T}$
for highly conductive NPs
$(k_{p}\gg k_{bf})$
, which is not in accord with experimental measurements involving nanoflows (Martin & Bou-Ali Reference Martin and Bou-Ali2011; Gargiulo et al.
Reference Gargiulo, Cerrota, Cortes, Violi and Stefani2016) and indicates that this model is not fit to represent the thermophoretic diffusivity in nanofluids. Recently, Michaelides (Reference Michaelides2015) proposed a model for the thermophoretic velocity of NPs in liquids, and suggested that it depends strongly on the size of the NPs. The model which is validated with experimental data suggests fairly considerable values of the thermophoretic velocity and accordingly thermophoretic diffusivity of NPs. The thermophoretic diffusion coefficient derived from his model is similar to (3.9) with
$\unicode[STIX]{x1D6FD}=A_{0}(r_{p}/r_{p_{0}})^{-B_{0}}$
, where
$r_{p}$
is the radius of the NPs in nanometres,
$r_{p_{0}}=1~\text{nm}$
, while
$A_{0}$
and
$B_{0}$
are experimental constants that depend on the NP/base fluid system (Michaelides Reference Michaelides2015).
Thermophoretic diffusion of particles is typically expressed in the literature by the thermophoretic mobility
$\hat{D_{T}}=D_{T}/C_{n}T\,[m^{2}/sK]$
or by
$\widetilde{D_{T}}=D_{T}/C_{n}\,[m^{2}/s]$
, known as the true thermophoretic diffusion coefficient (Piazza & Parola Reference Piazza and Parola2008). Although very diverse values can be found in the literature, in most investigated systems the thermophoretic mobility varies within the fairly limited range of
$10^{-12}<\hat{D_{T}}<10^{-11}m^{2}/sK$
(Piazza & Parola Reference Piazza and Parola2008). Accordingly, the value of the true thermophoretic diffusion coefficient is of the order of
$\widetilde{D_{T}}<O(10^{-9})m^{2}/s$
. However, larger values,
$O(10^{-10})m^{2}/sK$
for
$\hat{D_{T}}$
or
$O(10^{-8})m^{2}/s$
for
$\widetilde{D_{T}}$
, have also been reported in some studies dealing with metal NPs (Gargiulo et al.
Reference Gargiulo, Cerrota, Cortes, Violi and Stefani2016). Values of the true thermophoretic diffusion coefficient
$\widetilde{D_{T}}$
for some metal or metal-oxide NPs derived from the correlation introduced by Michaelides (Reference Michaelides2015) can even be
$O(10^{-7})m^{2}/s$
.
Regarding the relative importance of Brownian versus thermophoretic diffusions, small values of
$\unicode[STIX]{x1D6FF}=\widetilde{D_{T}}/D_{Br}$
, as low as zero (or
$\unicode[STIX]{x1D6FF}<0$
, which is not the focus of this study) as well as large ones, e.g.
$\unicode[STIX]{x1D6FF}=156$
(
$R=53~\text{nm}$
) and
$\unicode[STIX]{x1D6FF}=641$
(
$R=253~\text{nm}$
), for polystyrene (PS) spheres in aqueous systems have been reported (Braibanti, Vigolo & Piazza Reference Braibanti, Vigolo and Piazza2008). Accordingly, the effects of thermophoresis can be very important compared to those of Brownian diffusion and therefore cannot be ignored when modelling the transport of NPs.
In this study the Brownian and thermophoretic diffusion coefficients are assumed to be linear functions of
$T$
and
$C_{n}$
respectively, and to depend on the base fluid constant viscosity and density (Buongiorno Reference Buongiorno2006). Accordingly, one may express them as
$D_{Br}=D_{Br_{0}}T/T_{a}$
and
$D_{T}=D_{T_{0}}C_{n}/C_{n0}$
where
$D_{Br_{0}}$
and
$D_{T0}$
are constants. It must be noted that
$D_{T}$
must not be treated as constant, as this may lead to unphysical negative values of the NP concentrations, especially when the NP concentration gradient and temperature gradient are in opposite directions. This is a common misapplication encountered in the literature involving numerical simulations of nanoflows based on Buongiorno’s model, e.g. Zaraki et al. (Reference Zaraki, Ghalambaz, Chamkha, Ghalambaz and De Rossi2015).
The governing equations are made dimensionless and then formulated in a reference frame moving at the injection velocity
$U$
:
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle x^{\ast }/y^{\ast }=\frac{x/y}{\displaystyle \frac{D_{BA}\unicode[STIX]{x1D719}}{U}},\quad \boldsymbol{V}^{\ast }=\frac{\boldsymbol{V}_{D}}{U},\quad t^{\ast }=\frac{t}{\displaystyle \frac{D_{BA}\unicode[STIX]{x1D719}^{2}}{U^{2}}},\quad P^{\ast }=\frac{P}{\bar{\unicode[STIX]{x1D707}}_{a0}D_{BA}\unicode[STIX]{x1D719}},\\ \displaystyle \unicode[STIX]{x1D707}^{\ast }=\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D707}_{a0}},\quad C_{b}^{\ast }=\frac{C_{b}}{C_{b_{0}}},\quad C_{n}^{\ast }=\frac{C_{n}}{C_{n_{0}}},\quad \unicode[STIX]{x1D703}=\frac{T}{T_{a}}.\end{array}\right\}\end{eqnarray}$$
Note that since the permeability
$K$
is constant, it is incorporated in the definition of the viscosity such that, henceforth,
$\unicode[STIX]{x1D707}/K$
stands for
$\unicode[STIX]{x1D707}$
.
The dimensionless equations, where the asterisks have been dropped for convenience, are
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{V}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}P=-\unicode[STIX]{x1D707}(\boldsymbol{V}+\boldsymbol{i}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}C_{b}}{\unicode[STIX]{x2202}t}+\boldsymbol{V}\boldsymbol{\cdot }\unicode[STIX]{x1D735}C_{b}=\unicode[STIX]{x1D6FB}^{2}C_{b}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \!\!\frac{\unicode[STIX]{x2202}C_{n}}{\unicode[STIX]{x2202}t}+\boldsymbol{V}\boldsymbol{\cdot }\unicode[STIX]{x1D735}C_{n}=\unicode[STIX]{x1D6FF}_{B}\unicode[STIX]{x1D703}\unicode[STIX]{x1D6FB}^{2}C_{n}+\left(\!\unicode[STIX]{x1D6FF}_{B}+\frac{\unicode[STIX]{x1D6FF}_{T}}{\unicode[STIX]{x1D703}}\right)\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\boldsymbol{\cdot }\unicode[STIX]{x1D735}C_{n}+\left(\!\unicode[STIX]{x1D6FF}_{T}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D703}}\right)-Dep\!\right)C_{n},\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t}+(\unicode[STIX]{x1D706}-1)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x}+\unicode[STIX]{x1D706}(\boldsymbol{V}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703})=Le\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}+\unicode[STIX]{x1D716}\left(\!\unicode[STIX]{x1D6FF}_{B}\unicode[STIX]{x1D703}\unicode[STIX]{x1D735}C_{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}+\unicode[STIX]{x1D6FF}_{T}C_{n}\frac{\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}}{\unicode[STIX]{x1D703}}\right),\qquad & \displaystyle\end{eqnarray}$$
where
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D6FF}_{B}=\frac{D_{Br_{0}}}{D_{BA}},\quad \unicode[STIX]{x1D6FF}_{T}=\frac{D_{T_{0}}}{D_{BA}C_{n_{0}}},\quad Le=\frac{\unicode[STIX]{x1D6FC}}{D_{BA}},\\ \displaystyle Dep=\frac{k_{dep}D_{BA}\unicode[STIX]{x1D719}^{2}}{U^{2}},\quad \unicode[STIX]{x1D716}=\frac{C_{n_{0}}\unicode[STIX]{x1D719}(\unicode[STIX]{x1D70C}c_{p})_{p}}{(\unicode[STIX]{x1D70C}c_{p})_{m}}.\end{array}\right\}\end{eqnarray}$$
In the above equations
$Le$
is the Lewis number,
$Dep$
the deposition rate coefficient, while
$\unicode[STIX]{x1D6FF}_{B}$
and
$\unicode[STIX]{x1D6FF}_{T}$
are referred to as the Brownian and the thermophoretic diffusivity, respectively. It should be noted that the present study is limited to displacements with no NP deposition, and henceforth
$Dep=0$
. This is a reasonable assumption as NP deposition in porous media can indeed be avoided through the use of either surfactants or surface charge technology (Yu & Xie Reference Yu and Xie2012).
Dimensionless geometric quantities result from the domain boundaries which are
$(-Pe/2,Pe/2)$
in the
$x$
-direction and
$(0,Pe/A_{s})$
in the
$y$
-direction, where
$A_{s}=L/W$
is the domain aspect ratio and
$Pe=UL/(\unicode[STIX]{x1D719}D_{BA})$
is the Péclet number. The boundary conditions in this reference frame consist of zero flux for concentrations and temperature in the
$x$
-direction, and periodicity in the
$y$
-direction:
$$\begin{eqnarray}\displaystyle & \displaystyle (u,v)\left(-\frac{Pe}{2},y,t\right)=(u,v)\left(\frac{Pe}{2},y,t\right)=(0,0), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\unicode[STIX]{x2202}C_{b}}{\unicode[STIX]{x2202}x},\frac{\unicode[STIX]{x2202}C_{n}}{\unicode[STIX]{x2202}x},\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x}\right)\left(-\frac{Pe}{2},y,t\right)=\left(\frac{\unicode[STIX]{x2202}C_{b}}{\unicode[STIX]{x2202}x},\frac{\unicode[STIX]{x2202}C_{n}}{\unicode[STIX]{x2202}x},\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x}\right)\left(\frac{Pe}{2},y,t\right)=(0,0,0),\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle (u,v,C_{b},C_{n},\unicode[STIX]{x1D703})(x,0,t)=(u,v,C_{b},C_{n},\unicode[STIX]{x1D703})\left(x,\frac{Pe}{A_{s}},t\right). & \displaystyle\end{eqnarray}$$
$x$
-direction requires that
$\unicode[STIX]{x1D6FF}_{B}\unicode[STIX]{x1D703}(\unicode[STIX]{x2202}C_{n}/\unicode[STIX]{x2202}x)+(\unicode[STIX]{x1D6FF}_{T}C_{n}/\unicode[STIX]{x1D703})(\unicode[STIX]{x2202}\unicode[STIX]{x1D703}/\unicode[STIX]{x2202}x)$
must be zero, which is automatically satisfied by the conditions
$\unicode[STIX]{x2202}C_{n}/\unicode[STIX]{x2202}x=\unicode[STIX]{x2202}\unicode[STIX]{x1D703}/\unicode[STIX]{x2202}x=0$
.To complete the model a form for the dependence of the viscosity on concentration and temperature must be specified. Following previous studies, an exponential viscosity–concentration–temperature relationship is adopted (Tan & Homsy Reference Tan and Homsy1986; Pritchard Reference Pritchard2009):
$$\begin{eqnarray}\unicode[STIX]{x1D707}=\exp \left(R_{b}C_{b}+R_{\unicode[STIX]{x1D703}}\frac{1-\unicode[STIX]{x1D703}}{1-r}+R_{n}C_{n}\right),\end{eqnarray}$$
where
$r=T_{b}/T_{a}$
is the temperature ratio of the displaced fluid to the displacing one, while
$R_{b}$
,
$R_{\unicode[STIX]{x1D703}}$
and
$R_{n}$
are mobility ratios defined as
$$\begin{eqnarray}R_{b}=\ln \left(\frac{\unicode[STIX]{x1D707}_{b_{0}}}{\unicode[STIX]{x1D707}_{a_{0}}}\right)_{T_{a}},\quad R_{\unicode[STIX]{x1D703}}=\ln \left(\frac{\unicode[STIX]{x1D707}_{T_{b}}}{\unicode[STIX]{x1D707}_{T_{a}}}\right),\quad R_{n}=\ln \left(\frac{\unicode[STIX]{x1D707}_{n_{0}}}{\unicode[STIX]{x1D707}_{a_{0}}}\right)_{T_{a}}.\end{eqnarray}$$
In particular,
$R_{b}>0(R_{b}<0)$
,
$R_{\unicode[STIX]{x1D703}}=R_{n}=0$
represents an unstable (stable) isothermal NP-free system where a low (high) viscosity fluid displaces a high (low) viscosity one, while
$R_{\unicode[STIX]{x1D703}}>0(R_{\unicode[STIX]{x1D703}}<0)$
,
$R_{b}=R_{n}=0$
represents an unstable (stable) non-isothermal NP-free system where a hot (cold) fluid displaces a cold (hot) fluid. Flows with
$R_{n}>0$
represent the case where addition of NPs increases the viscosity of the displacing base fluid. Note that
$R_{\unicode[STIX]{x1D703}}$
and
$r$
enter the model as a group:
$R_{\unicode[STIX]{x1D703}}^{\prime }=R_{\unicode[STIX]{x1D703}}/1-r$
.
4 Numerical method
By taking the curl of Darcy’s law and expressing the velocities in terms of the streamfunction, the pressure is eliminated from the equations and the continuity equation is automatically satisfied. Then by using the expression for the viscosity (3.18), the governing equations are formulated in terms of the vorticity
$(\unicode[STIX]{x1D714})$
, streamfunction
$(\unicode[STIX]{x1D713})$
, concentrations
$(C_{i})$
and temperature
$(\unicode[STIX]{x1D703})$
. The resulting equations are solved numerically using the Hartley based pseudo-spectral method (Zimmerman & Homsy Reference Zimmerman and Homsy1991). To implement this method the boundary conditions (BCs) must be periodic. This condition is automatically satisfied in the transverse direction. Furthermore, since the equations are formulated in the moving reference frame,
$\unicode[STIX]{x1D714}$
and
$\unicode[STIX]{x1D713}$
are zero at the streamwise boundaries and accordingly periodic conditions can be implemented. In order to impose the longitudinal periodic BCs for concentrations and temperature, the model equations are formulated in terms of a base state
$(\bar{C_{i}}(x,t),\bar{\unicode[STIX]{x1D703}}(x,t))$
, that satisfies the BCs, and perturbations
$(C_{i}^{\prime }(x,y,t),\unicode[STIX]{x1D703}^{\prime }(x,y,t))$
which decay to zero at the domain’s streamwise boundaries. One can therefore solve the equations for the concentration and temperature perturbations, which then are added to the base state values to obtain the total
$C_{i}(x,y,t)$
and
$\unicode[STIX]{x1D703}(x,y,t)$
.
The base sate concentrations and temperature are derived from
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\bar{C}_{b}}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x2202}^{2}\bar{C}_{b}}{\unicode[STIX]{x2202}x^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\bar{C}_{n}}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D6FF}_{B}\bar{\unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}^{2}\bar{C}_{n}}{\unicode[STIX]{x2202}x^{2}}+\left(\unicode[STIX]{x1D6FF}_{B}+\frac{\unicode[STIX]{x1D6FF}_{T}}{\bar{\unicode[STIX]{x1D703}}}\right)\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}\bar{C}_{n}}{\unicode[STIX]{x2202}x}+\left(\unicode[STIX]{x1D6FF}_{T}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\frac{1}{\bar{\unicode[STIX]{x1D703}}}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}x}\right)-Dep\right)\bar{C}_{n}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}t}+(\unicode[STIX]{x1D706}-1)\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}x}=Le\frac{\unicode[STIX]{x2202}^{2}\bar{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}x^{2}}+\unicode[STIX]{x1D716}\left(\unicode[STIX]{x1D6FF}_{B}\bar{\unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}\bar{C}_{n}}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x1D6FF}_{T}\bar{C}_{n}}{\bar{\unicode[STIX]{x1D703}}}\left(\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}x}\right)^{2}\right). & \displaystyle\end{eqnarray}$$
It should be mentioned that the base state velocities are zero in the moving reference frame. Zero flux boundary conditions are used to solve the above-mentioned base state equations. The solution of (4.1) is determined analytically as
$$\begin{eqnarray}\bar{C_{b}}(x,t)=\frac{1}{2}\text{erfc}\left(\frac{-x}{2\sqrt{t}}\right).\end{eqnarray}$$
Furthermore, the coupled (4.2) and (4.3) are solved numerically with the method of lines.
The equations for the perturbation terms are
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}=-\unicode[STIX]{x1D714}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D714}=R_{b}N_{b}-\frac{R_{\unicode[STIX]{x1D703}}}{1-r}N_{\unicode[STIX]{x1D703}}+R_{n}N_{n}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}C_{b}^{\prime }}{\unicode[STIX]{x2202}t}=J_{b}+\frac{\unicode[STIX]{x2202}^{2}C_{b}^{\prime }}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}C_{b}^{\prime }}{\unicode[STIX]{x2202}y^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}C_{n}^{\prime }}{\unicode[STIX]{x2202}t}=J_{n}+G_{n}-DepC_{n}^{\prime }, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{\prime }}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D706}J_{\unicode[STIX]{x1D703}}+(1-\unicode[STIX]{x1D706})\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{\prime }}{\unicode[STIX]{x2202}x}+\unicode[STIX]{x1D716}G_{\unicode[STIX]{x1D703}}+Le\left(\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D703}^{\prime }}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D703}^{\prime }}{\unicode[STIX]{x2202}y^{2}}\right), & \displaystyle\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle & \displaystyle N_{i}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}x}\left(\frac{\unicode[STIX]{x2202}\bar{X_{i}}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}X_{i}^{\prime }}{\unicode[STIX]{x2202}x}\right)+\left(1+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}y}\right)\frac{\unicode[STIX]{x2202}X_{i}^{\prime }}{\unicode[STIX]{x2202}y}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle J_{i}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}X_{i}^{\prime }}{\unicode[STIX]{x2202}y}-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}y}\left(\frac{\unicode[STIX]{x2202}\bar{X_{i}}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}X_{i}^{\prime }}{\unicode[STIX]{x2202}x}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle G_{n} & = & \displaystyle \unicode[STIX]{x1D6FF}_{B}\left(\unicode[STIX]{x1D703}^{\prime }\frac{\unicode[STIX]{x2202}^{2}\bar{C_{n}}}{\unicode[STIX]{x2202}x^{2}}+(\bar{\unicode[STIX]{x1D703}}+\unicode[STIX]{x1D703}^{\prime })\left(\frac{\unicode[STIX]{x2202}^{2}C_{n}^{\prime }}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}C_{n}^{\prime }}{\unicode[STIX]{x2202}y^{2}}\right)\right)-\frac{\unicode[STIX]{x1D6FF}_{T}\unicode[STIX]{x1D703}^{\prime }}{\bar{\unicode[STIX]{x1D703}}(\bar{\unicode[STIX]{x1D703}}+\unicode[STIX]{x1D703}^{\prime })}\frac{\unicode[STIX]{x2202}\bar{C}_{n}}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}x}\nonumber\\ \displaystyle & & \displaystyle +\left(\unicode[STIX]{x1D6FF}_{B}+\frac{\unicode[STIX]{x1D6FF}_{T}}{\bar{\unicode[STIX]{x1D703}}+\unicode[STIX]{x1D703}^{\prime }}\right)\left(\frac{\unicode[STIX]{x2202}\bar{C_{n}}}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{\prime }}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}C_{n}^{\prime }}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}C_{n}^{\prime }}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{\prime }}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}C_{n}^{\prime }}{\unicode[STIX]{x2202}y}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{\prime }}{\unicode[STIX]{x2202}y}\right)\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FF}_{T}(\bar{C_{n}}+C_{n}^{\prime })\!\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\frac{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{\prime }}{\unicode[STIX]{x2202}x}+\displaystyle \frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}x}}{\bar{\unicode[STIX]{x1D703}}+\unicode[STIX]{x1D703}^{\prime }}\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\left(\frac{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{\prime }}{\unicode[STIX]{x2202}y}}{\bar{\unicode[STIX]{x1D703}}+\unicode[STIX]{x1D703}^{\prime }}\right)\right]\!-\unicode[STIX]{x1D6FF}_{T}\bar{C_{n}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\frac{\displaystyle \frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}x}}{\bar{\unicode[STIX]{x1D703}}}\right)\!,\qquad\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle G_{\unicode[STIX]{x1D703}} & = & \displaystyle \unicode[STIX]{x1D6FF}_{B}\unicode[STIX]{x1D703}^{\prime }\frac{\unicode[STIX]{x2202}\bar{C_{n}}}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x1D6FF}_{T}}{\bar{\unicode[STIX]{x1D703}}+\unicode[STIX]{x1D703}^{\prime }}\left(\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}x}\right)^{2}\left(C_{n}^{\prime }-\frac{\bar{C_{n}}\unicode[STIX]{x1D703}^{\prime }}{\bar{\unicode[STIX]{x1D703}}}\right)\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\unicode[STIX]{x1D6FF}_{T}(\bar{C_{n}}+C_{n}^{\prime })}{\bar{\unicode[STIX]{x1D703}}+\unicode[STIX]{x1D703}^{\prime }}\left(\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{\prime }}{\unicode[STIX]{x2202}x}\right)^{2}+\left(\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{\prime }}{\unicode[STIX]{x2202}y}\right)^{2}+2\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{\prime }}{\unicode[STIX]{x2202}x}\right)\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FF}_{B}(\bar{\unicode[STIX]{x1D703}}+\unicode[STIX]{x1D703}^{\prime })\left(\frac{\unicode[STIX]{x2202}\bar{C_{n}}}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{\prime }}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}C_{n}^{\prime }}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}C_{n}^{\prime }}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{\prime }}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}C_{n}^{\prime }}{\unicode[STIX]{x2202}y}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{\prime }}{\unicode[STIX]{x2202}y}\right),\end{eqnarray}$$
where
$X_{i}=(C_{i},\unicode[STIX]{x1D703})$
. As discussed above,
$\unicode[STIX]{x1D714}(x,y,t)$
,
$\unicode[STIX]{x1D713}(x,y,t)$
,
$C_{i}^{\prime }(x,y,t)$
and
$\unicode[STIX]{x1D703}^{\prime }(x,y,t)$
are zero at the streamwise boundaries and according to (3.17c
) they are periodic in the transverse direction. Moreover, zero values for the vorticity and streamfunction, and a random number distribution for the concentration and the temperature perturbations are used as initial conditions:
$$\begin{eqnarray}(C_{i}^{\prime },\unicode[STIX]{x1D703}^{\prime })(x,y,t_{0})=\unicode[STIX]{x1D6FF}\;\text{rand}(y)\exp \left(-\frac{x^{2}}{\unicode[STIX]{x1D70E}^{2}}\right).\end{eqnarray}$$
Where
$\text{rand}(y)$
is a random number between
$(-1,1)$
,
$\unicode[STIX]{x1D6FF}$
, a small number relative to unity, is the perturbation magnitude, and
$\unicode[STIX]{x1D70E}^{2}$
, also a very small number, represents the penetration of the disturbances in the domain. After transforming the equations to the Hartley space, the resulting ordinary differential equations are stepped in time using a semi-implicit algorithm based on the Adams–Bashforth and Adams–Moulton predictor–corrector scheme. One can find more details about this procedure in Islam & Azaiez (Reference Islam and Azaiez2005). Finally, in order to validate the code, the numerical convergence was checked, and it was found that a grid size of
$256\times 256$
is satisfactory. Furthermore, results were also validated by comparing with isothermal NP-laden and NP-free displacements. LSA was also carried out and the instability characteristic equation was derived for an arbitrary choice of parameters in an NP-laden non-isothermal system. The number of fingers predicted from the LSA was compared with that of the early growing fingers from the NLS. From the characteristic curve in figure 2, the critical wavenumber is
$k_{cr}\approx 0.2$
, resulting in a number of fingers
$NF=kPe/2A_{s}\unicode[STIX]{x03C0}\approx 16.3$
. This is in agreement with the results of NLS at early times where
$NF\approx 17$
. This further confirms the validity of the nonlinear code.
(a) Instability characteristic curve from LSA; (b) iso-contour of
$C_{b}=0.5$
at early times from NLS.
$R_{b}=3$
,
$R_{\unicode[STIX]{x1D703}}=1$
,
$R_{n}=1$
,
$\unicode[STIX]{x1D706}=0.8$
,
$\unicode[STIX]{x1D6FF}_{B}=1$
,
$\unicode[STIX]{x1D6FF}_{T}=3$
,
$Le=3$
,
$\unicode[STIX]{x1D716}=0.01$
,
$Dep=0$
,
$r=0.5$
,
$t_{0}=1$
,
$Pe=1024$
,
$A_{s}=2$
.
5 Results and discussion
Before presenting the results and discussion, a few adopted terminologies and reference parameters are defined. Throughout this study, an intrinsically unstable (stable) system corresponds to a system without NPs that is unstable (stable). This is to be contrasted with unstable (stable) systems that refer to unstable (stable) systems in the presence of NPs with or without heat transfer. Furthermore, this study focuses on NPs that increase the viscosity of the base fluid
$(R_{n}>0)$
.
Two types of non-isothermal displacements are considered: a hot fluid displacing a cold one (HDC) where
$R_{\unicode[STIX]{x1D703}}>0$
and
$0<r<1$
, and a cold fluid displacing a hot one (CDH) where
$R_{\unicode[STIX]{x1D703}}<0,r>1$
. The reference parameters for the HDC system are
$R_{\unicode[STIX]{x1D703}}=1,r=0.5$
, while they are
$R_{\unicode[STIX]{x1D703}}=-1,r=2$
for the CDH system. Furthermore, unless mentioned otherwise, the other parameters are set as
$R_{b}=3$
,
$R_{n}=1,\unicode[STIX]{x1D706}=1$
,
$\unicode[STIX]{x1D6FF}_{B}=1$
,
$\unicode[STIX]{x1D6FF}_{T}=3$
,
$Le=3$
,
$\unicode[STIX]{x1D716}=0.001$
,
$Dep=0$
,
$Pe=1024$
,
$A_{s}=2$
and
$t_{0}=1$
. The choice of these parameters represents unstable systems.
The qualitative characterization of the instability will be supported by a quantitative analysis based on transverse-averaged concentrations
$(C_{av})$
and the mixing length (ML). The mixing length is defined as the ratio of the length within concentration range
$C_{b,av}=0.01$
to
$C_{b,av}=0.99$
to the length of the whole domain, where
$$\begin{eqnarray}C_{b,av}(x,t)=\frac{A_{s}}{Pe}\int _{0}^{Pe/A_{s}}C_{b}(x,y,t)\,\text{d}y.\end{eqnarray}$$
Larger values of mixing length indicate a more unstable situation and vice versa.
In our previous related study involving isothermal systems (Dastvareh & Azaiez Reference Dastvareh and Azaiez2017) it was found that the addition of NPs to the displacing fluid
$(R_{n}>0)$
attenuates the instability of an intrinsically unstable solution as it decreases the viscosity contrast in the system. More specifically, it was found that in the absence of deposition, the system is unstable if
$R_{b}-R_{n}>0$
and is stable otherwise. However, in the presence of deposition, NPs are gradually removed from the system and their effects on the viscosity of the fluids are weakened. As a result the instability is determined mainly by the viscosity ratio of the base fluids. In other words, if
$R_{b}>0$
the system is unstable, although the addition of NPs mitigates the instability. It implies that NPs
$(R_{n}>0)$
cannot make an intrinsically unstable system stable in an isothermal case in the presence of deposition. It was also reported that Brownian diffusion increases the instability of an unstable system. However, it cannot make either stable or intrinsically stable systems, unstable. Unlike the isothermal flow, Brownian diffusion coefficient of the NPs
$(D_{Br})$
here, varies with temperature. We will therefore start the discussion by analysing its effects under non-isothermal conditions.
Contours of
$C_{b}$
for two values of Brownian diffusivity in both HDC and CDH systems.
5.1 Brownian diffusion
Figure 3 depicts concentration distributions of species B for the HDC and CDH case for two values of the Brownian diffusivity. Similar to the isothermal system, as
$\unicode[STIX]{x1D6FF}_{B}$
increases, both HDC and CDH systems exhibit more intricate finger configurations, reflecting the destabilizing effect of Brownian diffusion. Plots of the variation of the mixing length with time for different values of
$\unicode[STIX]{x1D6FF}_{B}$
are presented for both systems in figure 4. The trends of a systematic increase of the mixing length with Brownian diffusion confirm the previous conclusion quantitatively. One may also note that the effects of
$\unicode[STIX]{x1D6FF}_{B}$
are more prominent in the CDH system compared to the HDC one, and that the concentration distribution is more diffuse in the latter case. This points to an important role that the temperature probably plays in conjunction with Brownian diffusion in determining the distribution of NPs in the flow, and in turn in affecting the instability. In order to analyse these effects, plots of the transversely averaged concentration of the NPs
$(C_{n,av})$
are presented in figure 5 with an inset of NP concentration contours at an instant in time. An examination of the concentration contours in the insets reveals that the distributions of the NPs are different in the HDC and CDH systems, and that in the latter, changes in
$\unicode[STIX]{x1D6FF}_{B}$
induce a dramatic variation of the distribution at the front. In particular, NPs in the HDC system move downstream easily in the low resistance paths provided by the growing fingers. On the other hand, in the CDH system, there is a resistance to the downstream transport of the nanoparticles, which latter tend to accumulate at the front. For the small Brownian diffusivity, the NP concentration is characterized by a narrow band at the front with very large concentration values. This local accumulation of NPs leads to a local increase of the nanofluid viscosity which reduces the mobility ratio between the fluids and results in the observed attenuation of the instability. For
$\unicode[STIX]{x1D6FF}_{B}=1$
, however, the distribution is more spread out around the front with a wide region of small NP concentrations that contributes to the growth of the instability. These conclusions about the distribution of the NP concentrations are confirmed by the corresponding plots of
$C_{n,av}$
. Furthermore, it is worth noting from figure 5 that the maximum average NP accumulation in the CDH system starts to decrease after some time, but its location is systematically receding from the front and tends to move upstream.
Variation of the mixing length with time for different Brownian diffusivities: (a) HDC system and (b) CDH system.
Figure 5 also reveals the development of localized regions of high NP concentrations in the CDH case for
$\unicode[STIX]{x1D6FF}_{B}=1$
, where NPs in the upstream tend to accumulate in the concave regions between the backward-growing fingers (see figure 5(b) in the CDH case). A close examination of the corresponding streamlines presented in figure 6 reveals that as a result of the convective flow where counter-rotating dipoles develop, fingers grow backward and at the same time NPs are entrained into the space between the growing fingers. This results in the development of regions with large localized concentrations of NPs.
Variation of the transversely averaged NP concentration along with contours of
$C_{n}$
for both HDC and CDH systems: (a,c)
$\unicode[STIX]{x1D6FF}_{B}=0.1$
and (b,d)
$\unicode[STIX]{x1D6FF}_{B}=1.0$
.
Accumulation of the NPs in the region between the backward growing fingers as a result of the counter-spinning dipoles in the CDH system
$(\unicode[STIX]{x1D6FF}_{B}=1)$
.
The intensive accumulation of NPs in the narrow band at the front in the CDH system with
$\unicode[STIX]{x1D6FF}_{B}=0.1$
leads to the creation of stable regions (negative viscosity gradient) upstream of the front. The viscosity distribution depicted in figure 7 clearly illustrates the local regions of stable flow where a large viscosity fluid displaces a low viscosity one. The development of local stable spots leads to a situation where the instabilities cannot grow easily as already observed in figure 3(a).
To further investigate this, the viscosity distribution is presented along with the streamlines in figure 8. Four representative regions with specific features have been identified and are labelled in the figure. A relatively high viscosity region 1, which is the consequence of the accumulation of NPs, and region 2 sandwich a lower viscosity region 3. Due to this configuration, the less viscous (more mobile) fluid in region 3 experiences resistance from region 2 to flow downstream and from region 1 to flow upstream. As a result, region 2 acts as a barrier to the flow while the identified region 4 offers a low resistance path to the flow. Due to this barrier, the fluid in region 1 cannot easily flow downstream to feed both forward and backward fingers. The fluid in region 2 experiences the same resistance to flow upstream resulting in the suppression of backward extending fingers. Although part of the fluid can flow downstream through the low resistance region 4, the forward fingers are also unable to grow, since first the extent of the low resistance paths is very small as the stable spots are dispersed in the channel width, and second the mobility ratio between the fluids is reduced because of the presence of small diffusing NPs. This explains why the initiated instabilities are hardly able to grow in the case
$\unicode[STIX]{x1D6FF}_{B}=0.1$
. The larger value
$\unicode[STIX]{x1D6FF}_{B}=1$
results in the weakening or even disappearance of the upstream stable spots where the fluid does not experience resistance to flow downstream or to backward growth. Furthermore, due to the diffusion of NPs, the mobility ratio between the fluids increases, resulting in a relatively easier growth of forward fingers. Accordingly, finger growth is stronger compared to the case
$\unicode[STIX]{x1D6FF}_{B}=0.1$
.
Viscosity contours for two values of Brownian diffusivity in the CDH system: (a)
$\unicode[STIX]{x1D6FF}_{B}=0.1$
and (b)
$\unicode[STIX]{x1D6FF}_{B}=1.0$
.
Viscosity contours and streamlines in the CDH system for
$\unicode[STIX]{x1D6FF}_{B}=0.1$
.
Figure 9 depicts the effects of the temperature ratio
$r$
on the NP distribution and the corresponding variation of
$C_{n,av}$
along the channel. It is clear that increasing the temperature difference of the fluids at constant Brownian diffusivity leads to more downstream migration of NPs in the HDC system and higher upstream accumulation of the NPs in the CDH system, similar to the reported effects of
$\unicode[STIX]{x1D6FF}_{B}$
. Clearly temperature gradients in NP-laden systems have dramatic effects on the instability. One may therefore infer that Brownian diffusion is not the only mechanism responsible for the observed trends, and heat related mechanisms may be important. Since these trends were not reported in NP-free flows, it is natural to suspect that thermophoresis, which is responsible for the transport of NPs as a result of temperature gradients, is playing a major role in the dynamics of the flow.
Variation of the transversely averaged NP concentration with insets of NP concentration contours at different temperature ratios for (a) the HDC system
$(t=500)$
and (b) the CDH system
$(t=1500)$
.
5.2 Thermophoresis
Thermophoretic effects are examined by analysing results for different values of the thermophoretic diffusivity
$\unicode[STIX]{x1D6FF}_{T}$
at constant temperature ratio
$r$
. It should be noted that as
$\unicode[STIX]{x1D6FF}_{T}$
increases, the thermophoretic velocity of the NPs increases and vice versa.
Figure 10 shows contours of
$C_{b}$
for different values of the thermophoretic diffusivity in the HDC system. It is observed that higher values of
$\unicode[STIX]{x1D6FF}_{T}$
lead to a more unstable flow. The quantitative analysis also validates this conclusion, where for instance, according to figure 11, the system with higher values of
$\unicode[STIX]{x1D6FF}_{T}$
has a larger mixing length. This increase in the instability corresponds to the local enhancement in the viscosity contrast at the interface as a result of the migration of NPs downstream at higher
$\unicode[STIX]{x1D6FF}_{T}$
.
Figure 11 reveals that as
$\unicode[STIX]{x1D6FF}_{T}$
increases, the shift from the dispersion regime to the nonlinear fingering regime occurs at earlier times. In the nonlinear fingering regime, the fingers interact closely with each other and the ML grows linearly with time. Furthermore, for the set of parameters considered, the rate of growth of the ML is virtually independent of
$\unicode[STIX]{x1D6FF}_{T}$
. One must note that this linear growth with time in the convective regime is typically observed in constant viscosity ratios and large Péclet number flows where the ML growth rate does not depend on dispersion (Tan & Homsy Reference Tan and Homsy1986; Zimmerman & Homsy Reference Zimmerman and Homsy1991; Malhotra, Sharma & Lehman Reference Malhotra, Sharma and Lehman2015; Nijjer, Hewitt & Neufeld Reference Nijjer, Hewitt and Neufeld2018). However, since thermophoresis is not a pure diffusive phenomenon, a close analysis of the flow revealed another important factor that can affect the ML growth rate in the nonlinear regime, namely the development of non-monotonic viscosity distributions due to the NP migration. In particular we found that the ML growth deviates from the linear regime and also becomes dependent on the thermophoretic diffusivity for large enough NP viscosity ratios,
$R_{n}$
, as long as the Péclet number is not very large. This effect, which is due to the reverse fingering, may be observed in both HDC and CDH systems but it is more prominent in the latter case.
An examination of contours of NP concentration depicted in figure 12 reveals that, at large enough values of the thermophoretic diffusivity, localized regions with accumulation of NPs can actually develop in the flow which move upstream with the backward growing fingers (see result for
$\unicode[STIX]{x1D6FF}_{T}=10$
). This trend was found to be stronger for even larger values of the thermophoretic diffusivity and the nanofluid viscosity. This unexpected trend will be discussed further below.
Contours of
$C_{b}$
for different thermophoretic diffusivities in the HDC system.
Variation with time of the mixing length for different thermophoretic diffusivities in the HDC system.
Contours of
$C_{n}$
for different thermophoretic diffusivities in the HDC system.
Different trends were obtained in the CDH case. The concentration contours of figure 13 show that, unlike the HDC system, increasing
$\unicode[STIX]{x1D6FF}_{T}$
tends now to attenuate the instability. In particular backward-growing fingers that tend to extend upstream in the absence of thermophoretic effects are noticeably inhibited for
$\unicode[STIX]{x1D6FF}_{T}=10$
. To understand the reason for this attenuation, the distributions of the NP concentration and of the viscosity are analysed (figure 14). It is clear that increasing
$\unicode[STIX]{x1D6FF}_{T}$
leads to the accumulation of NPs in the upstream part of the flow, such that the larger the value of
$\unicode[STIX]{x1D6FF}_{T}$
the more the NP accumulation. The upstream NP accumulation also affects the presence of NPs downstream. Specifically, as
$\unicode[STIX]{x1D6FF}_{T}$
increases, the larger accumulation of NPs in the upstream part of the flow implies lower downstream concentration. At higher values of
$\unicode[STIX]{x1D6FF}_{T}$
, where the accumulation of NPs is more intense, a local increase in the viscosity is observed. The new distribution of NPs leads to the creation of locally stable upstream regions and hampers the growth of backward developing fingers. This stable position is clearly identified in the viscosity distribution for
$\unicode[STIX]{x1D6FF}_{T}=10$
(figure 14). Further analysis of the transverse-averaged NP concentration
$(C_{n,av})$
, which is not illustrated here, indicates that as time passes, the location of
$(C_{n,av})_{max}$
is receding from the front and its value starts to decreases after some time, similar to our discussion of this trend in figure 5 for the CDH system.
Contours of
$C_{b}$
for different thermophoretic diffusivities in the CDH system.
Contours of
$C_{n}$
and viscosity for different thermophoretic diffusivity in the CDH system.
Further increase of
$\unicode[STIX]{x1D6FF}_{T}$
results in another trend in terms of the forward-growing fingers. Figure 15 shows the contours of
$C_{b}$
,
$C_{n}$
and of the viscosity for
$\unicode[STIX]{x1D6FF}_{T}=20$
. As expected, the NP accumulation is intensified and the backward growth of fingers is more obstructed. However, counter-intuitively, forward developing fingers are actually growing faster, ultimately resulting in the increase of the mixing length (see figure 16). This reveals two distinct regimes, whereby initially the mixing length increases with increasing
$\unicode[STIX]{x1D6FF}_{T}$
, followed by a reverse of the trend beyond a critical value of
$\unicode[STIX]{x1D6FF}_{T}\approx 10$
. It must however be noted that the conclusions regarding higher accumulations of NPs at the front and the creation of strong stable regions leading to the suppression of backward instabilities, are still valid for higher
$\unicode[STIX]{x1D6FF}_{T}$
. Furthermore, the results indicate that, generally, the CDH system is less unstable in the presence of thermophoresis in comparison to the case where thermophoresis is absent
$(\unicode[STIX]{x1D6FF}_{T}=0)$
.
Contours of concentrations and viscosity in the CDH system for
$\unicode[STIX]{x1D6FF}_{T}=20$
.
Variation with time of the mixing length for different thermophoretic diffusivities in the CDH system.
To understand the tendency for a faster downstream growth of fingers at higher
$\unicode[STIX]{x1D6FF}_{T}$
, viscosity contours for
$\unicode[STIX]{x1D6FF}_{T}=20$
along with the streamlines are presented in figure 17. First, it must be mentioned that the reverse flow observed in all miscible displacements is driven by convection rather than dispersion, such that the fluid failing to flow through the low resistance path created by the forward fingers tends to flow backward (Manickam & Homsy Reference Manickam and Homsy1994). Via the same arguments used in the interpretation of the results in figure 8, at large values of
$\unicode[STIX]{x1D6FF}_{T}$
the backward fingers encounter a region of strong resistance to growth, that arises from the local accumulation of NPs. As a result, the fluid tends instead to return to flow downstream and contributes to the further growth of forward-developing fingers. The streaming of the fluid to the forward extending finger is coupled with another effect. As one can notice from figure 15, the concentration of NPs downstream is reduced as a result of their larger accumulation at the front. The ensuing depletion of NPs, especially inside the forward developing fingers, leads to an increase of the viscosity contrast inside and outside these fingers. These coupled effects explain the observed faster growth of forward-developing fingers and lead to the larger mixing length.
Viscosity contours and streamlines in the CDH system for
$\unicode[STIX]{x1D6FF}_{T}=20$
.
We will now focus on analysing the mechanisms of NP accumulation in both HDC and CDH systems. First, the distributions of the base state characteristics in the HDC case, where an unexpected trend was reported, are examined. Figure 18 depicts the base state concentration of NPs at different values of
$\unicode[STIX]{x1D6FF}_{T}$
along with a representative variation of
$\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x$
and
$\unicode[STIX]{x2202}^{2}\ln \bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x^{2}$
. In this figure,
$\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x$
varies from zero to negative values with its absolute maximum at the centre of the channel. Furthermore,
$\unicode[STIX]{x2202}^{2}\ln \bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x^{2}=0$
at both ends, with a local positive maximum value at
$x>0$
and a local negative minimum value at
$x<0$
. It is observed that the NP front advances in the channel as
$\unicode[STIX]{x1D6FF}_{T}$
increases. In addition, for small
$\unicode[STIX]{x1D6FF}_{T}$
, NP concentration decreases monotonically along the channel. However, as
$\unicode[STIX]{x1D6FF}_{T}$
increases, an inflection point appears in the profile of
$\bar{C_{n}}$
after a critical value of the thermophoretic diffusivity
$(\unicode[STIX]{x1D6FF}_{T_{cr}})$
. Further increase of
$\unicode[STIX]{x1D6FF}_{T}$
results in a non-monotonic concentration profile, exhibiting a local minimum and maximum. The local increase in the profile of
$\bar{C_{n}}$
represents the local accumulation of the NPs in the flow. The critical values of
$\unicode[STIX]{x1D6FF}_{T}$
associated with the development of an inflection point in the profile of
$\bar{C_{n}}$
are generally a function of all other parameters and time. However, as a general trend it was found that with all other parameters fixed,
$\unicode[STIX]{x1D6FF}_{T_{cr}}$
increases with increasing Brownian diffusivity
$\unicode[STIX]{x1D6FF}_{B}$
. This fact is illustrated in a representative plot in figure 19.
Base state NP concentration and a scaled variation of
$\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x$
and
$\unicode[STIX]{x2202}^{2}\ln \bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x^{2}$
along the channel in the HDC system at
$t=500$
.
A typical variation of
$\unicode[STIX]{x1D6FF}_{T_{cr}}$
with
$\unicode[STIX]{x1D6FF}_{B}$
in the HDC system at
$t=500$
and other default parameters.
These trends can be understood by looking at the equation for the NP base state concentration:
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\bar{C}_{n}}{\unicode[STIX]{x2202}t}=\unicode[STIX]{x1D6FF}_{B}\bar{\unicode[STIX]{x1D703}}\frac{\unicode[STIX]{x2202}^{2}\bar{C}_{n}}{\unicode[STIX]{x2202}x^{2}}+\unicode[STIX]{x1D6FF}_{T}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\frac{1}{\bar{\unicode[STIX]{x1D703}}}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}x}\right)\bar{C}_{n}+\left(\unicode[STIX]{x1D6FF}_{B}+\frac{\unicode[STIX]{x1D6FF}_{T}}{\bar{\unicode[STIX]{x1D703}}}\right)\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}x}\frac{\unicode[STIX]{x2202}\bar{C}_{n}}{\unicode[STIX]{x2202}x}.\end{eqnarray}$$
The rate of change of NP concentration is the result of three terms. The first one corresponds to a diffusion with a temperature-dependent diffusivity, while the second term, whose strength is directly proportional to the thermophoretic diffusivity, acts as either a source or a sink term depending on the sign of
$\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x((1/\bar{\unicode[STIX]{x1D703}})(\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x))$
. It is clear that since
$\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x((1/\bar{\unicode[STIX]{x1D703}})(\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x))=\unicode[STIX]{x2202}^{2}\ln \bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x^{2}$
is not constant, the strength of the source/sink term varies in the channel. The last term, can be considered as a convective term with an effective velocity that depends on the temperature gradient and the Brownian and thermophoretic diffusivities. It can also be regarded as the difference of the NP frontal speed compared to that of the solutal/thermal front in the case
$\unicode[STIX]{x1D706}=1$
and negligible thermophoretic and Brownian effects in the energy equation (see § 5.4). The expression of this difference in the fronts’ speeds at
$x=0$
is
$$\begin{eqnarray}u_{re}=\frac{1-r}{2\sqrt{\unicode[STIX]{x03C0}tLe}}\left(\unicode[STIX]{x1D6FF}_{B}+\frac{2\unicode[STIX]{x1D6FF}_{T}}{r+1}\right).\end{eqnarray}$$
Clearly
$u_{re}=0$
in the isothermal case
$(r=1,\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x=0)$
, implying that all fronts travel at the same speed. The expression of
$u_{re}$
also indicates that the higher
$Le$
, the smaller this speed difference and the ensuing effects on the instability.
This front speed is positive in the HDC case
$(r<1,\unicode[STIX]{x2202}\unicode[STIX]{x1D703}/\unicode[STIX]{x2202}x<0)$
, implying that the NP front travels ahead of the solutal/thermal front resulting in the downstream transport of nano-particles. This difference in the fronts’ speeds increases with increasing thermophoretic or Brownian diffusivity. As
$\unicode[STIX]{x1D6FF}_{T}$
increases, the downstream transport of NPs towards the cold fluid increases, bringing the NPs into a flow region where their convective velocity is essentially zero
$(\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x\rightarrow 0)$
. Therefore, the NPs in the region with larger absolute temperature gradients are flowing downstream faster
$(\boldsymbol{V}_{T}\propto -\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x)$
while those in the region of cold fluid are moving slowly, leading to the local accumulation of NPs. In fine, the NP front is convected further downstream to a flow region with a high local value of the source term
$\unicode[STIX]{x1D6FF}_{T}(\unicode[STIX]{x2202}^{2}\ln \bar{\unicode[STIX]{x1D703}})/\unicode[STIX]{x2202}x^{2}>0$
resulting in a local increase in the NP concentration, as already seen in figure 18 in the case
$\unicode[STIX]{x1D6FF}_{T}=20$
.
Base state NP concentration and scaled variations of
$\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x$
and
$\unicode[STIX]{x2202}^{2}\ln \bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x^{2}$
along the channel at
$t=1500$
in the CDH system.
To complete the discussion the NP accumulation mechanism in the CDH system is examined. Figure 20 depicts the base state NP concentration variation along the channel for different
$\unicode[STIX]{x1D6FF}_{T}$
with the representative scaled variation of
$\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x$
and
$\unicode[STIX]{x2202}^{2}\ln \bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x^{2}$
. Both
$\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x$
and
$\unicode[STIX]{x2202}^{2}\ln \bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x^{2}$
decay to zero at both ends of the domain. However,
$\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x$
has a positive local maximum at the centre of the channel, while the variation of
$\unicode[STIX]{x2202}^{2}\ln \bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x^{2}$
exhibits a positive local maximum at
$x<0$
and a negative local minimum at
$x>0$
. The figure shows that as
$\unicode[STIX]{x1D6FF}_{T}$
increases, NPs tend to accumulate at the front and the location of both NP front and
$(\bar{C_{n}})_{max}$
move upstream. As in the HDC case, these trends can be understood by examining the NP base state equation, equation (5.2). Since
$\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x>0$
in this case, the convective term
$(\unicode[STIX]{x1D6FF}_{B}+\unicode[STIX]{x1D6FF}_{T}/\bar{\unicode[STIX]{x1D703}})\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x$
is always positive (or
$u_{rel}<0$
), implying that the NP front convects more upstream as
$\unicode[STIX]{x1D6FF}_{T}$
increases and leads to a stronger presence of NPs in the upstream region
$x<0$
, and in turn a weaker presence in the downstream section;
$x>0$
. However, in
$x<0$
, NPs have different thermophoretic velocities, with stronger backward velocity due to the higher positive
$\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x$
(
$\boldsymbol{V}_{T}\propto -\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x$
) close to
$x=0$
, and an almost zero thermophoretic velocity far away from the front. As a result, the NPs tend to accumulate at the front and this accumulation intensifies more as
$\unicode[STIX]{x1D6FF}_{T}$
increases.
The accumulation is a result of thermophoresis and is strengthened with
$\unicode[STIX]{x1D6FF}_{T}$
, since as the thermophoretic diffusivity increases, the source term
$\unicode[STIX]{x1D6FF}_{T}(\unicode[STIX]{x2202}^{2}\ln \bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x^{2})$
gets larger at
$x<0$
and also more NPs are convected upstream
$((\unicode[STIX]{x1D6FF}_{B}+\unicode[STIX]{x1D6FF}_{T}/\bar{\unicode[STIX]{x1D703}})\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D703}}/\unicode[STIX]{x2202}x)$
. Although an increase of the Brownian diffusivity results in larger NP frontal speed
$u_{re}$
and consequently a stronger upstream transport of NPs, it is also accompanied by stronger Brownian diffusion which explains the observed trend towards smaller NP accumulation as the Brownian effects gets stronger (see figure 5 for CDH case).
The previous results indicate that Brownian and thermophoretic effects have similar destabilizing effects in the HDC system, where larger values of
$\unicode[STIX]{x1D6FF}_{B}$
or
$\unicode[STIX]{x1D6FF}_{T}$
lead to a more unstable displacement. The effects of these two mechanisms are however different in the CDH case, where a stronger Brownian diffusion tends to increase the instability while larger thermophoretic effects actually have a stabilizing effect, at least in some ranges. In particular, a large accumulation of NPs, leading to local stable regions upstream of the front, is observed in the CDH case for smaller
$\unicode[STIX]{x1D6FF}_{B}$
and larger
$\unicode[STIX]{x1D6FF}_{T}$
, which hinders the growth of backward-developing fingers. These results point to the close interplay of these two effects that can either act synergistically or have counter-effects. In the HDC case, Brownian diffusivity and thermophoresis act together towards the migration of the NPs in the downstream direction, while in the CDH case where the temperature gradient is opposite to that of the HDC case, Brownian diffusion still tends to transport NPs downstream of the flow towards the hot fluid while now thermophoresis acts in the opposite direction resisting the downstream migration of the NPs. Accordingly, the counter-effects of thermophoresis in the CDH case are more noticeable for smaller
$\unicode[STIX]{x1D6FF}_{B}$
or larger
$\unicode[STIX]{x1D6FF}_{T}$
. In the HDC case, on the other hand, the two mechanisms have synergistic effects, and larger values of
$\unicode[STIX]{x1D6FF}_{B}$
and/or
$\unicode[STIX]{x1D6FF}_{T}$
lead to the migration of NPs in the downstream direction towards the cold less viscous fluid, resulting in a more unstable displacement.
5.3 Thermophoretic effects in stable/intrinsically stable systems
The previous results and physical interpretations lead to the natural question of whether thermophoresis effects would destabilize a stable or intrinsically stable system. To answer this question, the instability of both HDC and CDH systems with large nanofluid viscosity ratio,
$R_{n}=5$
, is investigated in the presence and the absence of thermophoretic effects. Note that such a flow is unstable (
$R_{b}=3$
,
$R_{\unicode[STIX]{x1D703}}=1$
(HDC),
$R_{\unicode[STIX]{x1D703}}=-1$
(CDH)) in the absence of NPs.
Figure 21 shows the concentration contours in the HDC system in the presence and absence of thermophoretic effects. It is observed that the system is stable in the absence of thermophoresis
$(\unicode[STIX]{x1D6FF}_{T}=0)$
, while it is unstable otherwise
$(\unicode[STIX]{x1D6FF}_{T}=10)$
. The instability develops as a result of the advancement of the NP front downstream, which then leads to the non-monotonicity in the viscosity distribution. Therefore, one can extend the enhancement of instability in an already unstable flow to a destabilization one of a stable flow, when one is dealing with the HDC case. The analysis did not show this effect in the CDH system, such that the system was stable both in the presence and the absence of thermophoresis.
Contours of
$C_{b}$
in the HDC system with
$R_{n}=5$
.
The problem is now extended to an intrinsically stable system. The choice of
$R_{b}=-1$
(more viscous fluid displacing a less viscous one) represents the intrinsically stable HDC system
$(R_{\unicode[STIX]{x1D703}}=1)$
in the absence of NPs
$(R_{n}=0,\,\unicode[STIX]{x1D716}=0)$
. Note that the default values of the other parameters are still used. Figure 22 reveals that although the addition of NPs further increases the viscosity of the displacing fluid and may lead one to expect the flow stability to be maintained, it is found that the dynamics are changed and fingers start to develop in the flow as a result of thermophoretic effects (
$\unicode[STIX]{x1D6FF}_{T}=20$
). In other words, such a system became unstable only as a result of the addition of NPs and their migration due to the already present temperature gradient. This shows that thermophoresis can even make an intrinsically stable HDC system, unstable. The CDH system with
$R_{b}=-1$
did not show this effect.
Contours of
$C_{b}$
in the HDC system for
$R_{b}=-1$
in (a) the absence, (b) the presence of NPs. Note that
$R_{n}$
and
$\unicode[STIX]{x1D716}$
are set to zero to generate the case without NPs.
5.4 Thermophoresis and Brownian diffusion in the energy equation
The previous results and interpretation were based mainly on changes in the distribution of NPs as a result of Brownian and thermophoretic effects. However, these effects are also present in the energy equation, and it is legitimate to ask if the dynamics can be explained by an analysis of the contributions of different terms in the energy equation, (3.15). The last two terms on the right-hand side of this equation represent the Brownian and thermophoresis contributions. However, their effects on the instability depend directly on the value of
$\unicode[STIX]{x1D716}=C_{n_{0}}\unicode[STIX]{x1D719}(\unicode[STIX]{x1D70C}c_{p})_{p}/(\unicode[STIX]{x1D70C}c_{p})_{m}$
which has not been addressed yet.
Note that the volumetric heat capacity
$(\unicode[STIX]{x1D70C}c_{p})$
of both solids and liquids are of the order of
$10^{6}$
. Furthermore, the maximum order of magnitude of the NP mass/volume fraction and the porosity are typically
$10^{-2}$
and
$10^{-1}$
, respectively. This leads to the conclusion that
$O(\unicode[STIX]{x1D716})=O(C_{n_{0}}\unicode[STIX]{x1D719})$
, or
$\unicode[STIX]{x1D716}<0.01$
. Adopting
$\unicode[STIX]{x1D716}=0.01$
, the contributions of the Brownian motion and thermophoresis are analysed while other parameters are fixed. Examinations of contours and maximum and average values of the diffusive term
$Le\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}$
and of the combined Brownian–thermophoretic term
$\unicode[STIX]{x1D716}(\unicode[STIX]{x1D6FF}_{B}\unicode[STIX]{x1D703}\unicode[STIX]{x1D735}C_{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}+\unicode[STIX]{x1D6FF}_{T}C_{n}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}/\unicode[STIX]{x1D703}))$
revealed that the latter term is order of magnitudes smaller than the former, for both HDC and CDH cases. The analysis revealed almost the same conclusion even with
$\unicode[STIX]{x1D716}=0.1$
. This result was found systematically over a wide range of the other flow parameters. It clearly points to the negligible contribution of the Brownian motion and thermophoresis in the energy equation and to the fact that such mechanisms can actually be ignored in the modelling and analysis of the energy effects, at least for the system at hand. It also implies that the main effects of Brownian diffusion and thermophoresis arise from their contribution to the NP transport while their contribution through the energy equation can be safely disregarded:
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t}+(\unicode[STIX]{x1D706}-1)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}x}+\unicode[STIX]{x1D706}(\boldsymbol{V}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703})=Le\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D703}.\end{eqnarray}$$
6 Conclusion
Flow instability in homogeneous porous media is investigated in the presence of thermal effects and nanoparticles dispersed in the displacing fluid. The analysis is conducted for two representative systems: HDC, where a hot fluid displaces a cold fluid, and CDH, where a cold fluid displaces a hot fluid.
In the HDC system the NP concentration and temperature gradients are in the same direction. As a result, Brownian diffusion and thermophoresis act synergically to move NPs downstream towards the cold fluid in the low resistance paths created by the forward-growing fingers. It is found that increasing
$\unicode[STIX]{x1D6FF}_{T}$
makes an initially unstable NP-laden system more unstable, which is similar to the effects of
$\unicode[STIX]{x1D6FF}_{B}$
. The destabilizing effects of thermophoresis in the HDC flows are extended to initially stable and intrinsically stable systems. On the other hand, the gradient of NP concentration and temperature in the CDH system are in opposite directions. As a result, Brownian diffusion still tends to transport the NPs downstream towards the hot fluid while thermophoresis resists such migration. These counter-effects lead to a local upstream accumulation of NPs and their depletion downstream, at higher
$\unicode[STIX]{x1D6FF}_{T}$
but lower
$\unicode[STIX]{x1D6FF}_{B}$
. Because of the NP accumulation, local stable positions are generated which in turn affect the growth of the instabilities, especially those of backward-growing fingers. Accordingly, initially unstable NP-laden CDH systems are less unstable in the presence of thermophoresis compared to the case without thermophoresis, although Brownian diffusivity still has destabilizing effects.
Two counter-intuitive phenomena are observed and discussed in both HDC and CDH systems. First, although Brownian diffusion and thermophoresis transport NPs downstream synergically in the HDC system, local NP accumulation spots have been observed in the flow after a critical value of
$\unicode[STIX]{x1D6FF}_{T}$
. This was attributed to the advancement of the NP front downstream into the cold fluid region for large enough values of
$\unicode[STIX]{x1D6FF}_{T}$
. As a result, the viscosity profile becomes non-monotonic at large enough values of
$\unicode[STIX]{x1D6FF}_{T}$
and
$R_{n}$
which can cause the destabilization of initially stable HDC systems. Second, two regimes of instabilities are observed in the CDH case. For a range of small enough values of thermophoretic diffusivity, the flow instability was found to be attenuated with increasing
$\unicode[STIX]{x1D6FF}_{T}$
. These trends are however reversed beyond a certain critical value of
$\unicode[STIX]{x1D6FF}_{T}$
. Interestingly, these non-monotonic effects of
$\unicode[STIX]{x1D6FF}_{T}$
on the instability in the CDH case are systematically accompanied by a suppression of the growth of backward-developing fingers as
$\unicode[STIX]{x1D6FF}_{T}$
increases. This change in the trends was explained in terms of the development of regions resisting backward flow that can lead to an enhancement of the development of forward-growing fingers. One must note that the general trends towards an attenuation of the instability are still observed in the presence of thermophoresis compared to the case without thermophoresis.
It is also mentioned that the NP accumulation due to the mechanisms discussed is intensified for low
$\unicode[STIX]{x1D6FF}_{B}$
while it is mitigated under smooth temperature distributions, observed in particular at large Lewis numbers. Finally, it was found that the main effects of Brownian diffusion and thermophoresis arise from their contribution to the NP transport while their contribution to the energy distribution can be safely disregarded.
Acknowledgements
This work is supported by the University of Calgary’s Eyes High Scholarship programme and the Natural Sciences and Engineering Research Council of Canada (NSERC). The authors also acknowledge the use of computing resources of the West-Grid clusters.
