## Abstract

In this paper, the stability of unidirectional Couette flows of a spatially inhomogeneous complex fluid is investigated via a linear stability analysis. The employment of a classical normal-mode analysis results in a fourth-order polynomial eigenvalue problem that is solved numerically via a Chebyshev polynomial method. The results of this study suggest that the flow of interest is unconditionally unstable but nevertheless the growth rates remain small so that the observability of the instability requires a very large time window. The observed instability manifests itself through the formation of concentration waves that, due to their location at the bottom of the channel, closely resemble a series of superposed, high-concentration piles.

## 1. Introduction

Conglomerates of non-Brownian particles, dense blood suspensions and, more generally, partially ordered mixtures constitute paradigms of spatially inhomogeneous complex fluids [1,2]. Such fluids are associated with a free energy that depends on both the volume fraction and the volume fraction's gradient, and their behaviour is substantially intricate. For example, while they can respond to external forces in a fluid-like manner, they can support strain-independent shear stresses at zero shear rates. Moreover, upon attaining their maximum concentration, they experience jamming, i.e. a form of phase change during which only rigid-body motion is possible.

The explicit dependence of the free energy on the volume fraction and the volume fraction gradient has important consequences in the modelling of the materials of interest. More specifically, as observed by Goodman & Cowin [3], it dictates the inclusion of these quantities as independent thermodynamic variables in the state vector of the inhomogeneous complex fluid; see also the relevant discussions of Kirchner [4], Papalexandris [5] and Ván *et al.* [6]. In turn, this allows for the derivation of a rate equation for the volume fraction which, when supplemented with the balance laws of mass, momentum and energy, constitutes a closed system. Furthermore, the constitutive expression for the stress tensor of the inhomogeneous complex fluid comprises an additional, non-dissipative component, the so-called ‘configuration stress tensor’. This tensor involves volume fraction gradients that are associated in a non-Newtonian way and models ‘configuration stresses’. These are stresses induced by the rearrangement in the distribution of the interfacial density. At equilibrium, its off-diagonal components are identified with the strain-independent shear stresses mentioned above.

In general, flows of inhomogeneous complex fluids are manifestly dilatational. Indeed, even if the density is assumed to be constant, as in most cases of interest, dilatancy can be engendered by changes in the volume of the interstices and can affect the distribution of the volume fraction; see, for instance, Goodman & Cowin [3] and Passman *et al.* [7]. However, as the concentration increases, the effects of dilatancy become smaller. Then, as stressed by Málek & Rajagopal [2], sufficiently dense flows of inhomogeneous complex fluids can be approximated as isochoric ones.

Even if the flow is isochoric, the materials of interest can still support concentration inhomogeneities. Consequently, the configuration stresses can in principle be of the same order of magnitude as those of inertial, viscous or gravitational forces. In other words, independent of the strength of dilatancy, the configuration stresses can significantly modify the properties of the flow by either triggering or suppressing hydrodynamic instabilities. Nonetheless, there is a dearth of studies that aim to quantify the interplay between the configuration stresses, inertia, viscosity and gravity. As a consequence, the stability properties of incompressible flows of inhomogeneous complex fluids remain rudimentarily understood.

This is precisely what we are concerned with in this paper. More specifically, we investigate the stability of plane, unidirectional Couette flows of an inhomogeneous complex fluid by means of a linear stability analysis. Our study is based on a particular continuum model for the flows of interest, previously derived by Papalexandris [5], and follows the steps of a classical normal-mode analysis. Accordingly, our focus pertains to the asymptotic stability of Couette flows subject to infinitesimal perturbations.

The objective and novelty of this study is twofold. First, to systematically investigate the stability properties of incompressible Couette flows of inhomogeneous complex fluids. For this reason, extensive parametric studies, with respect to the dimensionless groups that appear in the governing equations, are also performed. The secondary objective of this study is to provide further insight into the role of the configuration stress tensor and, thus, the volume fraction's gradient, in transient flowing conditions. Thus far, relevant studies, such as the ones of Passman *et al.* [7], Wang & Hutter [8] and Massoudi & Boyle [9], have focused exclusively on steady-flow conditions.

