A Primer on CFD-DEM for Polymer-Filled Suspensions

Abstract

This work reports on an evaluation of the computational fluid dynamics–discrete element method (CFD-DEM) numerical approach to study the behavior of polymer-filled suspensions in a parallel-plate rheometer. For this purpose, an open-source CFD-DEM solver is used to model the behavior of such suspensions considering different particle volume fractions and different types of fluid rheology. We first validate the numerical approach for the single-phase flow of the continuum phase (fluid phase) by comparing the fluid’s azimuthal velocity and shear stress components obtained from the open-source solver against the analytical expressions given in cylindrical coordinates. In addition, we compare the numerical torque given by the numerical procedure with analytical expressions obtained for Newtonian and power law fluids. For both cases, there is a remarkable agreement between the numerical and analytical results. Subsequently, we investigated the effects of the particle volume fraction on the rheology of the suspension. The numerical results agree well with the experimentally measured ones and show a yield stress phenomenon with the increase of the particle volume fraction.

1. Introduction

Filled suspensions are commonly used in various industries, such as the chemical industry, mining, food processing, and cosmetics. The properties of the particles suspended in the base fluid (e.g., size, shape, and surface properties) are important for the development of new and better materials/products. For example, numerous polysaccharides and proteins can combine to generate colloidal particles that can create a viscous solution [1], an elastic gel, or a liquid crystalline phase (which are frequently utilized as structuring or thickening ingredients in meals because of this). There is also particle erosion behavior in viscoelastic surfactant abrasive slurry jetting [2], where viscoelastic surfactant fluids are used due to the advantages of low friction, environmental friendliness, and microstructures recoverability after mechanical degradation. Other applications could be cited here, but the focus should be on the viscous and elastic nature of the fluids. This system’s rheological behavior is extremely complicated and strongly non-linear, because the system’s microstructure changes as a result of the applied deformation. The most simple model that takes into account such changes is the power law model. Therefore, in this work, we aim at understanding the behavior of such a complex system, where viscosity changes with the rate of deformation. This understanding is crucial in order to study more complex flows where elasticity plays a fundamental role. For low and weak deformations/flows, the model can be used to study viscoelastic materials; for high and strong deformations/flows, a more complex viscoelastic model should be used instead. Before we proceed to such complex flows, it is imperative to study the controlled deformation of particle flows modeled by a power law fluid.

Therefore, to understand the behaviors of the materials under different deformations and to identify the influences of the particles on the final product, rheological characterization is often performed.

This characterization is done with the help of rheometers (see Figure 1), which allow the fluid properties to be measured under shear stress or shear rate conditions (considering controlled flows).

Figure 1 shows a schematic of a rheometer consisting of an upper plate ($z=H$), a lower plate ($z=0$), and a gap between the plates, which is filled with a sample of the material to be tested. Both plates have a radius of r. The lower plate is fixed and the upper plate can rotate with an angular velocity $\Omega$, resulting in a deformation of the material sample ($\dot{\gamma }=r\Omega /H$, with $\dot{\gamma }$ being the shear rate). This deformation in a controlled environment allows the properties of the material to be studied based on the material’s response to the imposed deformation. This characterization is commonly supported by numerical methods, which allow simulations of filled suspensions to be carried out. In this way, it is possible to determine the behaviors of the suspensions at any time and any point in the domain.

Regarding the vast scientific literature about filled suspensions, early research studies focused on the relationship between relative viscosity and particle volume fraction, such as those presented by Rutgers [3,4] and Thomas [5]. A review of the literature on experimental studies of the rheology of suspensions shows that only the apparent viscosity of the suspension ($\eta =\tau /\dot{\gamma }$, where $\tau$ is the shear stress) is presented. However, to capture the non-Newtonian rheology, such as shear-thinning and yield stress, this one-dimensional measure is not sufficient. Some theoretical and experimental studies have attributed shear-thinning to the dependence of viscosity on the Péclet number, particle Reynolds number, or Stokes number of the suspension [6,7]. Other works, such as [8] investigated the dependence of the shear-thinning degree on the particle volume fraction $\varphi$. The dependence of shear-thinning on $\varphi$ suggests that its origin lies in the micromechanics of interparticle interactions. However, there is a gap in the scientific literature regarding the dependence of the yield stress of the suspension on the particle volume fraction. For this reason, in this work we will develop a numerical approach to capture the influence of the particle volume fraction on the yield stress of a suspension flow.

To model the continuum and dispersed phases, one must use governing equations that describe the phase behavior on a particular time and length scale. Then a suitable coupling scheme must be used so that the influence of one phase on the other can be captured. Various coupling schemes can be used, each with its own advantages and disadvantages [9]. Particle-laden flows are widely described in the literature by two popular coupling schemes: the continuum-continuum approach at a macroscopic level represented by the so-called two-fluid model (TFM) [10], and the continuum-discrete approach at a multiscale level that results from the combination of the Computational Fluid Dynamics (CFD) and Discrete Element Method (DEM), the so-called CFD-DEM approach [11,12].

