Suspensions of deformable particles in Poiseuille flows at finite inertia

Particles suspended in a carrier fluid can be encountered in several domains, such as biology and geophysics, but also biomechanical, pharmaceutical and cement industries. Pyroclastic flows from volcanoes, fluidized beds, the blood flow in veins and arteries, sedimentation in sea beds and slurry flows are other examples of flows where the knowledge of the particle dynamics is relevant. Currently, it is still hard to estimate the pressure drop needed to drive particle laden flows, while it can be precisely forecasted in single-phase flows as a function of the Reynolds number Pope (2001) and the properties of the wall surface (e.g. porosity Breugem et al (2006), Rosti et al (2015), Rosti et al (2018) and elasticity Rosti and Brandt (2017)). The additional complexity of multiphase flows is due to the presence of extra parameters such as the size, shape and elasticity of particles which become relevant. The density difference with respect to the carrier fluid and the solid volume fraction, denoted here with Φ, are also crucial to describe their dynamics. In the present work, we analyze moderate dense suspensions where particles are deformable, fully resolving the fluid–structure interactions thanks to numerical simulations. In recent years, there has been a growth in the number of studies related to deformable objects, with the analysis of a single liquid droplet with constant surface tension, red blood cells enclosed by a biological membrane, and cells with stiff nuclei Loewenberg (1998), Freund et al (2012), Freund (2014), Takeishi et al (2019), Alizad Banaei et al (2017), Rosti and Brandt (2018). Here, we expand the literature of the field (in particular the work in Rosti and Brandt (2018)) by studying a full suspension of deformable particles and by considering the effect of finite inertia (although in the laminar regime), thus studying the interaction of elasticity and inertia for the first time. Furthermore, we study the problem in a Poiseuille flow, thus considering a non-uniform shear rate in the domain with the consequent particle migration and accumulation documented below.

Already in 1911, Einstein (1911) showed that the suspension viscosity increases linearly with the particle volume fraction Φ, when focusing in the limit of vanishing inertia and for dilute suspensions. Banaei et al (1977) and Batchelor and Green (1972) added a second-order correction in Φ, but for higher volume fractions, the viscosity increases faster than a second-order polynomial Stickel and Powell (2005). In these cases, it is necessary to turn to empirical fits such as the Eilers fit Ferrini et al (1979), Zarraga et al (2000), Singh and Nott (2003), Kulkarni and Morris (2008), which matches the rheology of rigid particle suspensions at small Reynolds numbers, for both low and high concentrations. When inertia becomes finite, deviations from the behavior predicted by the different empirical fits start to occur Lin et al (1970), Kulkarni and Morris (2008), an effect that can be related to an increase in the effective volume fraction at intermediate values of Φ Picano et al (2013). Furthermore, at high volume fractions, once friction forces become relevant, other interesting phenomena can also be observed Morris (2018).

The comprehension of the rheology of deformable objects has been a challenging task for many years Roscoe (1967). In this case, the capillary number Ca (the ratio between viscous and elastic forces) is the non-dimensional parameter governing particle deformability. When Ca is low, particles easily recover the equilibrium shape, and deformations are small. These configurations are dominated by elastic forces. In 1932, Taylor first assumed small deformations and showed that, for small Φ, the coefficient of the linear term in Einstein's relation is a function of the ratio between the particle and fluid viscosities. Using perturbative expansions, later analytical studies Cox (1969), Frankel and Acrivos (1970), Choi and Schowalter (1975), Pal (2003) tried to extend the result to higher orders in Φ and Ca, similarly to what done by Batchelor for rigid particles. Only recently, different authors used high-fidelity numerical simulations to analyze such problem Ii et al (2012), Krüger et al (2014), Oliveira and Cunha (2015), Srivastava et al (2016), Matsunaga et al (2016), Rosti and Brandt (2018).