This paper is organized as follows. In §2, we describe the governing equations of an inhomogeneous complex fluid and derive the base-flow profile. Section 3 is devoted to the linear stability analysis and numerical methodology that is employed for the treatment of the resulting eigenvalue problem. In §4, we present and discuss the numerical results and determine the stability diagrams for the problem in hand. Finally, §5 concludes.

## 2. Governing equations

Governing equations for incompressible flows of inhomogeneous complex fluids can be obtained by taking the single-phase limit of the continuum model for incompressible fluid-saturated granular flows of Varsakelis & Papalexandris [10]. This system of equations constitutes the zero-Mach-number limit of the mixture theory for coexisting and interacting continua of Papalexandris [5], whose free energy depends on both their bulk density and its spatial gradient.

The systematic reduction of two-phase flow models to their single-phase limit is a delicate task due to the singularity of this limit. Herein, we follow the steps of Bdzil *et al.* [11]. Then, upon doing so, we arrive at the following system of (dimensional) equations:

*Mass and momentum balance equations*

*Compaction equation*
*ϕ* and

The term *configuration stress tensor* discussed in the Introduction. Accordingly, the term *configuration stress coefficient*. The configuration stress tensor bears strong resemblances to the so-called Korteweg capillarity tensor. This is not surprising because the continuum theory that equations (2.1)–(2.3) have been built upon is to a large extent preceded by the theory of Korteweg [12].

For the rheology of the complex fluid *ϕ*_{M} represents the maximum concentration upon which the complex fluid experiences jamming. As a matter of fact, equation (2.4) constitutes an incarnation of the Krieger–Dougherty rheology which is widely exercised in the study of suspensions; see, for instance, the review article of Stickel & Powel [14]. For more information about particle jamming and its properties, the interested reader is referred to Liu & Nagel [15].

On the other hand, and to the best of our knowledge, systematic experimental measurements for *et al.* [7], who argued that, as

Finally, it is worth noting that the well-posedness of equations (2.1)–(2.3), endowed with a quite general volume-fraction-dependent rheology and configuration stress vector, has been established in various cases. For the equilibrium limit, we refer to Varsakelis & Papalexandris [16], whereas for the unsteady case we refer to Nakano & Tani [17].

### (a) Non-dimensionalization, base-flow profiles and boundary conditions

In the configuration of a plane, unidirectional Couette flow, a complex fluid is placed between two infinite parallel planes in the

For the purposes of a stability analysis, it is convenient to work with the non-dimensionalized version of the governing equations so that the effect of the various non-dimensional groups is apparent. In order to do so, we first determine the (dimensional) reference values that will be denoted with the subscript ‘r’. The reference length, *Re* and *Fr*, corresponding to the Reynolds and Froude numbers of the flow, are defined as usual,
*ϕ*. Then,

By performing the above non-dimensionalizations and by introducing the expressions for the rheology

We now proceed to the determination of the base-flow profiles. Following standard practices, we postulate that the base flow is time-independent and quasi-parallel, i.e. that its variation is located entirely in the direction normal to the streamwise direction. On this premise, we consider a base flow of the following form:

The introduction of (2.13) to (2.10)–(2.12) yields, upon rearrangement,

The postulate for a stationary, quasi-parallel flow asserts that the compaction equation (2.3) is identically satisfied. Consequently, the system of ODEs given by (2.14)–(2.15) is underdetermined. Málek & Rajagopal [2] employed arguments of kinematics to establish that this underdeterminacy is an inherent property of the materials of interest. More specifically, it results from the correlation between the initial configuration of the complex fluid and incompressibility. In other words, the solution to the above ODEs requires the *a priori* knowledge of either the volume fraction distribution, the ambient pressure or the velocity.