In this paper, we evaluate the CFD-DEM method by simulating suspensions in a parallel-plate rheometer. We consider the CFD-DEM solver [13,14] using OpenFOAM for CFD and LIGGGHTS for DEM. These general approaches aim to significantly extend and generalize the capabilities of numerical investigation of the behavior of particle-filled systems. The aim is to show that the numerical results agree either qualitatively or quantitatively (depending on the availability of experimental data for comparison) with the experimentally measured results [15,16].

The interested reader should consult references [17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33]. For a brief explanation on the evolution of CFD-DEM, please consult [34].

The remainder of this work is structured as follows. In Section 2, we describe the governing equations that model the incompressible single-phase flow of Newtonian and power law fluids inside a parallel-plate rheometer. In addition, the numerical implementation and validation of the single-phase algorithm are performed for torque calculations. Section 3 is devoted to particle-laden flows. We start by presenting the incompressible volume-averaged Navier–Stokes equations and we then explain the CFD-DEM numerical procedure. In Section 4, we assess the CFD-DEM numerical model for simulations of suspensions involving several particle volume fractions. Finally, in Section 5, we summarize the main conclusions of this work.

2. Single-Phase Fluid Flows

2.1. Governing Equations

First, we will consider that an incompressible fluid is confined between two circular plates, as shown in the schematic diagram of Figure 1. To take advantage of the symmetry of this flow, we can write the governing equations for the motion of a laminar, incompressible, and isothermal fluid in cylindrical coordinates. The equation describing the conservation of mass is as follows,

$$\begin{array}{c}\hfill \frac{1}{r}\frac{\partial \left(r{v}_r\right)}{\partial r}+\frac{1}{r}\frac{\partial v_{\theta }}{\partial \theta }+\frac{\partial v_z}{\partial z}=0,\end{array}$$

(1)

and the linear momentum equations are given by

$$\begin{array}{cc}\hfill & \rho \left(\frac{\partial v_r}{\partial t}+v_r\frac{\partial v_r}{\partial r}+\frac{v_{\theta }}{r}\frac{\partial v_r}{\partial \theta }-\frac{v_{\theta }^2}{r}+v_z\frac{\partial v_r}{\partial z}\right)=-\frac{\partial p}{\partial r}+\left(\frac{1}{r}\frac{\partial \left(r{\tau }_{rr}\right)}{\partial r}+\frac{1}{r}\frac{\partial {\tau }_{\theta r}}{\partial \theta }-\frac{{\tau }_{\theta \theta }}{r}+\frac{\partial {\tau }_{rz}}{\partial z}\right)+\rho g_r,\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill & \rho \left(\frac{\partial v_{\theta }}{\partial t}+v_r\frac{\partial v_{\theta }}{\partial r}+\frac{v_{\theta }}{r}\frac{\partial v_{\theta }}{\partial \theta }+\frac{v_{\theta }{v}_r}{r}+v_z\frac{\partial v_{\theta }}{\partial z}\right)=-\frac{1}{r}\frac{\partial p}{\partial \theta }+\left(\frac{1}{r^2}\frac{\partial \left(r^2{\tau }_{r\theta }\right)}{\partial r}+\frac{1}{r}\frac{\partial {\tau }_{\theta \theta }}{\partial \theta }+\frac{\partial {\tau }_{\theta z}}{\partial z}\right)+\rho g_{\theta },\hfill \end{array}$$

(2b)

$$\begin{array}{cc}\hfill & \rho \left(\frac{\partial v_z}{\partial t}+v_r\frac{\partial v_z}{\partial r}+\frac{v_{\theta }}{r}\frac{\partial v_z}{\partial \theta }+v_z\frac{\partial v_z}{\partial z}\right)=-\frac{\partial p}{\partial z}+\left(\frac{1}{r}\frac{\partial \left(r{\tau }_{rz}\right)}{\partial r}+\frac{1}{r}\frac{\partial {\tau }_{\theta z}}{\partial \theta }+\frac{\partial {\tau }_{zz}}{\partial z}\right)+\rho g_z.\hfill \end{array}$$

(2c)

Here, t is time, $\underline{v}=(v_r,v_{\theta },v_z)$ are the three components of fluid velocity vector $\underline{v}$ in cylindrical coordinates $(r,\theta ,z)$, ${\tau }_{rr}$, ${\tau }_{\theta \theta }$, ${\tau }_{zz}$, are the normal stress components, ${\tau }_{r\theta }$, ${\tau }_{rz}$, ${\tau }_{z\theta }$ are the shear stress components of the stress tensor $\underline{\underline{\tau }}$, $\rho$ is the fluid density and $\underline{g}=(g_r,g_{\theta },g_z)$ are the three components of the gravity vector $\underline{g}$.