The present work considers the inertial flow of viscous hyperelastic deformable particles, generally used to describe rubber-like substances. This class of constitutive relations can often show non-linear stress–strain curves and have a constitutive equation that is only a function of the current state of deformation. Moreover, if the work done by the stresses during a deformation process is dependent only on the initial and final configurations, the behavior of the material is path independent, and a stored strain energy function or elastic potential can be defined Bonet and Wood (1997). Many studies used this class of constitutive relations to describe particles, capsules, vesicles and even red blood cells Sugiyama et al (2011), Ii et al (2011), Ii et al (2012), Villone et al (2014), Villone et al (2016).

The objective of this paper is twofold: on one hand, we would like to have a better understanding of the physics involved in suspensions of deformable particles by considering the interactions of inertial and elastic effects, while on the other hand, we want to implement a robust numerical tool that can be applied in many applications involving elastic objects in inertial flows. Particle-resolved numerical simulations (full DNS) are becoming increasingly important for analyzing multiscale physical phenomena and providing detailed results, which can be exploited in modeling (i.e. turbulence or interaction modeling). In this context, we hope this work will open new pathways, so that more accurate and reliable closure equations will be derived in the near future. Another aspect to be highlighted is the intrinsic difficulty of the problem under investigation. First of all, the problem is a 4-way coupling interaction, meaning that every single object in the simulation can interact with its neighbors and with the surrounding fluid. Additionally, the resolution required to simulate accurately dense deformable particle suspensions is significantly higher than for rigid objects.

The obtained results are very promising. We show that suspensions of deformable particles exhibit a shear-thinning behavior that is much more sensitive to changes of Ca rather than variations in Re. In particular, deformable particles appear to move toward the center of the channel, and this phenomenon is particularly marked in dilute suspensions. Finally, our computations reveal that the more deformable particles tend to accumulate more in the central region of the channel, as they can modify their shape more easily.

2.1. Governing equations

We consider the flow of a suspension of deformable viscous hyperelastic particles in an incompressible Newtonian viscous fluid. The flow, as in the Poiseuille configuration, is driven by an imposed constant discharge, and the resulting pressure gradient in the stream-wise direction is thus computed accordingly. The number of particles in the domain is summarized by the total solid volume fraction Φ, which is the ratio between the volume occupied by the solid phase and the total volume. The solid suspension is neutrally buoyant, i.e. it has the same density of the carrier fluid $\rho^{s} = \rho^{f} = \rho$, as in many analyses for biological systems Torii et al (2008), Zhao et al (2008). The unstressed reference shape of a particle is a sphere and no external body forces are applied on the continua.

The fluid and solid phase motion is governed by the conservation of momentum, and the incompressibility constraint is enforced in each phase

$$\begin{eqnarray} \frac{\partial u_i^f}{\partial t} + \frac{\partial u_i^f u_j^f}{\partial x_j} = \frac{1}{\rho} \frac{\partial \sigma_{ij}^f}{\partial x_j} \ \ \ \ \ \ \ \ \ \ &\textrm{if }& x \in \Omega^{f}, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \frac{\partial u_i^f}{\partial x_i} = 0 \ \ \ \ \ \ \ \ \ \ &\textrm{if }& x \in \Omega^{f}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \frac{\partial u_i^s}{\partial t} + \frac{\partial u_i^s u_j^s}{\partial x_j} = \frac{1}{\rho} \frac{\partial \sigma_{ij}^s}{\partial x_j} \ \ \ \ \ \ \ \ \ \ &\textrm{if }& x \in \Omega^{s}, \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \frac{\partial u_i^s}{\partial x_i} = 0 \ \ \ \ \ \ \ \ \ \ &\textrm{if }& x \in \Omega^{s}, \end{eqnarray} \tag{ 4 }$$

where the superscripts ^f and ^s denote the fluid and solid phase fields. In the previous balance equations, u_i indicates the ith component of the velocity vector, ρ is the mass density and σ_ij is the Cauchy stress tensor. The latter has a different form in the two phases and characterizes the solid and liquid behavior.

The kinematic and dynamic interactions between the fluid and solid phases are determined by enforcing the continuity of the velocity and traction force at the fluid–structure interface, namely

