A numerical method for simulating concentrated rigid particle suspensions in an elongational flow using a fixed grid

Abstract

In this work a new numerical method for concentrated inertialess rigid particle suspensions in a planar elongational flow using a fixed mesh is presented. The main concept is to randomly relocate a particle on an inflow section of the domain when it crosses the outflow boundaries. A three-layer domain is considered in order to: (i) develop a small computational domain as the representative sample of the whole suspension, (ii) impose the elongational flow boundary conditions far from the particles, (iii) achieve a steady state (in a statistical meaning). Our scheme uses a time-independent fixed grid avoiding the difficulties involved in deforming meshes and remeshing of the domain. In this way, computations can proceed indefinitely and micro-structural fluctuations around a steady state can be studied.

A fictitious domain is implemented in order to easily manage the rigid-body motion. The particles are described by their boundaries only (rigid-ring description) and the rigid-body motion is imposed through Lagrange multipliers. The bulk properties are recovered by using an averaging procedure where the traction forces on the particle surface are recovered by the Lagrange multipliers.

The scheme has been combined with a standard velocity–pressure finite element formulation and 2D simulations of a large number (150 and 225) of particles in a Newtonian medium are performed. Local as well as bulk properties are evaluated and discussed. The results show very good agreement with dilute theory as well as with other numerical simulations in the literature for higher concentrations.

Our formulation is well suited for viscoelastic suspensions and can be easily extended to 3D simulations.

Introduction

In the last decades, direct numerical simulations (DNSs) techniques have been developed in order to predict and understand the complex flow of particle-filled fluids. The motion of the fluid is governed by the (Navier–)Stokes equations and the motion of the particles by the linear and angular momentum equations of rigid-body dynamics. The coupling of the fluid and the particles is achieved through the no-slip condition on the particle boundaries and the hydrodynamic forces and torques on the particles. The hydrodynamic forces and torques are, of course, those arising from the computed motion of the fluid, and therefore are not known in advance. It has to be pointed out that no approximation for these forces and torques is made. So, in DNS methods, hydrodynamic interactions are not modeled but computed.

An increasing interest in rigid particle suspensions can be observed. Typically, in all systems of practical interest, the concentration of particles is high so they are non-dilute or concentrated. In other words, a many-particle system should be considered and the hydrodynamic interactions play a crucial role by affecting the local flow fields, bulk properties and the final behavior of the material.

In order to manage the problem computationally, we need to develop a suited simulation scheme that uses the smallest domain that still has the same average properties as the whole suspension. Hence, by solving the flow problem in this domain, we are able to predict the average micro-structure and the bulk properties of the suspensions, with reduced CPU time and memory.

This idea has been used by Hwang et al. [1], [2] where the authors combine Lees–Edwards boundary conditions, i.e. a sliding bi-periodic domain, with a standard velocity–pressure finite element formulation for a Newtonian suspension as well as a DEVSS/DG (discrete elastic viscous split stress/discontinuous Galerkin) scheme for viscoelastic suspensions in simple shear flow. According to this scheme, each frame slides relatively to one another by an amount determined by a given shear rate. So, a frame can be considered as a sample of the whole suspension and transforms the many-particle suspension into a single unit cell.

Recently, the bi-periodic frame concept has been extended to planar extensional flow [3]. However, in order to deal with such a flow, a deformation in time of the bi-periodic frame is proposed. As a consequence, after a certain time, the frames cannot be deformed anymore since the smallest characteristic length of the frame is comparable with the characteristic dimension of the particles. Hence, it is difficult to achieve a steady state for this imposed flow field. This is especially true for a viscoelastic fluid at high Weissenberg number, where large strains are needed before a steady state is obtained. Finally, in the scheme described in [3] remeshing of the domain is also done in order to keep the aspect ratio of the elements close to one.