For a Newtonian fluid, the stress tensor is given by,

$$\underline{\underline{\tau }}=2\eta \underline{\underline{D}}=\eta \left(\nabla \underline{v}+{(\nabla \underline{v})}^T\right),$$

(3)

with $\eta$ being the fluid dynamic viscosity, $\underline{\underline{D}}$ is the rate of the deformation tensor and $\nabla \underline{v}$ is the gradient of velocity (both given in cylindrical coordinates).

For the power law fluid, the stress tensor is given by,

$$\underline{\underline{\tau }}=2\eta \left(\right|\underline{\underline{D}}\left|\right)\underline{\underline{D}},$$

(4)

where the viscosity is now a function of the magnitude of the rate of the deformation tensor ($|\underline{\underline{D}}|$),

$$\eta \left(\right|\underline{\underline{D}}\left|\right)=k|\underline{\underline{D}}{|}^{n-1},$$

(5)

with n being the power law constant and k being the flow consistency index.

2.2. Analytical Torque Profiles

Consider that the gap between the plates is H and the top plate rotates with a constant angular velocity $\Omega$. This flow can be considered an axisymmetric torsional flow because only the azimuthal velocity component is nonzero $\left(\mathrm{i}.\mathrm{e}.,v_{\theta }\ne 0,v_r=0,v_z=0\right)$. Thus, for this geometry, the velocity field satisfying Equations (1) and (2a)–(2c) is given by,

$$\begin{array}{c}\hfill \underline{v}={\left(\begin{array}{c}0\\ \frac{r\Omega z}{H}\\ 0\end{array}\right)}_{r\theta z}.\end{array}$$

(6)

We then calculate the total torque $\left(\underline{\mathcal{T}}\right)$ on the top fluid surface $\left(S\right)$ in the parallel-plate rheometer as [35],

$$\begin{array}{c}\hfill \underline{\mathcal{T}}=\underset{S}{\int \int }{[\underline{R}\times (\underline{\widehat{n}}·\underline{\underline{\tau }})]}_{\mathrm{at}\mathrm{surface}}dS.\end{array}$$

(7)

Note that to compute the torque with the numerical modeling data, we choose an infinitesimally small area $dS=rd\theta dr$, which is located at a distance of r from the axis of the rotation. The lever arm vector $\underline{R}$ is a vector from the axis of rotation to the axis where the torque should be computed. For this infinitesimal area, the lever arm vector is $\underline{R}=r\underline{{\widehat{e}}_r}$. The surface of the fluid in contact with the top plate has a unit normal vector $\underline{\widehat{n}}=\underline{{\widehat{e}}_z}$. For comparison purposes, the numerical torque obtained through finite volume simulations will be compared with the analytical torque expression on the top fluid surface in the parallel-plate rheometer for a Newtonian fluid,

$$\begin{array}{c}\hfill \underline{\mathcal{T}}=\frac{\pi \Omega \eta r^4}{2H}\underline{{\widehat{e}}_z}={\left(\begin{array}{c}0\\ 0\\ \frac{\pi \Omega \eta r^4}{2H}\end{array}\right)}_{xyz}={\left(\begin{array}{c}0\\ 0\\ \frac{\pi \Omega \eta r^4}{2H}\end{array}\right)}_{r\theta z},\end{array}$$

(8)

and for a power law fluid,

$$\begin{array}{c}\hfill \underline{\mathcal{T}}=2\pi k{\left(\frac{\Omega }{H}\right)}^n\left(\frac{r^{n+3}}{n+3}\right)\underline{{\widehat{e}}_z}={\left(\begin{array}{c}0\\ 0\\ 2\pi k{\left(\frac{\Omega }{H}\right)}^n\left(\frac{r^{n+3}}{n+3}\right)\end{array}\right)}_{xyz}={\left(\begin{array}{c}0\\ 0\\ 2\pi k{\left(\frac{\Omega }{H}\right)}^n\left(\frac{r^{n+3}}{n+3}\right)\end{array}\right)}_{r\theta z}.\end{array}$$

(9)

2.3. Numerical Procedure and Meshes

The governing equations for a single fluid phase are solved using the finite-volume method together with the SIMPLE algorithm (a procedure for the calculation of pressure and velocity, by guessing and correcting pressure and velocity fields) [36]. For the case of power law fluid, we have to use an iterative procedure to update the viscosity and stress components.

The velocity gradient is computed using a second-order least-squares method, and the diffusive term in the momentum balance is discretized using second-order linear interpolation. The advective terms in the momentum and constitutive equations are discretized using the second-order linear-upwind scheme, which follows a component-wise and deferred correction approach to improve numerical stability. The time derivatives are discretized using the bounded implicit second-order Crank–Nicolson scheme. Therefore, the overall accuracy of the numerical method is equal to two for both spatial and temporal discretizations.