In the absence of any further guidance, we choose to prescribe the concentration. In particular, we consider a linear, stable stratification
*ϕ*_{1}<*ϕ*_{2}<1 are constants, so that the concentration decreases with height. We note that if the base-flow concentration is homogeneous, then in the context of a linear stability analysis the configuration stress tensor is rendered redundant as a higher-order term. Moreover, such a choice, combined with the quasi-parallel flow assumption, completely decouples the compaction equation from the momentum equation, thus rendering the concentration a passive scalar. It is interesting to note that experimental measurements on Couette dry granular flows, in the presence of gravity, have documented that, in steady state, the concentration decreases monotonically from the top of the channel up until a borderline which lies away from the bottom of the channel and whose location depends on the material under study and the configuration of the flow. Below this borderline the concentration remains constant [18]. Consequently, opting for a linear profile, which constitutes a first approximation to the curved profile observed in GDR-MIDI [18], amounts to zooming in on the area of the flow where non-negligible concentration gradients are developed.

Next, by introducing (2.16) to (2.14) and after some algebra, we arrive at the following profile for the streamwise component of the base-flow velocity:
*c*_{1} and *c*_{2} being integration constants whose values are fixed once the boundary conditions have been prescribed. Even though the above expression appears to be quite involved, the base-flow velocity is actually a simple convex function, reminiscent of a half parabola. Nevertheless, due to the nonlinear correlation between *μ*_{s} and *ϕ*_{s}, it is not possible to express the base-flow velocity in a more succinct form.

Subsequently, upon insertion of (2.16) and (2.17) into (2.15), one can directly acquire the expression for the base pressure *p*^{0}(*y*), which reads

In order to formulate a well-posed boundary-value problem, the governing equations need to be endowed with boundary conditions that are consistent with the geometry of the Couette flow. Owing to the quasi-parallel assumptions, boundary conditions are required only for the upper and the lower planes. Herein, we assume that the complex fluid velocity **u** obeys the wall-boundary condition at both the upper and the lower planes, namely **u**(0)=(0,0) and **u**(1)=(1,0). On the other hand, Dirichlet boundary conditions are prescribed for the volume fraction at both planes and, in particular, we set *ϕ*(0)=*ϕ*_{2} and *ϕ*(1)=*ϕ*_{1}. Moreover, we postulate that **u**^{0} and *ϕ*^{0} satisfy the same boundary conditions as **u** and *ϕ*, respectively, i.e. we set
*c*_{1} and *c*_{2} in the base-flow velocity **u**^{0} are uniquely determined.

It is worth noting that one may not assign arbitrary boundary conditions to the pressure field *p*. In fact, as shown in Varsakelis & Papalexandris [19], appropriate boundary conditions for *p* are induced by the Helmholtz decomposition and Ladyzhenskaya's decomposition theorem. Consequently, imposing additional boundary conditions leads to an overdetermined problem. The situation is completely analogous to the Navier–Stokes equations; see, for example, Ladynzhenskaya [20] for a thorough exposition.

## 3. Linear stability analysis

In this section, we study the linear stability of the base flow profile, given by relations (2.17) and (2.18), subject to infinitesimal disturbances. In order to do so, we seek solutions expressed as the superposition of the base flow and a disturbance **u**^{1}=**u**^{1}(*x*,*y*,*t*), *p*^{1}=*p*^{1}(*x*,*y*,*t*) and *ϕ*^{1}=*ϕ*^{1}(*x*,*y*,*t*), i.e.
**L** is an algebraic operator that depends also on the base flow. The exact expression of **L** is too lengthy to reproduce here, and offers little insight, but is available to the interested reader upon request.

The system of the linearized equations given by (3.2) is valid for small but otherwise arbitrary disturbances. Herein, we opt for a classical normal-mode analysis, as detailed in Drazin & Reid [21], and assume that disturbances are travelling waves,
*α* is the prescribed real wavenumber and *c*=*c*_{r}+i*c*_{I} is the complex frequency. Instability is associated with a wave speed *c* such that *c*_{I}>0. Also, as we are concerned with two-dimensional incompressible flows, it is advantageous to employ the stream-function representation of the disturbance of the velocities

