2 Mathematical formulation

The liquid flows through a channel between two parallel plates and then is extruded into the ambient atmosphere. A sketch of the geometry of the problem is shown in figure 1, where $H$ is half of the channel gap. The free-surface height $h$ is unknown a priori and depends on the flow parameters. The lengths of the inlet channel and free-surface region were chosen to be $75H$ and $25H$ , respectively. The flow is assumed to be symmetric with respect to the $x_{2}=0$ plane. Therefore, only the upper half of the channel is considered. It is worth mentioning that particle migration does not reach a fully developed condition at the end of the upstream channel. As discussed by Nott & Brady (Reference Nott and Brady1994), the channel length necessary to obtain a fully developed particle concentration distribution is very high. However, the geometry used here is long enough to lead to a strongly non-uniform particle concentration field and consequently to viscosity gradients that change the flow behaviour, as discussed in the results section.

Figure 1. A schematic diagram of the geometry of the problem (not to scale).

Because of the typical scales involved in the problem, i.e. the particle diameter is much smaller than the channel height, the suspension can be treated as an equivalent continuum liquid in the flow domain. Under this circumstance, the flow is governed by incompressible mass conservation and momentum conservation equations, namely

and

In (2.1) and (2.2), $\boldsymbol{u}$ is the velocity field vector, $\unicode[STIX]{x1D70C}$ is the suspension density, $\unicode[STIX]{x1D72E}$ is the suspension stress tensor and $\boldsymbol{g}$ is the gravitational field vector. The stress tensor is split as $\unicode[STIX]{x1D72E}=-p\unicode[STIX]{x1D644}+\unicode[STIX]{x1D749}$ , where $p$ is the pressure field, $\unicode[STIX]{x1D644}$ is the unit tensor and $\unicode[STIX]{x1D749}$ is the extra-stress tensor field.

We assume that the extra-stress is a purely viscous stress that obeys Newton’s law with a concentration-dependent viscosity, such that $\unicode[STIX]{x1D749}=\unicode[STIX]{x1D702}\dot{\unicode[STIX]{x1D738}}$ , where $\unicode[STIX]{x1D702}=\unicode[STIX]{x1D702}(\unicode[STIX]{x1D719})$ is the suspension viscosity, $\unicode[STIX]{x1D719}$ is the local particle concentration and $\dot{\unicode[STIX]{x1D738}}=\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u}^{\text{T}}$ is the rate-of-strain tensor, in which the superscript T denotes the transpose operation. The suspension viscosity is given by a differential-effective-medium-based model recently proposed by Santamaría-Holek & Mendoza (Reference Santamaría-Holek and Mendoza2010), such that

where $\unicode[STIX]{x1D702}_{0}$ is the viscosity of the Newtonian solvent liquid, $[\unicode[STIX]{x1D702}]$ is the so-called intrinsic viscosity and $c=(1-\unicode[STIX]{x1D719}_{c})/\unicode[STIX]{x1D719}_{c}$ is a crowding factor that guarantees that the particles cannot occupy all the volume of the sample due to geometric restrictions, in which $\unicode[STIX]{x1D719}_{c}$ is the critical particle concentration at which the suspension loses its fluidity, i.e.  $\unicode[STIX]{x1D702}\rightarrow \infty$ as $\unicode[STIX]{x1D719}\rightarrow \unicode[STIX]{x1D719}_{c}$ . For suspensions of hard spheres, it follows that $[\unicode[STIX]{x1D702}]=2.5$ and $\unicode[STIX]{x1D719}_{c}=0.637$ . It is worth mentioning that (2.3) recovers the classical result of Einstein (Reference Einstein1911) at the limit of infinitely dilute suspensions, and provides good results when compared to experimental data over the entire concentration range.