The numerical torque calculations were assessed using the forces postProcessing utility [37] available in OpenFOAM by comparison with the analytical profiles. This tool computes the forces and torques by integrating the stress tensor effect in a given list of patches (sets of boundary faces). For moment/torque calculations, the center and axis of rotation need to be defined.

As shown in Figure 2, a O-block mesh was used to discretize the geometrical domain represented in Figure 1. Three different mesh refinement levels were created using the blockMesh application, being defined as M1: 15 × 30 × 4, M2: 30 × 60 × 8 and M3: 60 × 120 × 16, where M#: $R\times \theta \times H$ represents the spatial discretizations on the radial, azimuthal and longitudinal directions, respectively. The accuracies of the numerical torque $\underline{\mathcal{T}}$ calculations were accessed through the comparison with the analytical results given by Equations (8) or (9) via the application of Richardson’s extrapolation [38] to the limit.

For the numerical calculations the following parameters were used: the cylinder has a height of $H=0.5\mathrm{mm}$ and a radius of $r=15\mathrm{mm}$. The angular velocity was imposed on the top plane using the rotatingWallVelocity boundary condition available on the OpenFOAM CFD framework. Simulations were done for shear rates between 0 and $100{\mathrm{s}}^{-1}$. The Newtonian fluid has a dynamic viscosity of $\eta =0.001$ Pa.s and density of $\rho =1000$ kg/m³. For the power law fluid, the flow consistency index k adopted the same value of the Newtonian viscosity and the power law exponent n was varied between 0 and 1.

2.4. Single-Phase Fluid Model Assessment

In this section, we start by validating the methodology for the calculation of the fluid velocity, stress tensor components, and torque expressions derived in Section 2.2. In Figure 3, we show the fluid azimuthal velocity component $v_{\theta }$ and the shear stress component ${\tau }_{\theta z}$ at several radius locations. As expected from Equation (6), the azimuthal fluid velocity component profile is linear and the shear stress component is constant along the z-direction for each radial location.

Figure 4 shows the contour of the magnitude of the fluid velocity for a Newtonian fluid moving inside a parallel-plate rheometer. As expected, the magnitude of the fluid velocity $\underline{v}$ varies from zero at the bottom plate to a maximum value at the top plate. In addition, for a fixed height in the parallel-plate rheometer, the velocity profile is linear from the center to the edge of the plate.

The numerically calculated torque for a Newtonian fluid is shown in Figure 5 and is compared with the exact value given by Equation (8). The relative error, in percentage, between the numerical torque results and the values of the analytical solution is also shown in Figure 5. For all shear rates tested ($0\le \dot{\gamma }\le 100$), the relative errors are less than $0.15\%$. Note that the radius is fixed; thus, the dependence of the torque profile on the radius is not observed.

The torque results for a power law fluid were also numerically calculated and compared with the analytical solution from Equation (9) for $n=0.3,0.6$, and $0.9$ (see Figure 6). Again, the numerical results and the values of the analytical solution agree remarkably well. We no longer have a linear profile. Introducing a nonlinearity into the constitutive equation leads to a nonlinear torque profile that depends on the power of n. For higher values of n ($n\stackrel{}{\to }1$) the Newtonian profile is restored and we observe an almost linear profile. For low values of n, the non-Newtonian flow behavior is more pronounced and the fluid has more difficulty in increasing the torque for higher shear rates. This means that the stresses on the wall are lower due to the strong shear-thinning effect leading to lower viscosity.

We can, therefore, conclude that the numerical implementation for the single-phase fluid reproduces well the physics of both Newtonian and non-Newtonian flows.

3. Particle-Laden Flows

A coupled CFD-DEM approach was used for modeling particle-laden flows [13]. The so-called unresolved approach [13,39,40] was employed in this study because only the bulk behavior of the suspension is of interest. In this formulation, the particle size is smaller than the computational grid size. Thus, it is assumed that the particles do not completely fill a computational cell. Consequently, the particles are not resolved in the CFD simulation, as in the case of the immersed boundary method [41], but only their interaction with the fluid phase is considered in terms of momentum transfer. Hereafter, the equations for the dispersed phase (DEM) and continuum phase (CFD) are described.

3.1. Governing Equations for the Dispersed Phase (DEM)

To model the dispersed phase, the discrete element method (DEM) is employed. Following [42], the DEM is used to solve Newton’s second law equations for the motion of a particle i in the suspension

$$\begin{array}{c}\hfill m_i\frac{d{\underline{v}}_i^p}{dt}=\sum _{j=1}^{n_i^c}{\underline{f}}_{ij}^c+\sum _{k=1}^{n_i^{nc}}{\underline{f}}_{ik}^{nc}+{\underline{f}}_i^{pf}+{\underline{f}}_i^g,\end{array}$$