$$\begin{eqnarray} u_i^f = u_i^s \ \ \ \ \ \ \ \ \ \ &\textrm{on }& \partial \Omega^{f,s}, \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} \sigma_{ij}^f n_j = \sigma_{ij}^s n_j \ \ \ \ \ \ \ \ \ \ &\textrm{on }& \partial \Omega^{f,s}, \end{eqnarray} \tag{ 6 }$$

where n_j denotes the jth component of the normal vector at the interface.

To numerically solve the fluid–structure interaction, it is useful to introduce the so-called one-continuum formulation Tryggvason et al (2007), where only one set of equations is solved over the whole domain $\Omega = \Omega^{f} \cup \Omega^{s}$. This is achieved by introducing a monolithic velocity vector field u valid everywhere. We also introduce a local solid volume fraction function φ, applying the volume averaging procedure Takeuchi et al (2010), Quintard and Whitaker (1994). The latter behaves as an indication function: it is zero in the fluid phase and equal to one in the solid phase. Due to the finite representation of our discrete grid, we have $0 \lt \phi \lt 1$ close to the fluid–solid interface. In this light, we could interpret φ as a smoothed Heaviside function at the grid scale. In particular, the isoline φ = 0.5 represents the interface. Moreover, if we integrate φ over the entire volume, we obtain the global solid volume fraction in the domain, Φ. The one-fluid formulation then reduces to

$$\begin{equation} \sigma_{ij} = \left( 1 - \phi \right) \sigma_{ij}^f + \phi \sigma_{ij}^s. \end{equation} \tag{ 7 }$$

The fluid is assumed to be Newtonian, so that the stress in the fluid phase can be written as

$$\begin{equation} \sigma_{ij}^f = -p \delta_{ij} + 2 \mu^f D_{ij}, \end{equation} \tag{ 8 }$$

where p is the pressure, µ^f the fluid dynamic viscosity, δ_ij the Kronecker delta and D_ij the strain rate tensor, defined as $D_{ij} = \left( \partial u_i/\partial x_j + \partial u_j/\partial x_i \right)/2$. The solid is modeled as an incompressible viscous hyperelastic material, undergoing only the isochoric motion, with constitutive equation

$$\begin{equation} \sigma_{ij}^s = -p \delta_{ij} + 2 \mu^s D_{ij} + G B_{ij}, \end{equation} \tag{ 9 }$$

where the last term is the hyperelastic contribution modeled as a neo-Hookean material, thus satisfying the incompressible Mooney–Rivlin law. In equation (9), µ^s is the solid dynamic viscosity, B_ij the left Cauchy–Green deformation tensor, and G is the modulus of transverse elasticity. In the following, we will assume that the solid and liquid phases have the same dynamic viscosity, so that $\mu^f = \mu^s$.

The set of equations for the solid material can be closed in a purely Eulerian manner by introducing a transport equation for the volume fraction φ:

$$\begin{equation} \frac{\partial \phi}{\partial t} + \frac{\partial u_k \phi}{\partial x_k} = 0, \end{equation} \tag{ 10 }$$

and updating the left Cauchy–Green deformation tensor with the following transport equation

$$\begin{equation} \frac{\partial B_{ij}}{\partial t} + \frac{\partial u_k B_{ij}}{\partial x_k} = B_{kj}\frac{\partial u_i}{\partial x_k} + B_{ik}\frac{\partial u_j}{\partial x_k}. \end{equation} \tag{ 11 }$$

Combining all the previous equations we finally obtain

$$\begin{eqnarray} &&\rho \Big( \frac{\partial u_i}{\partial t} + \frac{\partial u_i u_j}{\partial x_j} \Big) = - \frac{\partial p}{\partial x_i} + 2 \mu^f \frac{\partial D_{ij}}{\partial x_j} + G \frac{\partial B_{ij}}{\partial x_j}, \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} &&\frac{\partial u_i}{\partial x_i} = 0, \end{eqnarray} \tag{ 13 }$$

where B_ij is different from zero only in the presence of a solid particle. When Φ = 0, the governing equations reduce to those for Newtonian fluids.