Because of irreversible particle–particle interactions that occur in non-uniform shear flows, the suspended particles tend to migrate to preferred regions in the flow domain. Shear-induced particle migration is described by the diffusive flux model developed by Phillips et al. (Reference Phillips, Armstrong, Brown, Graham and Abbott1992). The model states a typical convection–diffusion equation accounting for particle transport in the flow, such that

where $\boldsymbol{N}_{\unicode[STIX]{x1D719}}$ is the total diffusive flux of particles due to different mechanisms. Initially, Phillips et al. (Reference Phillips, Armstrong, Brown, Graham and Abbott1992) recognized that particle migration arises from two different mechanisms, namely gradients in shear rate and gradients in suspension viscosity. They used the theory of Leighton & Acrivos (Reference Leighton and Acrivos1987) to develop mathematical expressions for these two fluxes, as follows:

and

where $\boldsymbol{N}_{c}$ and $\boldsymbol{N}_{\unicode[STIX]{x1D702}}$ are the diffusive flux of particles because of gradients in shear rate and gradients in suspension viscosity, respectively. In (2.5) and (2.6), $\dot{\unicode[STIX]{x1D6FE}}=|\dot{\unicode[STIX]{x1D738}}|$ is the local shear rate (i.e. the second invariant of the rate-of-deformation tensor), $a$ is the particle radius, and the parameters $K_{c}$ and $K_{\unicode[STIX]{x1D702}}$ are diffusion-like coefficients inherent to the model that should be found from experimental results. We set $K_{c}=0.41$ and $K_{\unicode[STIX]{x1D702}}=0.62$ , as found by Phillips et al. (Reference Phillips, Armstrong, Brown, Graham and Abbott1992) in pressure-driven flows through tubes and cylindrical Couette flows. We also consider an additional diffusive flux of particles because of the curvature of the streamlines in the flow, as introduced by Krishnan et al. (Reference Krishnan, Beimfohr and Leighton1996) and later adapted by Kim et al. (Reference Kim, Lee and Kim2008) to a frame-independent formulation. The mathematical expression for this flux is

where $\unicode[STIX]{x1D705}$ is the local radius of curvature of a streamline and $\boldsymbol{n}_{s}$ is the unit normal vector in the radially outward direction of a curved streamline. Following Kim et al. (Reference Kim, Lee and Kim2008), the additional diffusion-like coefficient is defined as $K_{\unicode[STIX]{x1D705}}=K_{c}$ . It is worth mentioning that we are considering suspensions of neutrally buoyant particles, so that the effects related to sedimentation because of gravity are neglected, and the suspension density is equal to the density of the Newtonian solvent liquid.

The boundary conditions used to solve (2.1), (2.2) and (2.4) are given below, in which $\boldsymbol{n}$ and $\boldsymbol{t}$ are local unit normal and tangent vectors to the boundaries, respectively.

1. (i) Synthetic inflow plane: At the inflow plane, we assume a uniform particle concentration, $\unicode[STIX]{x1D719}=\bar{\unicode[STIX]{x1D719}}$ , where $\bar{\unicode[STIX]{x1D719}}$ is the suspension bulk concentration, and a fully developed, parabolic velocity profile, so that

   (2.8) $$\begin{eqnarray}\boldsymbol{u}=\frac{3}{2}V\left[1-\left(\frac{x_{2}}{H}\right)^{2}\right]\hat{\boldsymbol{e}}_{1},\end{eqnarray}$$
   where $V=Q/H$ is the flow mean velocity and $Q$ is the imposed flow rate (per unit width).

2. (ii) Symmetry plane: Because of symmetry, we have that $\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{u}=0$ , $\boldsymbol{t}\boldsymbol{\cdot }(\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D72E})=0$ and $\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{N}_{\unicode[STIX]{x1D719}}=0$ along the flow centreline.

3. (iii) Synthetic outflow plane: At the outflow plane, we assume a fully developed velocity profile with an imposed pressure, $\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=0$ and $p=0$ . Moreover, there is no diffusive flux of particles, such that $\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{N}_{\unicode[STIX]{x1D719}}=0$ .