(10)

$$\begin{array}{c}\hfill I_i\frac{d{\underline{\omega }}_i^p}{dt}=\sum _{j=1}^{n_i^c}{\underline{M}}_{ij},\end{array}$$

(11)

where $m_i$ and $I_i$ are particle’s i’s mass and moment of inertia, respectively, ${\underline{v}}_i^p$ and ${\underline{\omega }}_i^p$ are particle i’s linear and angular velocities, respectively, ${\underline{f}}_{ij}^c$ and ${\underline{M}}_{ij}$ are the contact force and torque promoted by the j contact (caused by a particle or a wall) with particle i, respectively; $n_i^c$ is the total number of contacts experienced by particle i, ${\underline{f}}_{ik}^{nc}$ are the non-contact (long-range) forces promoted by particle k or other sources on particle i; $n_i^{nc}$ is the total number of non-contacts experienced by particle i, ${\underline{f}}_i^{pf}$ are the particle–fluid interaction forces, and ${\underline{F}}_i^g$ is the force due to the gravitational acceleration vector.

In this work, the contact force due to the interaction between two rigid particles, ${\underline{f}}_{ij}^c$, was computed by the Hertzian spring dashpot model [43]. The normal component contribution, ${\underline{f}}_{nij}^c$, resulting from the contact of particle j with particle i, is given by

$$\begin{array}{c}\hfill {\underline{f}}_{nij}^c=(-k_n|{\delta }_n|-{\eta }_n{\underline{v}}_{ij}^p·\underline{n})\underline{n},\end{array}$$

(12)

and, if the contact is between particle i and a wall, ${\underline{f}}_{niw}^c$, then it is calculated as

$$\begin{array}{c}\hfill {\underline{f}}_{niw}^c=(-k_{nw}|{\delta }_{nw}|-{\eta }_{nw}{\underline{v}}_{iw}^p·{\underline{n}}_w){\underline{n}}_w.\end{array}$$

(13)

where $\underline{n}$ is the unit vector pointing from the center of particle i to that of particle j and ${\underline{n}}_w$ is the unit vector of the surface normal pointing outwards. Here, we should know beforehand the stiffness and the dam** coefficients of particles in the normal directions, $k_n$ and ${\eta }_n$, respectively, and also the ones due to the contact of a particle with the walls, $k_{nw}$ and ${\eta }_{nw}$. In addition, for the spring-dashpot model, the particles are allowed to overlap an amount of ${\delta }_n$ and ${\delta }_{nw}$ normal displacements when in contact with another particle or with a wall, respectively. The magnitude of these displacements read as follows

$$\begin{array}{c}\hfill |{\delta }_n|=r_i+r_j-|{\underline{p}}_j-{\underline{p}}_i|,\end{array}$$

(14)

and

$$\begin{array}{c}\hfill |{\delta }_{nw}|=r_i-|{\underline{p}}_i-{\underline{p}}_c|,\end{array}$$

(15)

where $r_i$ and $r_j$ are the radii of particles i and j, respectively, ${\underline{p}}_i$, ${\underline{p}}_j$, and ${\underline{p}}_c$ are the vectors with the coordinates that give the locations of particle i, particle j, and the contact point on the wall, respectively. Subsequently, due to the particle–particle interaction, we need to define the velocity of particle i relative to particle j, ${\underline{v}}_{ij}^p$,

$$\begin{array}{c}\hfill {\underline{v}}_{ij}^p={\underline{v}}_i^p-{\underline{v}}_j^p,\end{array}$$

(16)

and, if the contact is between particle i and a wall,

$$\begin{array}{c}\hfill {\underline{v}}_{iw}^p={\underline{v}}_i^p-{\underline{v}}_w^p,\end{array}$$

(17)

where ${\underline{v}}_i^p$ and ${\underline{v}}_j^p$ are particle i and particle j linear velocities, respectively, and ${\underline{v}}_w^p$ is the slip velocity of the sphere-wall contact point.

On the other hand, the tangential component contribution, ${\underline{f}}_{tij}^c$, resulting from the contact of particle j with particle i, reads as follows

$$\begin{array}{c}\hfill {\underline{f}}_{tij}^c=-k_t{\delta }_t-{\eta }_t{\underline{v}}_{tij}^p,\end{array}$$

(18)

and, if the contact is between particle i and a wall,

$$\begin{array}{c}\hfill {\underline{f}}_{tiw}^c=-k_{tw}{\delta }_{tw}-{\eta }_{tw}{\underline{v}}_{tiw}^p.\end{array}$$

(19)