Zoom In Zoom Out Reset image size

The volume of fluid method in general does not include any repulsive force for particle–particle and particle–wall interactions. Without any kind of repulsive force, if two particles are close enough, they will eventually touch and merge. While this is possible for certain fluids in certain conditions (i.e. bubbles), we do not expect this behavior from solid particles where instead the solid stress will grow abrouptly. Thus, in order to avoid these anomalies, we introduce a subgrid repulsive force, as previously done in Bolotnov et al (2011). The following expression for the force is used Clift et al (2005), both for particle–particle and particle–wall interactions:

$$\begin{equation} f_i = \mu^f V_b r \Big( \frac{a_1}{d} + \frac{a_2}{d^2} \Big) n_i, \end{equation} \tag{ 14 }$$

where µ^f is the fluid dynamic viscosity, V_b the imposed bulk velocity, r the particle radius, d the distance between the two objects, n_i the ith component of the local normal to the surface and a₁ = 550.0, a₂ = 35.0 are two model coefficients. This sub-grid force is only applied when two objects are close enough, i.e. dx ≤ d ≤ 3dx, where dx is the uniform cell size used in the computations. The range and intensity of these interactions are displayed in figure 2. The force is always computed in the middle of the computational cells situated near the interface of two interacting particles, and directed toward the nearest surface. It is worth noting that these forces are applied only in the fluid phase and appear in the momentum balance as localized volume forces. They may also be interpreted as additional local pressure terms, which cannot be directly computed because of the grid resolution, which inevitably becomes inadequate when two particles are about to collide. The same approach has been used to prevent droplet coalescence in s De Vita et al (2019), De Vita et al (2020).

Zoom In Zoom Out Reset image size

2.2. Numerical implementation

2.3. Numerical setup

We first show in figure 3(a) tridimensional visualization of the suspension at Φ = 0.15, Re = 3.75 and Ca = 0.075. It can be noticed how dense the suspension appears already at a volume fraction of 15%. As expected, the particle deformation is a function of the vertical coordinate and it is more evident near the walls where the local shear is higher.

Zoom In Zoom Out Reset image size

The typical behavior of the particle distribution and deformation for the same case can be observed looking at the upper panel of figure 4, where the solid hyperelastic shear stress is shown in a stream-wise-wall-normal cut plane. The hyperelastic shear stress, always null (blue) outside the solid phase (red), appears high when two or more particles interact and are deformed from the unperturbed reference spherical shape. At this capillary number, the shear rate accounts for a modest deformation, which however becomes considerable near the walls where the velocity gradient is maximum. In the bottom panel of figure 4 the span-wise vorticity component is illustrated by contours together with the velocity field (black arrows) in the same plane as in the top panel. The velocity field appears close to the parabolic Poiseuille solution, with the addition of the perturbation induced locally by the particles. The vorticity appears intense near the wall and in the gap regions between interacting particles.

Zoom In Zoom Out Reset image size

The mean stream-wise velocity profiles are plotted in figure 5. They do not differ significantly from the quadratic solution of the Navier–Stokes equation in the absence of particles, and show a slight shear-thinning behavior with a reduction of the centerline velocity and a slight increase near the walls.

Zoom In Zoom Out Reset image size

We display the mean particle volume fraction $\langle \phi \rangle$ (averaged over x, y and time) versus the wall-normal direction z for Ca = 0.3 and different Reynolds numbers in figure 6(a). It can be noted that for the smallest volume fraction considered here, e.g. Φ = 0.025, particles avoid the near-wall region, which appears depleted. Indeed, single deformable particles are known to migrate toward the pipe/channel center in Poiseuille flows, away from high shear regions Alghalibi et al (2190). In dilute cases, particle mutual interactions are rare and particles can migrate toward the channel center more freely, resulting in an empty portion of the channel near the wall. On the contrary, the volume fraction profiles of denser configurations, i.e. Φ≥ 0.05, show that the particles are distributed through the whole channel, with a clear layer near the wall. This is due to the mutual interactions between particles, increasing in number as Φ increases, forcing them to occupy the whole channel. In the figure, it is also evident that the maximum of the volume fraction profiles, denoting wall layering, moves closer to the wall as Φ increases because of the augmented mutual interaction between the particles. The particle distribution appears weakly affected by changes in the Reynolds number in the range considered here. A near-wall region depleted of particles is noteworthy. Indeed, this result implies that, not only the wall-region has to be properly modeled to account for particle layering, but also the bulk of the channel can have a non-uniform distribution of particles, thus making the problem fully dependent not only on the local shear-rate as for rigid particles, but also on the local volume fraction.