4. (iv) Free surface: At the free surface, the shear stress vanishes and there is no liquid flux across the air–liquid interface, so that $\boldsymbol{t}\boldsymbol{\cdot }(\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D72E})=0$ and $\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{u}=0$ . Moreover, the liquid traction should balance the sum of the pressure in the external gas and the capillary pressure induced by the curvature of the free surface. The stress jump across the interface is written as $\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D72E}=(-p_{g}+\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705})\boldsymbol{n}$ , where $p_{g}=0$ is the external gas pressure, $\unicode[STIX]{x1D70E}$ is the suspension surface tension and $\unicode[STIX]{x1D705}$ is the curvature of the free surface. The liquid surface tension is assumed to be constant and independent of the local particle concentration at the interface, such that there are no Marangoni effects (Chandrasekhar Reference Chandrasekhar1961). Finally, we also assume that there is no adsorption or desorption of particles in the free surface, which leads to $\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{N}_{\unicode[STIX]{x1D719}}=0$ .

5. (v) Static contact line: The contact line is pinned at the corner of the channel exit.

A dimensional analysis of the problem suggests that it is convenient to introduce three dimensionless numbers, which are combinations of the various macroscopic model parameters. These dimensionless numbers are: (i) the macroscopic Reynolds number, $Re=\unicode[STIX]{x1D70C}Q/\bar{\unicode[STIX]{x1D702}}$ ; (ii) the capillary number, $Ca=\bar{\unicode[STIX]{x1D702}}V/\unicode[STIX]{x1D70E}$ ; and (iii) the ratio of particle radius to channel gap, $\unicode[STIX]{x1D701}=a/H$ . Here, $\bar{\unicode[STIX]{x1D702}}=\unicode[STIX]{x1D702}(\bar{\unicode[STIX]{x1D719}})$ is the average suspension viscosity evaluated with the bulk concentration $\bar{\unicode[STIX]{x1D719}}$ , i.e. neglecting particle migration. All simulations were performed at $Re=0$ , $a=100~\unicode[STIX]{x03BC}$ m and $\unicode[STIX]{x1D701}=10^{-1}$ , so that the only parameters that are varied in this work are $Ca$ and $\bar{\unicode[STIX]{x1D719}}$ .

3 Numerical method

Because of the free surface, the domain of the extrudate swell flow analysed here is not known a priori. The domain is part of the solution and depends on the flow parameters. In order to solve this free-boundary flow by means of standard techniques for fixed domain problems, the set of equations and boundary conditions posed in the unknown domain (i.e. the physical domain) have to be transformed to an equivalent problem defined in a known domain (i.e. a reference computational domain). This transformation is made by a mapping $\boldsymbol{x}=\boldsymbol{x}(\unicode[STIX]{x1D743})$ that connects the two domains. The mapping procedure used here is based on the elliptic mesh generation method presented by Christodoulou (Reference Christodoulou1990) and de Santos (Reference de Santos1991). In this case, the inverse mapping satisfies an elliptic equation identical to those encountered in diffusion transport with variable diffusion coefficients,