Again, we should know beforehand the stiffness and the dam** coefficients of particles in the tangential direction, $k_t$ and ${\eta }_t$, respectively, and also the ones due to the contact of a particle with the walls, $k_{tw}$ and ${\eta }_{tw}$. In addition, the tangential displacements, ${\delta }_t$ and ${\delta }_{tw}$, resultant from the particle overlap with another particle or with a wall, respectively, are computed accordingly to the expressions given in [44]. In addition, we need to define the tangential slip velocity of particle i relative to particle j, ${\underline{v}}_{tij}^p$,

$$\begin{array}{c}\hfill {\underline{v}}_{tij}^p={\underline{v}}_{ij}^p-({\underline{v}}_{ij}^p·\underline{n})\underline{n}+(r_i{\underline{\omega }}_i^p+r_j{\underline{\omega }}_j^p)\times \underline{n},\end{array}$$

(20)

and if the contact is between particle i and a wall,

$$\begin{array}{c}\hfill {\underline{v}}_{tiw}^p={\underline{v}}_{iw}^p-({\underline{v}}_{iw}^p·\underline{n}){\underline{n}}_w+r_i{\underline{\omega }}_i^p\times {\underline{n}}_w.\end{array}$$

(21)

Notice that if the following relation is satisfied

$$\begin{array}{c}\hfill |{\underline{f}}_{tij}^c|\ge {\mu }_f\left|{\underline{f}}_{nij}^c\right|,\end{array}$$

(22)

then the Coulomb-type friction law should be applied because particle i will slide over particle j or at the wall. In this case, the tangential force is given by

$$\begin{array}{c}\hfill {\underline{f}}_{tij}^c=-{\mu }_f\left|{\underline{f}}_{nij}^c\right|\frac{{\delta }_t}{|{\delta }_t|},\end{array}$$

(23)

where ${\mu }_f$ is the friction coefficient.

Notice also that, in this work, the particles contacts follow the Hertzian contact theory [45]. This way, the normal component of the contact force, ${\underline{f}}_{nij}^c$, given by Equation (12), needs to be reformulated as

$$\begin{array}{c}\hfill {\underline{f}}_{nij}^c=(-K_n|{\delta }_n{|}^{3/2}-{\eta }_n{\underline{v}}_{ij}^p·\underline{n})\underline{n},\end{array}$$

(24)

where the stiffness coefficient, $K_n$, is evaluated using the Hertzian contact theory, meaning that physical properties such as Young’s modulus and Poisson ratio need to be known [43]. For the calculation of the other physical parameters, such as dam** and friction coefficients, the details can also be found in [43]. The equations for the calculation of the particle–particle interaction forces and torques can also be found in detail in [46].

Regarding non-contact forces, such as van der Waals or electrostatic forces, they are neglected in this work since they are much smaller in magnitude than the hydrodynamics or contact forces for the particles considered herein.

The particle–fluid interaction force, ${\underline{f}}_i^{pf}$, is the sum of all types of force contributions acting on individual particles by the fluid, including drag force ${\underline{f}}_i^d$, pressure gradient force ${\underline{f}}_i^{\nabla p}$, viscous force ${\underline{f}}_i^{\nabla ·\underset{=}{\tau }}$, virtual mass force ${\underline{f}}_i^{vm}$, Basset force ${\underline{f}}_i^{Bass}$, and lift forces, such as the Saffman force ${\underline{f}}_i^{Saff}$ and Magnus force ${\underline{f}}_i^{Mag}$. This way, the total particle–fluid interaction force on an individual particle i is given by

$$\begin{array}{c}\hfill {\underline{f}}_i^{pf}={\underline{f}}_i^d+{\underline{f}}_i^{\nabla p}+{\underline{f}}_i^{\nabla ·\underset{=}{\tau }}+{\underline{f}}_i^{vm}+{\underline{f}}_i^{Bass}+{\underline{f}}_i^{Saff}+{\underline{f}}_i^{Mag}.\end{array}$$

(25)

The expressions for the particle–fluid forces used in this work will be given in the next section.

3.2. Governing Equations for the Continuum Phase (CFD)

In this work, set II (or form A in [10]) of the incompressible volume-averaged Navier–Stokes equation is considered for the fluid phase [47]. The volume-averaged continuity equation reads as

$$\begin{array}{c}\hfill \frac{\partial {ϵ}_f}{\partial t}+\nabla ·\left({ϵ}_f\underline{v}\right)=0,\end{array}$$

(26)

and the volume-averaged linear momentum equations are given by

$$\begin{array}{c}\hfill \frac{\partial \left({ϵ}_f\underline{v}\right)}{\partial t}+\nabla ·\left({ϵ}_f\underline{v}\underline{v}\right)=-{ϵ}_f\nabla P-{\underline{F}}^{pf}+{ϵ}_f\nabla ·\underset{=}{\tau }+{ϵ}_f\underline{g},\end{array}$$

(27)