Zoom In Zoom Out Reset image size

Figure 6(b) illustrates the local mean volume fraction at a fixed Re = 3.75 and different Ca. As Ca increases, it can be seen that the particles tend to show a peak slightly further away from the wall. This is because the more deformable particles are more sensitive to the shear rate variations (shear-induced migration), and they can modify their shape more easily to accommodate more particles in the channel center, thus leaving empty a bigger portion of the channel near the wall.

Stream-wise and wall-normal root-mean-square velocity profiles are plotted in figure 7 for different volume fractions at a fixed capillary number Ca = 0.009. It can be observed that the former have higher intensities and present some peaks in correspondence with the particle layers, while the latter are less intense and smoother. Both the fluctuation components increase with the particle volume fractions, slightly increase with Re and are only marginally affected by Ca (not shown here).

Zoom In Zoom Out Reset image size

In order to investigate the rheology of the suspension, in figure 8 we select the simulation parameters to mimic experiments in a fixed geometry and same particle material, while changing the bulk shear rate $\dot{\gamma} = 3V_b/2~{\textrm h}$ (i.e. varying the flow rate in the channel). In terms of dimensionless parameters, this corresponds to consider cases at constant Re/Ca ratio and varying the non-dimensional shear rate $\widetilde{\dot{\gamma}} = \sqrt{\textrm{We}} = \sqrt{\textrm{Re} \, \textrm{Ca}}$. The parameter $\widetilde{\dot{\gamma}}$ can be thought of as the ratio between the convective time scale in the flow and the material deformation time scale. In the figure we show the suspension relative viscosity µ/µ^f , computed as a cross-sectional average, versus the non-dimensional shear rate $\widetilde{\dot{\gamma}}$. The relative viscosity proves to be a decreasing function of the shear rate, which confirms that the suspension has a shear-thinning behavior, as expected for suspensions of this type of deformable particles. The denser the suspension, the more intense this behavior.

Zoom In Zoom Out Reset image size

According to the definition of the Weber number We = Re Ca, both the capillary and Reynolds numbers may contribute to this shear-thinning. We analyze therefore the dependency of the suspension relative viscosity on these two parameters by fixing one and varying the other. In figure 9 we plot the relative viscosity versus Re at constant Ca = 0.3 in the left panel, and versus Ca at constant Re = 3.75 in the right panel, for different volume fractions Φ. The particle Reynolds number varies between 0.45 and 7.5, while the capillary number between 0.009 and 0.3. The data indicate that the relative viscosity weakly depends on Re, while it appears a decreasing function of Ca. Hence, for the range of parameters investigated in this study, the viscosity depends more on the capillary than on the Reynolds number.

Zoom In Zoom Out Reset image size

To gain further insight into the suspension behavior, we plot the same data as explicit functions of the volume fractions in figure 10. It can be noted that cases at different Reynolds numbers (left panel) cluster together and do not show any significant difference, while the effective viscosity curves vary as Ca changes (right panel). In particular, the effective viscosity decreases as Ca increases and this effect is more pronounced at high Φ. In figure 10 we also plot the Eilers fit, which is valid for rigid spherical particles in the inertialess limit, both for dense and dilute suspensions:

$$\begin{equation} \frac{\mu}{\mu^{f}} = \Big[ 1 + B_E \frac{\Phi}{1 - \Phi / \Phi_{m}} \Big]^2. \end{equation} \tag{ 15 }$$