where $\unicode[STIX]{x1D743}$ is the position vector in the computational domain and $\boldsymbol{{\mathcal{D}}}$ is a symmetric positive definite second-order tensor of diffusion-like coefficients used to control element spacing in the reference domain (Benjamin Reference Benjamin1994). Equation (3.1) is usually referred to as the mesh generation equations. Boundary conditions are also needed to solve the second-order partial differential equations of mesh generation. The solid boundaries and the synthetic inflow and outflow planes are defined by functions of the coordinates in the physical domain, and stretching functions are used to distribute the nodes along the boundaries. At the static contact line, either the position of the contact line is specified, or the angle between the unit normal vector outwards from the free surface and the unit normal vector outwards from the solid boundary is fixed, and the contact line lies on the solid boundary. Here, we kept the position of the contact line fixed at $\boldsymbol{x}/H=(0,1)$ , so that the contact angle is part of the solution of the problem. The position of the free surface is implicity defined by the kinematic condition $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}=0$ , such that there is no liquid flux across the air–liquid interface. This version of the elliptic mesh generation method has been widely used in numerical simulations of free-surface flows, such as coating flows of Newtonian liquids (Carvalho & Scriven Reference Carvalho and Scriven1997; Carvalho & Kheshgi Reference Carvalho and Kheshgi2000; Tjiptowidjojo & Carvalho Reference Tjiptowidjojo and Carvalho2011), non-Newtonian polymer solutions (Pasquali & Scriven Reference Pasquali and Scriven2002; Romero et al. Reference Romero, Suszynski, Scriven and Carvalho2004; Romero, Scriven & Carvalho Reference Romero, Scriven and Carvalho2006; Bajaj, Prakash & Pasquali Reference Bajaj, Prakash and Pasquali2008) and particle suspensions (Campana et al. Reference Campana, Silva and Carvalho2017; Siqueira et al. Reference Siqueira, Rebouças and Carvalho2017b ), stretching of liquid bridges (Dodds, Carvalho & Kumar Reference Dodds, Carvalho and Kumar2009, Reference Dodds, Carvalho and Kumar2011, Reference Dodds, Carvalho and Kumar2012) and predictions of dynamic wetting failure (Vandre, Carvalho & Kumar Reference Vandre, Carvalho and Kumar2012, Reference Vandre, Carvalho and Kumar2013, Reference Vandre, Carvalho and Kumar2014; Liu et al. Reference Liu, Vandre, Carvalho and Kumar2016).

Numerical approximations of the transport equations and the mesh generation equations were obtained by using the DEVSS-TG/SUPG mixed finite element method (Guénette & Fortin Reference Guénette and Fortin1995; Szady et al. Reference Szady, Salamon, Liu, Armstrong and Brown1995; Pasquali & Scriven Reference Pasquali and Scriven2002). The DEVSS-TG formulation treats the velocity gradient as an independent variable with a traceless continuous interpolation. The interpolated velocity gradient is then used to compute the viscous stress with a discrete adaptive split. The independent variables of the problem are written as a linear combination of a finite number of basis functions. Lagrangian biquadratic functions represent the velocity field, particle concentration and mesh position; Lagrangian bilinear functions represent the interpolated velocity gradient; and linear discontinuous functions are used for the pressure field. Galerkin weighting functions are used in the residual equations of mass conservation, momentum conservation, mesh position and interpolated velocity gradient, and the streamline-upwinding Petrov–Galerkin formulation is applied in the particle transport equation. In addition to that, the singularity of the diffusive flux model was treated with a modified version of the non-local stress contribution to the shear rate presented by Rebouças et al. (Reference Rebouças, Siqueira, de Souza Mendes and Carvalho2016). The result is a large, sparse set of coupled, nonlinear algebraic equations, which was solved using Newton’s method with a numerical Jacobian matrix obtained with a second-order central difference scheme. The tolerance on the $L^{2}$ -norm of the global residual vector and Newton’s update was set to $10^{-6}$ . At each iteration, the linear system was solved with a frontal solver based on the LU factorization method (Duff, Erisman & Reid Reference Duff, Erisman and Reid1989), and solutions at different parameters were obtained with a first-order arclength continuation (Bolstad & Keller Reference Bolstad and Keller1986).

We performed a grid-independence study for a representative combination of the governing parameters to ensure that the numerical results are independent of the mesh refinement. Three different meshes were examined, and the number of elements, nodes and degrees of freedom (DoFs) of each one of them is specified in table 1. The accuracy of our approximations was assured by computing the shape of the free surface and the particle concentration distribution at the free-surface outflow plane ( $x_{1}/H=25$ plane). These results are presented in figure 2. In this test, we have used $Ca=2.5$ and $\bar{\unicode[STIX]{x1D719}}=0.20$ . As can be seen, the three meshes yield quite close results. Mesh 3, illustrated in figure 3, was selected and employed for all cases investigated.