where ${ϵ}_f$ is the fluid volume fraction (also denoted as the void fraction), P is the modified pressure $(p/\rho )$ and the fluid-phase viscous stress tensor $\underset{=}{\tau }$ is given by:

$$\begin{array}{c}\hfill \underset{=}{\tau }=\frac{\eta }{\rho }\left[(\nabla \underline{v})+{(\nabla \underline{v})}^T-\frac{2}{3}\left(\nabla ·\underline{v}\right){\delta }_K\right],\end{array}$$

(28)

where ${\delta }_K$ is the Kronecker delta. Recall that for a power law fluid the viscosity is a function of the magnitude of the rate of deformation tensor $\eta \left(\right|\underline{\underline{D}}\left|\right)=k|\underline{\underline{D}}{|}^{n-1}$.

During the CFD-DEM calculations several particle–fluid momentum exchange correlations can be applied [47]. In this work, the momentum exchange term from the particles to the fluid is enforced via the source term ${\underline{F}}^{pf}$ in the volume-averaged momentum equations (see the second term on the right-hand side of Equation (27)). Again, following the formulation of set II in [47], the momentum exchange term ${\underline{F}}^{pf}$ reads as follows

$$\begin{array}{c}\hfill {\underline{F}}^{pf}=\frac{1}{\Delta V_i}\sum _{i=1}^{n_p}\left({\underline{f}}_i^d+{\underline{f}}_i^{vm}+{\underline{f}}_i^{Bass}+{\underline{f}}_i^{Saff}+{\underline{f}}_i^{Mag}\right),\end{array}$$

(29)

where $\Delta V_i$ is the volume of the fluid cell where particle i is located and $n_p$ is the number of particles in that cell.

Taking into account Equations (25) and (29) we would need to find expressions for each force term. However, depending on the physics of the problem under consideration, it can be possible to neglect some of the force terms. In this work, we only consider the drag, pressure, and viscous forces.

The drag force ${\underline{f}}_i^d$ on a particle in the fluid is calculated here by the Di Felice [48] correlation, which is valid for a wide variety of both fixed-bed and suspended-particle systems, being given by

$$\begin{array}{c}\hfill {\underline{f}}_i^d=\frac{1}{8}{C}_d\rho \pi d_p^2\left(\underline{v}-{\underline{v}}_i^p\right)|\underline{v}-{\underline{v}}_i^p|{\epsilon }_f^{-\chi },\end{array}$$

(30)

where $d_p$ is the particle diameter, ${\underline{v}}_i^p$ is the translational velocity of particle i and $C_d$ is the drag coefficient expressed as

$$\begin{array}{c}\hfill C_d={\left(0.63+\frac{4.8}{\sqrt{R{e}_p}}\right)}^2,\end{array}$$

(31)

with the particle Reynolds number, $R{e}_p$, defined by

$$\begin{array}{c}\hfill R{e}_p=\frac{{\epsilon }_f\rho d_p|\underline{v}-{\underline{v}}_i^p|}{\eta },\end{array}$$

(32)

and the empirical expression for the exponent $\chi$ is given by

$$\begin{array}{c}\hfill \chi =3.7-0.65\exp \left[-\frac{{\left(1.5-\log \left(R{e}_p\right)\right)}^2}{2}\right].\end{array}$$

(33)

The pressure and viscous forces are used because, with the unresolved CFD-DEM approach, the particles are not discretized explicitly in the CFD part [49,50]. The pressure gradient force, ${\underline{f}}_i^{\nabla p}$, which results from the difference in pressure across the fluid and particle’s surface, reads as follows [47]

$$\begin{array}{c}\hfill {\underline{f}}_i^{\nabla p}=-V_i\nabla p,\end{array}$$

(34)

where $V_i$ is the volume of particle i. Another dominant fluid–particle momentum exchange force is the viscous force, ${\underline{f}}_i^{\nabla ·\underset{=}{\tau }}$, due to the fluid shear stress or deviatoric stress tensor [47]

$$\begin{array}{c}\hfill {\underline{f}}_i^{\nabla ·\underset{=}{\tau }}=-V_i(\nabla ·\underset{=}{\tau }).\end{array}$$

(35)

Virtual mass, Basset force, and lift forces (Saffman and Magnus) are neglected due to the very low particle relaxation time $\left(t_p=d_p^2{\rho }_p/18\eta \right)$, and small relative velocity between the viscous fluid and the particles [50].

4. Rheology of Filled Systems