In the previous equation $\Phi_m = 0.58-0.63$ is the geometrical maximum packing, and $B_E = 1.25-1.7$ a fitting coefficient. Φ_m  = 0.58 and B_E  = 1.7 appear to be good values for our simulations. It can be remarked in figure 10(b) that the most rigid case essentially follows the Eilers fit, especially at high volume fractions. However, interestingly for small volume fraction, e.g. Φ = 0.025, the Eilers fit does not provide a good approximation and the suspension viscosity is closer to that of the unladen case. We explain this peculiar behavior by the shear-induced migration of deformable particles. In particular, deformable particles tend to migrate toward the middle of the channel, where the shear rate vanishes. In other words, to avoid deformation, particles tend to avoid the high shear rate region near the wall and stay near the centerline. This is a direct consequence of what previously observed in figure 6, further showing that simple correlations extracted from a Couette geometry may be not enough to fully describe more complex flow configurations such as the present one.

Zoom In Zoom Out Reset image size

The elastic nature of the particles is evaluated here in terms of the first and second normal stress differences, defined as:

$$\begin{eqnarray} \mathcal{N}_1 & = & \langle \langle \sigma_{22} - \sigma_{33} \rangle \rangle = \langle \langle \phi G B_{22} + 2\mu^f \frac{\partial v}{\partial y} - \phi G B_{33} - 2\mu^f \frac{\partial w}{\partial z} \rangle \rangle, \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} \mathcal{N}_2 & = & \langle \langle \sigma_{33} - \sigma_{11} \rangle \rangle = \langle \langle \phi G B_{33} + 2\mu^f \frac{\partial w}{\partial z} - \phi G B_{11} - 2\mu^f \frac{\partial u}{\partial x} \rangle \rangle. \end{eqnarray} \tag{ 17 }$$

These two quantities are a measure of the anisotropy of the flow and are important when characterizing the non-Newtonian behavior of fluids, since together with the effective viscosity can provide information of the full stress tensor. Furthermore, when considering particle suspensions, they are usually associated to the particle migration. $\mathcal{N}_1$ and $\mathcal{N}_2$ are reported in figure 11 as a function of the Reynolds number for different volume fractions at the largest value of Ca considered, and in figure 12 as a function of the volume fraction for different values of Ca. When the capillary number is fixed and sufficiently large, $\mathcal{N}_1$ is positive, $\mathcal{N}_2$ negative, and their ratio $-\mathcal{N}_1/\mathcal{N}_2$ larger than unity. This is consistent with what found for other elastic fluids, such as polymeric solution Shahmardi et al (2019), capsule Matsunaga et al (2016) and fiber suspensions Banaei et al (2020). We observe that, $\mathcal{N}_1$ grows as the Reynolds number increases, while $\mathcal{N}_2$ reduces, thus enhancing their difference. Very different results are observed when the capillary number is varied. Indeed, when Ca is small (almost rigid particles) $\mathcal{N}_1$ is actually negative and becomes positive only for Ca≈ 0.1. This is consistent with what typically found when considering rigid particle suspensions Kulkarni and Morris (2008). Also, $\mathcal{N}_2$ assumes an opposite sign when the capillary number is small, in the limit of rigid particles. In general, a negative value for $\mathcal{N}_1$ can be interpreted as a signature of the importance of particle interactions and contacts, which is dominant for the rigid ones but negligible for the deformable ones, when $\mathcal{N}_1$ even changes sign. Also, the increase of the first normal stress difference with Re and Φ for the deformable cases indicates a stronger tendency of the fluid to displace the particles towards the centerline Banaei et al (2020).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The increase of the relative viscosity is linked to an increase in the stress across the channel. The total shear stress σ₂₃ can be decomposed into the sum of the fluid $\sigma_{23}^{f}$ and particle stress $\sigma_{23}^{s}$, as in equation (7). The mean fluid stress is the sum of the viscous and Reynolds stresses, while the solid one is the sum of the viscous, Reynolds, hyperelastic and collision stresses Rosti and Brandt (2018), Rosti et al (2019), Picano et al (2015)