Figure 2. Mesh independence test: (a) position of the free surface; (b) particle concentration distribution at the outflow plane ( $x_{1}/H=25$ plane). Simulations at $Ca=2.5$ and $\bar{\unicode[STIX]{x1D719}}=0.20$ .

Figure 3. Representative mesh used in the simulations (not to scale).

Figure 4. Particle concentration along the flow centreline ( $x_{2}/H=0$ plane) for: (a)  $\bar{\unicode[STIX]{x1D719}}=0.10$ ; (b) $\bar{\unicode[STIX]{x1D719}}=0.20$ ; (c) $\bar{\unicode[STIX]{x1D719}}=0.30$ ; (d) $\bar{\unicode[STIX]{x1D719}}=0.40$ . Simulations at $Ca=1$ .

4 Results and discussion

First, we investigate particle migration by examining the particle concentration along the flow centreline. Figure 4 shows the results for suspensions at different bulk concentrations at $Ca=1$ . As expected, particle concentration along the centreline increases as the suspension flows inside the channel. Indeed, particles migrate from regions of high shear rate near the wall to regions of low shear rate near the symmetry plane. Particle concentration reaches a maximum value close to the end of the channel. For $\bar{\unicode[STIX]{x1D719}}=0.10$ , the maximum concentration is $\unicode[STIX]{x1D719}_{m}\approx 0.35$ and occurs at $x_{1}/H\approx -1.5$ , and for $\bar{\unicode[STIX]{x1D719}}=0.40$ it is $\unicode[STIX]{x1D719}_{m}\approx 0.58$ and occurs at $x_{1}/H\approx -8.1$ . Just before the exit of the channel, the centreline concentration shows a sudden decrease, which is a consequence of the transition from the bounded flow inside the channel to the free-surface flow outside the channel. At the intersection between the channel and the free surface, the shear rate undergoes a rapid variation from a maximum value at the channel wall to zero at the air–liquid interface. Therefore, the free-surface induces a redistribution in the particle concentration field near the channel exit, pushing the particles towards the air–liquid interface and thereby decreasing the centreline concentration. However, for $x_{1}/H\gtrsim 1$ , the centreline concentration along the free surface remains constant and higher than the suspension bulk concentration. For $\bar{\unicode[STIX]{x1D719}}=0.10$ , the concentration along the free-surface centreline is $\unicode[STIX]{x1D719}_{cfs}\approx 0.19$ , and for $\bar{\unicode[STIX]{x1D719}}=0.40$ it is $\unicode[STIX]{x1D719}_{cfs}\approx 0.50$ .

The constant concentration along the free-surface centreline shows that shear-induced particle migration is almost negligible in the flow under the free surface. The flow in that region is bounded by the flow symmetry plane and the air–liquid interface, and since the shear rate vanishes along these two boundaries, the gradient in shear rate is almost zero, i.e. the flow approaches a plug flow condition. The constant centreline concentration under the free surface is greater than the suspension bulk concentration, $\unicode[STIX]{x1D719}_{cfs}>\bar{\unicode[STIX]{x1D719}}$ , which indicates the existence of gradients in particle concentration in the extrudate. Thus, particle migration in this region arises only from gradients in suspension viscosity and curved streamlines. However, as discussed by Nott & Brady (Reference Nott and Brady1994), the diffusive mechanisms are relatively weak, so that particle migration under the free surface is very slow. Figure 5 displays the particle concentration distribution at the outflow plane. In all cases, particle concentration is higher near the flow centreline, and drastically decreases towards the air–liquid interface. For $\bar{\unicode[STIX]{x1D719}}=0.10$ , the interface concentration is $\unicode[STIX]{x1D719}_{i}\approx 0.06$ , and for $\bar{\unicode[STIX]{x1D719}}=0.40$ it is $\unicode[STIX]{x1D719}_{i}\approx 0.27$ . Figure 6 shows the particle concentration field in the entire flow domain, highlighting the non-uniform distribution of particle concentration inside the channel and under the free-surface region.