After validating the methodology used for torque calculation with a single-phase flow of the Newtonian and power law fluids, a CFD-DEM numerical model for multiphase flows of suspensions was assessed. The main coupled model parameters used were: spherical particles of radius $r=100$ ηm and density ${\rho }_p=1000$ kg/m³, a Newtonian fluid with the same rheological parameters of Section 2.3, DEM timestep = $1\times {10}^{-7}$ s, CFD timestep = $1\times {10}^{-5}$ s, coupling interval = 100, form A for the coupled flow [10,47], Di Felicie’s drag model, pressure gradient force model, and viscous force model. The Young’s modulus used was equal to $3\times {10}^9$ Pa and the Poisson ratio was equal to 0.33. The dam** and friction coefficient values were 0.9 and 0.3, respectively.

It should be remarked that we obtained the steady-state results by considering both the evolution of the velocity and pressure residuals, and, the evolution of the torque.

Figure 7 shows the comparison between numerical torque results for suspensions without particles and those with 10% and 30% of particle volume fractions for a Newtonian fluid. The numerical results shown in Figure 7 were fitted by a first-order polynomial. We can see that the ordinate of the intersection point of those polynomials with the y-axis increases with the increase of the particles volume fraction, which indicates that yield stress phenomena occur when the fluid is in the presence of particles, in accordance with the results published in the literature [8,15]. When increasing the shear rate, the fluid seems to show solid-like behavior at first (at really low shear rates), i.e., it does not deform any further, but at higher shear rates the fluid deforms. This could be due to the fact that the presence of solid particles gives the suspension a more solid-like behavior at low deformations. At higher deformations, the liquid-like behavior prevails and the suspension flows.

Figure 7 also shows that the torque increases with the increase in the number of solid particles in suspension. This result was to be expected as the presence of particles leads to more resistance within the suspension. This resistance is further enhanced by the interaction between the particles, the drag of the particles, and gravity. Note also that at higher shear rates, the influences of different particle percentages on the torque are more pronounced, as expected.

Figure 8 shows the particle distribution of a suspension inside a parallel-plate rheometer with a gap of 0.5 mm and 10% particle volume fraction. In Figure 8a we see the view from below (fixed plate) and in Figure 8b the view from above (moving plate).

It is interesting to see the transient dynamics of the suspension flow. For $t=1$, the number of particles seen from above is small, but for $t=4$ the number of particles is much higher, giving the impression that the suspended particles are moving toward the upper wall. However, looking at the view from below, the number of particles appears to be more or less the same over time. This means that the particles that migrate to the upper wall are far away from the lower and upper plates. The reason for this behavior can be attributed to the fact that the mixed interaction between particles, particles and walls, the rotation of the upper plate, and gravity lead to a pulling force that can attract the middle particles (in the z-direction).

Note also the particle formation in the view from below (Figure 8a). Due to the particle–particle interaction and the resistance that the fluid exerts on the particles, a radial distribution of particles is observed for $t=1$. As time increases, the particles become more random and dispersed, as expected. This result is not seen in the view from above and could therefore also be related to gravity.

5. Conclusions

In this work, a CFD-DEM methodology was used to simulate suspensions in a parallel-plate rheometer. The CFD model was validated by comparisons between numerical and analytical torque calculations for both Newtonian and power law fluids.

• The numerically calculated torque for a Newtonian fluid and a fixed radius shows a linear profile, with the torque increasing with the increasing shear rate;
• The numerically calculated torque for a power law fluid and a fixed radius shows a non-linear profile, with the torque also increasing with the increasing shear rate. In this case, the torque depends on the shear rate to the power of n (power law constant), i.e., ${\left(\frac{\Omega }{H}\right)}^n$. For highly thinning fluids, the torque shows a slow increase;
• The torque depends on the radius (r) of the plates. It depends on $r^4$ for the Newtonian fluid and on $r^{n+3}/(n+3)$ for the power law fluid. This means that the torque suffers drastic changes when the fluid goes from shear thinning to shear thickening (or vice-versa).

Subsequently, a CFD-DEM model was used to simulate particle suspensions for two different volume fractions, $10\%$ and $30\%$.

• The numerical results agree with those measured experimentally and show a yield stress phenomenon with the increase of the volume fraction of the particles;
• The torque increases with the number of solid particles in the suspension because the presence of particles leads to more flow resistance. This resistance is further increased by the interaction between the particles, the drag of the particles, and gravity;
• At higher shear rates, the influences of different particle percentages on the torque are more pronounced. For a shear rate of 100 $s^{-1}$, we see that the increase in torque (compared to the case of the fluid with no particles) is 5% for $\varphi =10\%$ and 35% for $\varphi =30\%$.

Note also that the unsteady dynamics of the suspension flow show that for $t=1$ s the number of particles seen from above is small, but for $t=4$ s the number of particles is much higher, giving the impression that the suspended particles are moving towards the upper plate. A radial distribution of particles is also observed for $t=1$ s near the bottom plate. With increasing time, the particles become more random and dispersed, as expected.

To further understand this phenomenon, in the future, we will consider flows with and without viscoelasticity. The methodology of CFD-DEM will also be used to study and optimize complex industrial particle fluid processes.