$$\begin{align} \langle\sigma_{23}\rangle = &\mathop{\underbrace{\langle\left( 1 - \phi \right) \left( \mu^f \frac{\partial v}{\partial z} - \rho v^{^{\prime}}w^{^{\prime}} \right)}\rangle}\limits_\textrm{fluid~stress}\nonumber\\ & + \mathop{\underbrace{\langle\phi \left( \mu^f \frac{\partial v}{\partial z} - \rho v^{^{\prime}}w^{^{\prime}} + G B_{23} \right) + \int_{0}^{z} f_{y} dz}\rangle}\limits_\textrm{solid~stress}. \end{align} \tag{ 18 }$$

If we consider the monolithic formulation independent of the actual phase, the total shear stress can also be decomposed into the sum of the total viscous, total Reynolds, hyperelastic and collision stresses in the following way

$$\begin{equation} \langle\sigma_{23}\rangle = \mathop{\underbrace{\langle\mu^f \frac{\partial v}{\partial z}}\rangle}\limits_\textrm{viscous~stress} + \mathop{\underbrace{\langle- \rho v^{^{\prime}}w^{^{\prime}}}\rangle}\limits_\textrm{Reynolds~stress} + \mathop{\underbrace{\langle G B_{23}}\rangle}\limits_\textrm{elastic~stress} + \mathop{\underbrace{\langle\int_{0}^{z} f_{y} dz}\rangle}\limits_\textrm{collision~stress}. \end{equation} \tag{ 19 }$$

In particular, the total shear stress $\langle\sigma_{23}\rangle$ in a plane channel flow is a linear function, equal to the total wall shear stress at the wall and null in the centerline. All the shear stress contributions are averaged and displayed in figure 13 as functions of the wall-normal distance z for the cases at Re = 3.75, Ca = 0.3 and volume fractions Φ = 0.05, 0.1 and 0.15. The stresses are normalized by the total wall value. In the figures, the light gray, dark gray and black color represent the viscous, Reynolds, and particle stresses respectively. The collision stress contribution is negligible in all the cases considered here, thus not shown in the plots for greater clarity.

Zoom In Zoom Out Reset image size

As expected, the viscous stress is the only not vanishing component at the wall. For all cases, since we are dealing with laminar suspension flows, the viscous stress is the dominant contribution. It however decreases in correspondence of the location of the maximum particle concentrations (see figure 6), where the particle stress reaches its maximum. As the volume fraction is increased, the viscous stress becomes smaller, which is compensated by an increase in the particle stress contribution, which eventually accounts for up to 50% of the total shear stress in the central region of the channel at Φ = 0.15. Concerning the Reynolds stresses, it can be observed that they are almost negligible in every case shown here, rising slightly as Φ and Ca are increased, with a maximum located at around z = 0.3 h.

In figure 14 we summarize the stress budgets for fixed Re = 3.75 and Ca = 0.009 by showing the volume-averaged contributions of all the components of the total shear stress. Viscous (light gray), Reynolds (dark gray), and particle stresses (black) are plotted both as percentages of the total stress and as bulk values. Comparing the budgets at different volume fractions, we observe that as the volume fraction increases, the bulk value of the viscous stress remains essentially the same, while the particle-induced stress increases, becoming of the same order of magnitude as the viscous stress at Φ = 0.15. The Reynolds shear stress, although slightly increasing, remains secondary.

Zoom In Zoom Out Reset image size

Finally, we display in figure 15 the budgets for different Re at fixed Ca and vice versa. At fixed Ca, when increasing Re we note that the particle stress increases in importance with respect to the viscous stress when inertia becomes more relevant. Concerning the increase of Ca at fixed Re, the total amount of the viscous stress is essentially unchanged, but the particle stress amount reduces, so that, in relative terms, the viscous stress becomes more dominant.

Zoom In Zoom Out Reset image size