Figure 5. Particle concentration distribution at the free-surface outflow plane ( $x_{1}/H=25$ plane) for: (a) $\bar{\unicode[STIX]{x1D719}}=0.10$ ; (b) $\bar{\unicode[STIX]{x1D719}}=0.20$ ; (c) $\bar{\unicode[STIX]{x1D719}}=0.30$ ; (d) $\bar{\unicode[STIX]{x1D719}}=0.40$ . Simulations at $Ca=1$ .

Figure 6. Particle concentration field in the flow domain for: (a)  $\bar{\unicode[STIX]{x1D719}}=0.10$ ; (b) $\bar{\unicode[STIX]{x1D719}}=0.20$ ; (c) $\bar{\unicode[STIX]{x1D719}}=0.30$ ; (d) $\bar{\unicode[STIX]{x1D719}}=0.40$ . Simulations at $Ca=1$ .

In all plots presented in figure 5 the $x_{2}$ coordinate in the free jet was normalized by the free jet final height, $h$ . However, the free-surface height and the shape of the free surface itself depend strongly on the flow parameters, as can be seen in figure 6. For a circular die at the limit of small Reynolds number, the interface swells to approximately 13 % of the tube diameter in the case of Newtonian liquids. In the case of polymeric solutions, the extrudate swell can be as high as 300 %. This difference in the extrudate diameter is related to the presence of normal stress differences in polymeric liquids (Tanner Reference Tanner1985; Bird, Armstrong & Hassager Reference Bird, Armstrong and Hassager1987). Normal stress differences in the flow of non-colloidal suspensions can be computed if the suspension balance model of Nott & Brady (Reference Nott and Brady1994) and Morris & Boulay (Reference Morris and Boulay1999) is used. However, the diffusive flux model together with the constitutive equation used here does not show normal stress difference, i.e.  $N_{1}=N_{2}=0$ . Therefore, the predicted variation of the extrudate swell is related to the non-uniform particle distribution in the flow due to shear-induced migration. Figure 7 compares the shape of the free surface of suspensions with different bulk concentrations (but with the same average bulk viscosity) and a Newtonian liquid with the same average bulk properties. The extrudate final height decreases considerably as the suspension concentration increases. At low concentrations, the extrudate height is smaller than that of the Newtonian flow. For the Newtonian liquid, $h/H\approx 1.13$ , and for $\bar{\unicode[STIX]{x1D719}}=0.10$ and $\bar{\unicode[STIX]{x1D719}}=0.20$ , we found $h/H\approx 1.11$ and $h/H\approx 1.05$ , respectively. In turn, at high concentrations, the extrudate height is smaller than the channel gap; at $\bar{\unicode[STIX]{x1D719}}=0.30$ , $h/H\approx 0.99$ , and at $\bar{\unicode[STIX]{x1D719}}=0.40$ , $h/H\approx 0.95$ .

Figure 7. Shape of free surface. Simulations at $Ca=1$ .