where we have employed the notation convention *A*_{ij}, *i*,*j*=0,…,4 are quite lengthy and complex. Nevertheless, for the sake of completeness, they have been included in appendix A. We stress that their computation and validation has been performed via the symbolic package Maxima [22].

The boundary conditions for

Equation (3.5) constitutes a non-self-adjoint eigenvalue problem with *c* and *ϕ*_{s} and the presence of the configuration stress tensor, the resulting problem belongs to the class of fourth-order polynomial eigenvalue problems. Before proceeding to the solution of the eigenvalue problem, it is pertinent to exploit the properties of equation (3.5), especially as it has not appeared in the literature before.

First, equation (3.5) can be seen as the equivalent to the Orr–Sommerfeld equation for the flows of interest. In fact, it can be shown that if *ϕ*^{0}≡*cte*, everywhere in the domain, and the rheology law is taken equal to *μ*_{s}=1, then the standard Orr–Sommerfeld equation is recovered. However, this result is neither coincidental nor surprising. Indeed, for these parameter values, the governing equations (2.1)–(2.3) reduce to the Navier–Stokes equations.

Further, we claim that the spectrum of the operator *y* is a complex variable. Then, the base flow velocity *u*^{0}(*y*) and volume fraction *ϕ*^{0}(*y*) are analytic functions of the variable. Now, equation (3.5) is regular in *y*, *c*, *a*, *Re*, *Fr* and *c*, form a totally disconnected set.

### (a) Numerical method

Owing to the complexity of equation (3.5), the corresponding boundary value problem (3.5) and (3.6) is integrated numerically as follows. First, the fourth-order polynomial eigenvalue problem (3.5) is rewritten as an equivalent generalized eigenvalue problem. Then, every function is evaluated numerically via a Chebyshev-polynomial interpolation while operators are approximated via a pseudo-spectral collocation method. The discretized eigenvalue problem is then solved by QR-decomposition on an adaptive grid so that oscillatory, grid-dependent eigenvalues are ruled out.

The employed numerical methodology is standard and has already been implemented in ready-to-use Matlab codes [24]. For our numerical simulations, we have used the Chebfun package [25], the accuracy and robustness of which have been assessed in numerous tests, including the efficient computation of the spectrum of the Orr–Sommerfeld operator.

Validation of the code has been performed as follows. We have numerically solved the Orr–Sommerfeld equation for various values of the Reynolds number. Subsequently, we have approximated numerically the reduced version of the eigenvalue problem (3.5), which coincides with the Orr–Sommerfeld one, and have verified the agreement between the computed spectra.

## 4. Numerical results

The numerical solution of the eigenvalue problem in hand requires that the various physical parameters and coefficients, appearing in equations (2.1)–(2.3), are assigned suitable values. Moreover, these values should ideally be obtained through systematic experimental measurements.

Herein, we choose to investigate the linear stability of a generic inhomogeneous complex fluid whose physical parameters resemble those of natural, angular beach sand because reliable values for most parameters have been provided in Passman *et al.* [7] and Wang & Hutter [8]. Then, the density of the complex fluid is set equal to *ϕ*_{M}=0.74. The value of the coefficient

As regards *ϕ*_{1} and *ϕ*_{2} in equations (2.16)–(2.18), we first recall that the isochoric approximation is valid only for very high concentrations. On the other hand, if the values of the concentration are in the vicinity of *ϕ*_{M}, the rheology of the complex fluid practically blows up. In turn, this significantly decreases the condition number of the discretized operators and thus deteriorates the performance of the numerical method. In view of these constraints, we set *ϕ*_{2}=0.7 and *ϕ*_{1}=0.6 so that an approximately 15% variation in the concentration is imposed.

The remaining values that need to be determined are the velocity of the upper plane *Re*=0.0031, *Fr*=1.19×10^{−4} and *α*, we initially set it equal to unity.