In this work, we propose a new simulation scheme that circumvents these problems. The main concept is to relocate a particle on the inflow boundary of the domain when it crosses the outflow sections. So, no periodic boundary condition is imposed. In particular, the computational domain is divided into three concentric square regions: in the internal one the particles move, the micro-structural and bulk properties are evaluated in this region. So, this region can be considered as a sample of the whole suspension. In the intermediate region the particles can move as well and, when they cross the outflow boundaries of that region, they are relocated randomly on one of the two inflow sections of the same region. Finally, the outer region only contains fluid since particles cannot enter. The elongational flow boundary conditions are imposed on the external boundaries of the outer region: so the particles feel the presence of the elongational flow boundary conditions only as an imposed ‘far field’.

According to this scheme, no deformation of the domain occurs and a time-independent fixed grid can be used (and no remeshing of the domain is needed). Furthermore, an average steady state can be achieved: we do not need to stop the simulation since the domain dimensions do not change. Finally, this scheme is suited for the simulation of viscoelastic suspensions. Indeed, after the relocation of the particles in the intermediate region, the stress has time to develop before particles enter the internal region where the properties are calculated.

In this work, we consider a concentrated suspension of rigid, non-Brownian disks in a planar elongational flow, where the particle and fluid inertia can be neglected. Indeed, a vast literature is based on the inertialess assumption (see for example the method of Stokesian dynamics [4] for concentrated suspensions in Newtonian fluids). Our final goal is to compare the results for a Newtonian suspension with the ones where the suspending fluid is a melt and can be represented by a viscoelastic model having a very high viscosity. So, the Newtonian suspensions under investigation are characterized by a high viscosity as well and therefore the inertialess assumption is appropriate here.

The analysis is carried out for a Newtonian medium. The particle–fluid interactions are taken into account by implementing a Lagrange multiplier/fictitious domain method (LM/FDM) [5], [6]. The force-free, torque-free rigid body motion of the particles is described by a rigid-ring problem [1], [2]. So, a fixed mesh is used for the computation and the particles are described by their boundaries only, through collocation points. This description is possible because inertia is neglected. Finally, the rigid-body motion constraints are imposed through Lagrange multipliers that can be identified as traction forces on the particle surfaces (with a correction due to the fluid stress inside the object).

Another difference with the works of Hwang et al. [1], [2], [3] is that with our scheme a particle is not splitted into parts since it never crosses the boundary of the whole domain. However, since a particle can cross the sample internal region, a slight modification of the bulk stress formula is required.

We limit the simulations to two dimensions in order to show the feasibility of the new method based on a fixed grid. Therefore, the simulations only represent the planar elongational flow of fluids filled with particles having a long aspect ratio, such as fibers, which are aligned normal to the plane of flow. For more general flows, like spherical particle suspensions, we need to extend the method to three dimensions. This will be part of future research and will require iterative solvers and parallel calculations.

Numerical simulations are performed and the local flow fields are presented for a many-particle problem. The bulk stress is recovered by using a standard averaging procedure [7]. Finally, the bulk rheological properties are discussed and a comparison with the results of Hwang and Hulsen [3] is carried out. Our results on the bulk viscosity of the suspension are in very good agreement. Moreover, an anisotropic structure is also found even if no transient behavior as in [3] is observed.

The paper is organized as follows: in Section 2, the problem definition is presented. The governing equations for fluid, particles and hydrodynamic interactions are given as well. In Section 3, the weak form for the whole domain is derived. Moreover, the spatial implementation and time integration algorithms are discussed. In Section 4, the bulk stress formula is given. In Section 5, the method is validated. A comparison between the Lagrange multipliers/fictitious domain method and a boundary fitted method is carried out. A simple test problem is chosen. In particular, local flow fields and bulk stress are exploited. The influence of the number of collocation points on the accuracy of the solution is also analyzed. Moreover, the relationship between Lagrange multipliers and traction forces on the particles is discussed. In Section 6, the simulation procedure is introduced. The computational scheme is presented and particle area fraction and bulk stress formulas are given. In Section 7, the results for planar extensional flow are presented. A many-particle problem (150 and 225 particles) is simulated. Local velocity, pressure, stress fields are analyzed and discussed, by means of snapshots of the simulations. Finally, bulk properties (stress tensor and viscosity) are evaluated.