Inside the channel, particles are concentrated near the symmetry plane and the wall layer is depleted of particles, which leads to a very intense viscosity gradient in the flow. Therefore, the velocity profile of the suspension flow at the channel exit is blunted when compared to the Newtonian liquid. As the suspension concentration increases, the viscosity near the centreline increases and thereby the centreline velocity decreases; in turn, the viscosity near the wall decreases, leading to a higher velocity gradient in this region. This result is shown in figure 8(a). Figure 8(b) shows the fluid velocity along the symmetry plane, and figure 9 displays the horizontal velocity field in the entire flow domain. For the Newtonian case, the centreline velocity inside the channel is constant and equal to $u/V=1.5$ , and experiences a sudden decrease to $u/V\approx 0.88$ as the free surface approaches. Since $u/V<1$ under the free surface, the flow decelerates and we observe a die expansion with respect to the channel gap in order to satisfy mass conservation. When particle migration is considered in the suspension flow, the particles are more concentrated near the flow centreline, such that the high viscosity in this region leads to a decrease in the centreline velocity. At $\bar{\unicode[STIX]{x1D719}}=0.10$ , the centreline velocity decreases from $u/V\approx 1.49$ inside the channel to $u/V\approx 0.91$ along the free surface. So, there is also a die expansion, but with an extrudate height smaller than the Newtonian case. At $\bar{\unicode[STIX]{x1D719}}=0.40$ , the centreline concentration is so high that the velocity falls from $u/V\approx 1.37$ inside the channel to $u/V\approx 1.05$ in the free surface. Since $u/V>1$ under the free surface, the extrudate velocity is higher than the flow mean velocity and a die contraction appears in order to satisfy mass balance. There are a few experimental reports on extrudate contraction in the flow of concentrated suspensions. Aral & Kalyon (Reference Aral and Kalyon1997) presented evidence on a variation in extrudate shape in the flow of dense suspensions of non-colloidal spheres. Nicolas (Reference Nicolas2002) observed a die contraction in gravity-driven jets of concentrated suspensions. A qualitatively similar result was reported by Furbank & Morris (Reference Furbank and Morris2004) on the formation of drops of particulate suspensions. However, in the works of Nicolas (Reference Nicolas2002) and Furbank & Morris (Reference Furbank and Morris2004), both inertial and gravitational effects were not negligible, so that we cannot be sure whether the die contraction they observed was caused by the same mechanism as discussed here. Extrudate contraction was also reported by Kharchenko et al. (Reference Kharchenko, Douglas, Obrzut, Grulke and Migler2004) in the flow of concentrated suspensions of carbon nanotubes, which was assumed to be a consequence of particle alignment in the flow direction.

Figure 8. Comparison of the velocity field in the flow of suspensions and a Newtonian liquid: (a) velocity profile at the channel exit ( $x_{1}/H=0$ plane); (b) centreline velocity at the free surface ( $x_{2}/H=0$ plane). Simulations at $Ca=1$ .

Figure 9. Horizontal velocity field in the entire flow domain for: (a) Newtonian liquid; (b) $\bar{\unicode[STIX]{x1D719}}=0.10$ ; (c) $\bar{\unicode[STIX]{x1D719}}=0.40$ . Simulations at $Ca=1$ .

The variations in suspension viscosity and flow velocity caused by particle migration imply a strong change in the suspension stress field near the channel exit and the entrance of the free-surface region. Figure 10 compares the stress field in this region for a Newtonian fluid, a dilute suspension at $\bar{\unicode[STIX]{x1D719}}=0.10$ and a concentrated suspension at $\bar{\unicode[STIX]{x1D719}}=0.40$ . In all cases, the fluids present the same average bulk viscosity. In these plots, the stress components are evaluated in units of $\unicode[STIX]{x1D70F}_{c}=\bar{\unicode[STIX]{x1D702}}V/H$ , the characteristic shear stress of the flow. The normal stress $\unicode[STIX]{x1D6F4}_{22}$ along the $x_{2}=0$ plane rises as the suspension bulk concentration increases. The resulting vertical net force along the plane is compressive for the Newtonian fluid. Therefore, the vertical component of the surface tension force acting at the contact line needs to be negative in order to satisfy momentum balance; the free-surface height is larger than the channel gap. On the other hand, the resulting vertical net force along the plane pulls the liquid downwards for the concentrated particle suspension at $\bar{\unicode[STIX]{x1D719}}=0.40$ . The vertical component of the surface tension force at the contact line needs to be positive, resulting in an extrudate contraction. This stress-based argument explains the change in the free-surface curvature from a Newtonian liquid and dilute suspensions, where a die expansion is observed, to concentrated suspensions, where a die contraction appears.