Figure 1*a*,*b* depicts the dispersion relation for the reference case, i.e. it portrays the growth rate *σ*=*c*_{I} as a function of the wavenumber *α*. It can be readily inferred that the linear stability analysis predicts the existence of unstable modes for all values of *α*∈(0,1]; even the value that seems to be a neutral stability point is actually positive. Still, the maximum growth rate, being of the order of 10^{−5}, is relatively small. The statistical trend of the dispersion curve suggests that the growth rate increases with the wavenumber, so that short-wave instabilities lead to higher growth rates. However, the observed correlation is not monotonic but rather follows an oscillatory pattern whose amplitude is negatively correlated with *α*. The subtle dependence of the growth rate *σ* on the wavenumber *α* is attributed to the fact that the eigenvalue problem (3.5) contains powers of *α* up to seventh order. Then, since *α*<1, as *α* decreases so does the power of the monomial of *α* that dominates the overall expression.

Additional information about the properties of the predicted instability can be obtained by the eigenfunction of the most ‘dangerous’ eigenvalue that constitutes the amplitude of the stream-function's perturbation, *α*. The profiles of *α* and thus the steepness of the profiles also becomes milder as *α* increases.

In order to gain further insight into the manifestation of the instability, in figures 3 and 4 we have plotted contour plots of the concentration and the vorticity field, respectively, at *t*=0.1. Both figures correspond to a perturbation with wavenumber *α*=1.

Upon the commencement of the instability, the motion of complex fluid quickly becomes two-dimensional and, similarly, non-trivial concentration gradients form. As a consequence, the profile of the concentration deviates from the initial linearly stratified one and a series of elongated, plane-parallel concentration waves is formed. The (lower concentration) trough of the waves is located in the middle of the channel and precedes and supersedes two (higher concentration) crests. Owing to the positioning of the waves, the onset of which is from the bottom of the channel where the concentration is higher, they bear a resemblance to half-piles, which, in view of the periodicity in the streamwise direction, are actually full ones if a periodic extension is performed. In this respect, the height of these piles is approximately 70% of the Couette gap.

According to figure 4, as the flow evolves two pairs of counter-rotating vortices that are antisymmetric both in the streamwise and in the normal direction are developed. These vortical structures are stretched in the streamwise direction with a ratio of diameters of approximately 7:1. Moreover, it is interesting to note that the area of minimum vorticity is located in the cavity, where the concentration is lower. Consequently, the motion of the complex fluid is characterized by a form of a preferential concentration.

### (a) Parametric study

In order to identify the mechanism of the predicted instability and additionally quantify the interplay between viscous forces, gravitational forces and the configuration stresses, a series of parametric studies have been carried with respect to the three dimensionless groups *Re*, *Fr* and

First, we recall that for the reference case we have taken *Re* varies from approximately 3×10^{−5} to 0.1, whereas the Froude number *Fr* varies from 10^{−5} to 0.01. On the other hand, as regards ^{−5} to 10^{−2}, so that ^{−5} to 0.1.

Figures 5–7 depict the results of the parametric study with respect to the Reynolds number *Re*, the Froude number *Fr* and the parameter ^{−5}≤*Re*≤0.1, 10^{−5}≤*Fr*≤0.01 and ^{−5}, suggesting that the observability of the predicted instability requires a quite large time window.

Besides asserting the unconditional instability of the flow under study, figures 5–7 also illustrate the correlation between the growth rate *σ* and the dimensionless groups. More specifically, we observe that *σ* is predicted to correlate negatively with *Re* and *Fr* and positively with

The flow under study belongs to the class of concentration and viscosity-stratified flows, i.e. both the partial density *ρϕ* and the viscosity *μ* support variations in the normal direction. As shown, among others, in Yih [26] and Hickox [27], viscosity-stratified Couette flows admit unstable modes. Moreover, as asserted by Ern *et al.* [28], such modes are present even if viscosity varies in an abrupt, albeit continuous, manner. In other words, for such flows, viscosity has been well documented to constitute a destabilizing mechanism.