We studied suspensions of deformable particles to characterize the effect of elasticity in the presence of finite inertia and non-uniform shear rate. The flow is characterized by three non-dimensional parameters: the solid volume fraction, the Reynolds and the capillary numbers. The suspensions proved to have a shear-thinning behavior as the relative viscosity happens to be a decreasing function of the shear rate even for finite Reynolds numbers. The denser the suspension, the more this behavior is emphasized.

We started examining the averaged mean velocity, particle distribution, and velocity fluctuation profiles along the wall-normal coordinate. To begin with, the mean velocity is very close to the parabolic solution of the Navier–Stokes equations but exhibits a blunt shape typical of shear-thinning fluids. It does not vary considerably in the cases under investigation. The particle distribution profiles show some concentration peaks at off-center locations, which point to a layered structure. These profiles depend mostly on the total volume fraction and the capillary number. On the other hand, the velocity fluctuations are more sensitive to changes of Φ and Re rather than Ca.

A separate analysis of the dependency of the suspension relative viscosity on the Reynolds and capillary numbers was carried out by fixing one parameter and varying the other. The data indicate that the relative viscosity weakly depends on the Reynolds number, while it appears a decreasing function of the deformability parameter. Hence, for the range of parameters investigated in this study, the viscosity is much more dependent on the capillary than on the Reynolds number.

We showed that for small volume fraction the Eilers fit does not provide a good approximation of the suspension viscosity, which is closer to that of the unladen case. We explained this peculiar behavior by the fact that deformable particles tend to migrate toward the middle of the channel, where the shear rate is low. In other words, to avoid deformation, particles tend to stay away from the high shear rate region near the wall and cluster toward the centerline of the channel. Because of the small value of the shear rate in the core region, the suspended phase does not contribute to the rise in the effective viscosity as in the rigid case, where no migration occurs. This phenomenon is particularly marked for very dilute suspensions, where a whole region near the wall is empty of particles at stationarity. When the volume fraction is increased, this near-wall region is gradually occupied. This is due to the mutual interactions between particles, increasing in number as the volume fraction increases, forcing them to occupy the whole channel in dense configurations. Moreover, as the capillary number increases, this deformation-driven migration is more noticeable. Indeed, the more deformable particles are more sensitive to shear rate variations, and they can modify their shape more easily to accommodate a greater number of particles in the channel center. This behavior makes the suspension highly non-uniform across the domain, demonstrating the difference with what usually measured and observed in a more standard Couette rheological configuration.

The analysis of the first and second normal stress differences allowed to better evaluate the elastic properties of the particles, and their dependency on the Reynolds and capillary numbers was investigated. We showed that when the capillary number is fixed and sufficiently large, $\mathcal{N}_1$ is positive, $\mathcal{N}_2$ negative, while when Ca is small $\mathcal{N}_1$ is actually negative and becomes positive only for Ca≈ 0.1. Also, $\mathcal{N}_2$ assumes an opposite sign when the capillary number is small, in the limit of rigid particles. The change of sign of the normal stress differences marks the change of behavior of the suspension, from a rigid-like behavior dominated by particle–particle interactions to the deformable one dominated by the particle deformation and consequent tendency to migrate towards the centerline.

Finally, the total stress budgets revealed that the relative particle-induced stress contribution grows with the volume fraction and Reynolds number, and reduces as a function of the particle deformability. For all cases, since we are dealing with laminar suspension flows, the viscous stress is the dominant contribution and, as expected, is the only not vanishing component at the wall. It becomes smaller as the volume fraction is increased, but this is compensated by a larger particle stress contribution. Concerning the Reynolds stresses, it can be noted that they are almost negligible in every case under investigation.

The present analysis demonstrated the importance of studying suspension of deformable objects at finite inertia and for non-uniform shear rates. Indeed, we show that results can be different from what observed in Couette flows due to non-uniformity of the shear which, coupled with particle defomation, brings to a suspension behavior which is dependent not only on the shear rate, but also on the local particle concentration. This should be carefully considered when proposing macroscopic models for suspensions of deformable objects.

The authors acknowledge computer time provided by SNIC (Swedish National Infrastructure for Computing) and by the Scientific Computing section of Research Support Division at OIST.