Figure 10. Stress field near the entrance of the free surface. (a–c) For $\unicode[STIX]{x1D6F4}_{22}/\unicode[STIX]{x1D70F}_{c}$ : (a) Newtonian liquid; (b) $\bar{\unicode[STIX]{x1D719}}=0.10$ ; (c) $\bar{\unicode[STIX]{x1D719}}=0.40$ . (d–f) For $\unicode[STIX]{x1D6F4}_{12}/\unicode[STIX]{x1D70F}_{c}$ : (d) Newtonian liquid; (e) $\bar{\unicode[STIX]{x1D719}}=0.10$ ; (f) $\bar{\unicode[STIX]{x1D719}}=0.40$ . Simulations at $Ca=1$ .

In the absence of inertia and gravity, the shape of the free surface is dictated by a balance between viscous and interfacial forces near the air–liquid interface, which is represented by the flow capillary number. Figure 11 shows the effect of the capillary number on the free jet final height. In all cases, the flow rate was kept fixed and the capillary number was varied only by changing the liquid surface tension. At very low capillary number, the flow is dominated by surface-tension forces, which act in the sense to avoid the curvature of the free surface. As a consequence, the extrudate final height is almost equal to the channel gap, i.e.  $h/H\approx 1$ , without a considerable die expansion or contraction. At $Ca=10^{-2}$ , $h/H\approx 1.005$ for the Newtonian liquid, $h/H\approx 1.002$ for suspensions at $\bar{\unicode[STIX]{x1D719}}=0.10$ and $h/H\approx 0.995$ for suspensions at $\bar{\unicode[STIX]{x1D719}}=0.40$ . As the capillary number increases and viscous effects become dominant, the free-surface curvature increases in order to balance the viscous force at the air–liquid interface. Therefore, the extrudate final height and the die expansion increase in the case of a Newtonian liquid and suspensions at low concentrations. In turn, the free jet final height decreases and the die contraction becomes more pronounced for suspensions at high concentrations. At the limit of very high capillary number, the relation $h/H$ tends to a constant asymptotic value. At $Ca=10^{2}$ , $h/H\approx 1.189$ for the Newtonian liquid, $h/H\approx 1.149$ for suspensions at $\bar{\unicode[STIX]{x1D719}}=0.10$ and $h/H\approx 0.936$ for suspensions at $\bar{\unicode[STIX]{x1D719}}=0.40$ . Although the capillary number has a strong effect on the extrudate shape, we verified that it has a negligible effect on the final particle distribution at the free-surface outflow plane when the $x_{2}$ coordinate is normalized by the corresponding free jet final height.

Figure 11. Effect of the capillary number on the free jet height.

5 Final remarks

Particle migration in planar extrudate flows was studied. The suspension was described as a Newtonian liquid with a concentration-dependent viscosity, and shear-induced particle migration was modelled according to the diffusive flux model. The set of fully coupled, nonlinear differential equations was solved by the DEVSS-TG/SUPG finite element method together with the elliptic mesh generation method to capture the position of the free surface.

The results show that particle concentration distribution at the free-surface outflow is highly non-uniform, presenting higher concentration near the flow symmetry plane and lower concentration near the air–liquid interface. Particle migration affects the local viscosity and flow velocity, and thereby changes the suspension stress field near the channel exit. For a Newtonian liquid and dilute suspensions, die swell was observed. On the other hand, at high concentrations, a die contraction was observed. Here, this change in the free-surface curvature is related to the stress field variation caused by particle migration. The curvature of the free surface increases with the capillary number in order to balance the viscous force near the air–liquid interface. Therefore, for a Newtonian liquid and dilute suspensions, extrudate height increases with capillary number, whereas for concentrated suspensions, the extrudate height falls as capillary number increases. These results may explain the die-swell suppression observed in experiments of extrusion of particle-filled suspensions.