According to our findings, for the flows of interest, viscous forces are also found to constitute a mechanism of instability; hence, the negative correlation between the magnitude of the Reynolds number and the growth. This result appears ostensibly to not be surprising and in line with the aforementioned ones; however, some caution has to be exercised. All the afore-cited studies, as well as those antedating and succeeding them, are typically based on the utilization of regularized, piecewise constant expressions for the density or the viscosity. On the other hand, as inferred from (2.4), the correlation employed for the rheology of the complex fluid is markedly more complex, being a nonlinear function in the volume fraction and having a singularity at the maximum concentration. In this respect, although one might conjecture that the results of Yih [26], Hickox [27] and Ern *et al.* [28] can be extrapolated to more involved expressions a formal confirmation or refutation is required. In this respect, the analysis in the present study provides this formal verification.

As mentioned in the Introduction, a primary impetus of this study has been the quantification of the role of the configuration stresses and more generally the effects of the free energy of a complex fluid being a function of the volume fraction and its gradient. As figure 7 illustrates, there is a positive correlation between *σ* and

In this study, the Froude number is quite small. Therefore, given the absence of concentration gradients in a direction opposite to gravity, the negative positive correlation between *Fr* and *σ* portrayed in figure 6 is expected and in line with previous findings in stably-stratified shear flows [30]. In other words, gravitational forces provide a stabilizing mechanism in the flow. However, for the parameter values considered herein, gravitational forces are not large enough to offset the destabilization induced by viscosity stratification and the configuration stresses.

## 5. Conclusion

In this paper, the asymptotic stability of unidirectional, plane Couette flows of an inhomogeneous complex fluid, as predicted by the continuum model of Varsakelis & Papalexandris [10], has been investigated via a linear stability analysis. The base flow comprises a stable, linearly stratified concentration, a semi-parabolic velocity and a pressure that is different from the standard hydrostatic one in order to accommodate the effects of the configuration stresses. By utilizing a classical normal-mode analysis, the stability modes have been identified with the solutions of a polynomial eigenvalue problem that have been approximated numerically.

According to our numerical predictions, the flows of interest are unconditionally unstable for all values of the Reynolds number, the Froude number and the novel dimensionless number ^{−5}. Consequently, the time window required for the observability of the predicted instability is quite large. The instability manifests itself through the formation of a series of plane-parallel concentration waves that, due to their location, bear strong resemblances to high-concentration piles.

Both the viscosity and the configuration stresses are found to be the governing mechanisms that drive the predicted instability. By contrast, due to the initial concentration being stably stratified, combined with the fact that perturbations do not induce volume fraction gradients pointing to the opposite of the gravity vector, gravity acts as a partial stabilizer. However, even though its magnitude is inadequate to completely stabilize the flow it is adequate to maintain growth rates of the order of 10^{−5}.

The stability analysis carried in this study constitutes the first step towards understanding the role of the volume fraction gradients in the stability of transient flows of spatially inhomogeneous complex fluids. However, given the complexity of the governing equations, nonlinear effects are also expected to play a prominent role in the determination of the evolution of the disturbances. In this respect, an important direction of research concerns the investigation of the nonlinear stability of the flows of interest via direct numerical simulations. Besides quantifying the role of the nonlinearities, such a numerical study will additionally serve as a tool for assessing the domain of validity of the linear stability analysis reported in this study. We intend to pursue this research in an upcoming study.

## Authors' contributions

C.V. and M.V.P. conceived and carried out the stability analysis, interpreted the numerical predictions and wrote the manuscript. C.V. implemented the code and conducted the numerical experiments in close collaboration with M.V.P. Both authors gave final approval for the publication of this manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

Financial support for C.V. has been provided by the National Research Fund of Belgium (FNRS) under a GRANMIX Projet de Recherche Grant.

## Acknowledgements

The authors thank the anonymous referees whose comments and suggestions have improved the quality of this paper.

## Appendix A

In this appendix, we provide the expressions for the coefficients of the eigenvalue problem (3.5)

- Received July 31, 2015.
- Accepted September 23, 2015.

- © 2015 The Author(s)