## Abstract

A new method for two-way fluid–particle coupling on an unstructured mesoscopically coarse mesh is presented. In this approach, we combine a (higher order) finite-element method (FEM) on the moving mesh for the fluid with a soft sphere discrete-element method for the particles. The novel feature of the proposed scheme is that the FEM mesh is a dynamic Delaunay triangulation based on the positions of the moving particles. Thus, the mesh can be multi-purpose: it provides (i) a framework for the discretization of the Navier–Stokes equations, (ii) a simple tool for detecting contacts between moving particles, (iii) a basis for coarse-graining or upscaling, and (iv) coupling with other physical fields (temperature, electromagnetic, etc.). This approach is suitable for a wide range of dilute and dense particulate flows, because the mesh resolution adapts with particle density in a given region. Two-way momentum exchange is implemented using semi-empirical drag laws akin to other popular approaches; for example, the discrete particle method, where a finite-volume solver on a coarser, fixed grid is used. We validate the methodology with several basic test cases, including single- and double-particle settling with analytical and empirical expectations, and flow through ordered and random porous media, when compared against finely resolved FEM simulations of flow through fixed arrays of particles.

## 1. Introduction

Fluid flow through particulate media is pivotal in many industrial processes, for example in fluidized beds, granular storage, industrial filtration and medical aerosols. Flow in these types of media is inherently complex and challenging to simulate, especially when the particulate phase is mobile. For the past two decades, particulate flows have been an active area of research and two widely used approaches are now considered state of the art. The first approach is based on an Eulerian continuum model of two phase flows, which only describes the averaged behaviour of the multi-phase media [1]. The second approach is based on an Eulerian–Lagrangian approach using finite-volume/finite-difference methods on a fixed grid as a fluid solver and either an immersed boundary [2], fictitious domain (FD) [3], marker and cell [4] or discrete-element method (DEM) [5] for the particles. Both one-way and two-way couplings have been explored using these methods. While many fluid solvers are based on a staggered grid finite-difference method, others, for example Ladd [6,7], Han *et al.* [8] and Feng & Michaelides [2], have successfully used the lattice Boltzmann method (LBM) as a fluid solver for particle–fluid suspensions. The LBM is an attractive alternative owing to its ease of implementation and parallelization; however, the selection of the lattice resolution, the collision kernel and accurate numerical implementation of various boundary conditions are still challenging [9].

A detailed description of flow through particulate media and accurate particle tracking can be obtained using discrete particle modelling (DPM), as proposed by Tsuji *et al.* [10], Kuipers *et al.* [11], Xu & Yu [5] and Wu *et al.* [12]. In DPM, individual particles are tracked using Newton's laws of motion and particle–particle/wall interactions are taken into account on the (smallest) scale of the contacts. These models invariably couple a continuum solver for fluid with the DEM, as originally proposed by Cundall & Strack [13], for particles. The coupling between fluid and particles is explicit and is achieved using semi-empirical drag laws or closure relations of fluid–particle interactions [14–17]. Yazdchi *et al.* [18,19] proposed modified closure relations (drag laws) applicable to a wider range of porosities for both ordered and random porous media, valid for creeping two-dimensional flows. DPM with hard sphere particle–particle interactions has been successfully applied to, for example, fluidized beds and slug formation in bubbly flows [20]. Those approaches represent a multi-scale method ranging from the smallest scale of the mechanical contacts via the microstructure and particle–fluid interaction (drag laws) to the macroscopic scale of industrial applications. Today's challenges involve developing methods that work for all ranges of densities (from dilute, Knudsen-like regimes to dense sediments or even compressed particle packings) and keeping as much as possible the information of smaller scales (e.g. the microstructure) available on the larger scales. For this the multi-scale, multi-phase model presented in this study features an intrinsic resolution of the fluid that adapts to the particle density and is of the order of the particle–particle distance.

Furthermore, for dense particulate flows, efficient contact detection in a DPM approach requires additional data structures and specialized algorithms adding to its computational overhead. On the other hand, the grid size for flow resolution is often coarse, i.e. they are several times bigger than the mean particle diameter. Thus, most discrete particle models ignore the subgrid scale flow characteristics and this affects the small-scale particle dynamics. Xu *et al.* [21] have recently proposed including subgrid-scale features to better capture the particle dynamics. However, fully resolved methods are often too slow to be able to reach macroscopic scales, so in this study we propose as a compromise using the resolution of the (dynamically changing) particle distance. Note that this, like all the preceding methods [1–5], is based on explicit coupling between fluid- and particle-solvers through empirical drag relations. By contrast, in many finely resolved approaches, an implicit coupling is present. For example, the distributed Lagrange multiplier (DLM) method of Glowinski *et al.* [22] has been successfully applied to simulate fluid–particle interaction in porous media and fluidized beds. Owing to an additional set of Lagrange multipliers, the DLM is more computationally expensive than DPM. However, similar to DPM, the particles are not modelled geometrically in this approach, but the flow in the vicinity of particles is better resolved. Using the DLM, Pan *et al.* [23] simulated the behaviour of fluidized beds; however, they ignored particle–particle interactions to keep the computational costs low. More recently, Kanarska *et al.* [24] have coupled the DLM with the DEM for particle–particle interactions. Fully resolved simulations of particle-laden flows using FD by Avci & Wriggers [3] is in spirit similar to the DLM, except that coupling forces are computed by integrating the stress field at the surface of the particles. In essence, the two methods are exact as no drag correlations are required to couple the two phases.

Interest in using a deforming mesh for fluid–structure/particle interactions has persisted for some time now. Tezduyar *et al.* [25] developed the so-called deforming spatial domain/deforming space–time (DSD/DST)–finite-element method (FEM) for flow problems with deforming interfaces using the so-called arbitrary Lagrangian–Eulerian (ALE) methods and space–time FEM. In this approach, particles are geometrically modelled in the mesh and the flow is fully resolved around each particle, hence leaving it as computationally expensive for dense flows [2].

In this paper, we introduce a new method for fluid–particle interaction based on a two-way coupling between a higher-order FEM and a soft sphere DEM approach on a deforming unstructured mesh. The main feature of our approach is a deforming Delaunay triangulation, which is used as an efficient contact detection tool for the moving particles as well as a finite-element mesh for discretizing the Navier–Stokes equation. It is known that the nearest-neighbour property of the Delaunay edges renders it an attractive algorithm for contact detection (see Ferrez *et al.* [26] and references therein). To better resolve the flow around the particles, we apply the interaction forces as point forces at the particle locations or at their surface (for different version see §3*a*). To the best of our knowledge, this study is the first attempt to apply a moving Delaunay triangulation (particle based) for both contact detection and a finite-element fluid solver. Coupling with the FEM as a fluid solver has advantages; it allows us to build upon very well developed numerical and mechanical methods and thus may provide the leverage of higher-order interpolations for simulating flow to the desired accuracy and scales, even when the mesh is relatively coarse. Another motivation to use the FEM is that, for packed beds and dense particulate flows, the mesh can also be used as a coarse-graining or upscaling tool for stress and strain fields, which are quantities of interest for a macroscale description of the particle–fluid system.

Despite the advantages of using the FEM for higher accuracies, Wu *et al.* [12] have pointed out several issues associated with implementing fluid–particle coupling on an unstructured mesh. The most restrictive one pertains to computing the particle volume fraction in a given cell, because particles may be shared between neighbouring cells and thereby add to the computational complexity. We circumvent this issue in the present methodology by resorting to a moving mesh and an ALE formulation. In this way, the particles have finite radius and always lie at the element vertices and consequently, owing to their repulsive forces and excluded volume, a coarse mesh in a moderately dense system generally remains robust with respect to element degeneration; in other cases, remeshing is deployed whenever necessary, or additional features have to be added to the method. Like the particle contacts, the particle–particle/wall interactions are modelled using a linear spring contact model with dissipation and friction, but the triangulation needs special attention in the vicinity of the walls in order to generate a ‘healthy’ grid.

This paper is organized as follows. We start with an introduction to the mathematical model applicable to viscous, incompressible flow through an isotropic random or ordered porous medium in §2. In addition, the drag force model used for coupling the FEM and DEM and the contact force model used in the DEM are discussed in detail. In §3, we detail the underlying finite-element formulation and discuss the methodology for approximating the porosity field and its impact on numerical computations. This is followed by numerical examples in §4, demonstrating the most basic flow situations in static and moving particulate media. Finally, §5 presents conclusions drawn in this paper and an outlook for future studies.

## 2. Mathematical model

The governing equation for the multi-phase flow is a set of porosity scaled Navier–Stokes equations, which define the flow of fluid in a particulate porous medium (see Anderson and Jackson [27], Deen *et al.* [28] and Xu & Yu [5]). Considering an incompressible fluid (i.e. the density, *ρ* is constant) in an Eulerian flow domain, *Ω*, we can write the equations of both fluid and solid phase as

*Fluid phase*
2.1
*Solid phase*
2.2
where *ε*, *μ*, ** u**,

*p*,

*τ*and

**are the porosity, viscosity, fluid velocity vector, pressure, shear stress and the acceleration owing to gravity, respectively. For the particles**

*g**m*,

*I*

_{i},

*r*

_{i},

*V*

_{i},

*u*_{i}and

*ω*_{i}represent particle mass, moment of inertia, radius, volume, translational and angular velocity, respectively. The represents the inter-particle/wall contact force and is the unit vector pointing from the centre of the particle to the contact point (with particle

*j*). Finally, and represent the drag force per unit volume on the fluid owing to interaction with the

*i*th particle and the total drag force acting on the

*i*th particle, defined in the following section. In the angular momentum equation, represents the torque experienced by the

*i*th particle owing to fluid drag when flow around the particle becomes asymmetric, as shown in §4

*b*. The pressure gradient term in equation (2.2) accounts for the net buoyancy force on each particle acting on its centre of mass. Because equation (2.2) is a system of ordinary differential equations in time, for good accuracy and conservation properties, we use the velocity-Verlet time integrator, which is second-order accurate in time.

^{1}Note that the indices

*i*and

*j*do not represent the tensorial components of respective fields in the above equation, instead

*i*represents particle number and

*j*is the index for the contacts of the

*i*th particle. In the rest of this section, we introduce the model for the drag force density, used to explicitly couple the fluid and particle dynamics.

### (a) Drag force model

The drag force accounts for the resistance to the flow through a porous medium, and is inversely related to its permeability, *K*. The permeability is the proportionality constant in Darcy's equation
2.3
where *U* is the *superficial* (discharge) fluid velocity, relative to the particle speed, defined as
2.4
where *V* and *V* _{f} are total available volume and the volume of fluid. On the other hand, the *intrinsic* average flow velocity is defined only over the fluid volume. Following Yazdchi *et al.* [18], the permeability, *K*, is related to the friction coefficient
2.5
where λ=*K*/*d*^{2} represents the non-dimensional permeability and is often used instead of *K* in the literature. Several existing correlations for λ are listed in table 1. Henceforth, the drag force density in the fluid, , is defined at a point *x*_{e} (see next section for details). The force density is modelled as
2.6
where *u*_{i} is the instantaneous velocity of the *i*th particle and *ψ* is a smoothing function describing the influence of the force density on its neighbourhood. While for *ψ* several possibilities exist, e.g. a Gaussian function, in this paper, we restrict ourselves to *ψ* (*x*−*x*_{e})=*δ*(*x*−*x*_{e}), i.e. the Dirac delta function, for reasons that are discussed in the next section. Equation (2.6) is a model of the drag force density in the fluid in the neighbourhood of the particle. The drag from fluid to particle is proportional to the relative velocity between particle and fluid. In other words, a particle moving in the direction of the flow in its neighbourhood with the average fluid velocity does not experience any drag.

Ergun's equation^{2} is a commonly used drag law, which is a nonlinear function of porosity, fluid velocity and particle size. It accurately predicts the total drag force for a limited range of porosities in three dimensions. Using this relation, one can derive the macroscopic permeability of the medium and use Darcy's equation to determine the average flow velocity through the medium. An aptly modified version of this equation applicable in two dimensions is deployed as suggested in [19]. Accordingly, a more general form of *β*, taken directly from Ergun [14], applicable towards inertial regimes can be written as
2.7
but we do not follow this further and use instead the laminar relation by Yazdchi *et al*. [18] owing to its widest regime of valid porosities. Having specified the drag law, in the following, we introduce a simple contact force model to account for inter-particle/wall forces.

### (b) Contact force model

We take into account the particle–particle/wall interactions and, therefore, the contact forces are essential in order to integrate the particle equations of motion. As elsewhere [31], we use a linear spring-dashpot model for the contact force
2.8
where *κ*, *η*, *v*_{ij} and are contact stiffness, viscous damping coefficient, relative velocity between particle *i* and *j*, and the overlap, respectively. A similar model can also be implemented in the tangential direction along with a sliding spring based on tangential overlap, for cases where rotation and friction are relevant (but is not used in this paper). The contact stiffness, *κ*, and overlap, , set a limit value for the DEM time step as Δ*t*_{DEM}≅1/50*π*/*ϖ* and for numerical simulations. A particle may also have more than one contact at any given time; in this case, the total contact force is found by summing over all the contacts. For further details and state of the art in DEM contact models, see the review paper by Luding [31] and references therein.

## 3. Finite-element formulation

Let us assume we have suitably defined discrete finite-element (polynomial) spaces *V* ^{h}, *S*^{h} for trial and test solutions and let *u*^{h}, *p*^{h} denote the trial solution of equation (2.1). Further, we divide our domain *Ω* into non-overlapping triangles *Ω*_{k} such that ∪_{k}*Ω*_{k}=*Ω*. The weak form is obtained by multiplying equation (2.1) by appropriate test functions (*v*^{h},*q*^{h}) and performing integrating by parts on the diffusion term. This yields a mixed Galerkin formulation for (*u*^{h},*p*), which reads as follows.

Find (*u*^{h},*p*^{h})∈*V* ^{h}×*S*^{h} such that ∀ (*v*^{h},*q*^{h})∈*V* ^{h}×*S*^{h},
3.1
where denotes the standard inner product on *V* ^{h} and *S*^{h}. Note that for the moving mesh we replace the convection velocity *u*^{h} by (*u*^{h}−*u*_{M}), where *u*_{M} is the mesh velocity in the ALE formulation and essentially takes into account the convection of the fluid momentum owing to mesh motion. To compute *u*_{M} at quadrature points inside a triangle, we interpolate the velocities of the nodal particles of that triangle. Using nodal velocities ensures that the geometrical conservation law for the ALE formulation is satisfied [32,33], because a constant solution is reproduced trivially. The above formulation requires *a priori* knowledge of the porosity field at every point inside the domain, which can be computed using smooth particle hydrodynamics (SPH) interpolation (see §3*b*).

For additional robustness and stability in our formulation, we add streamline-upwind/Petrov Galerkin (SUPG), pressure-stabilized/Petrov Galerkin (PSPG) and other terms similar to the least-squares incompressibility constraint (LSIC), as discussed in [34], to the above variational formulation. Förster *et al*. [35] have investigated that such stabilization is also effective when simulating on distorted meshes. Henceforth, following [25] we add residual based stabilization terms (STs) given by
3.2
where ** r** (

*u*^{h},

*p*

^{h}) denotes the residual of the continuity equation, equation (2.1), 3.3 The stabilization parameters are fixed using the following expressions: 3.4 where

*h*

_{e}is the length of the smallest edge of the element. Equation (3.4) has been shown as a convenient choice in computations [34].

Stable discretization of equation (3.1) can be difficult to construct and solutions are well studied in the literature; see [36] for a detailed theory. Classical methods forbid equal interpolation of both velocity and pressure variables in the above setting. Stable solutions can usually be obtained if *P*_{p}⊂*P*_{u} (polynomial spaces for *p* and *u*) in numerical approximations. Furthermore, it is known that the incompressibility constraint is not strongly enforced when using a continuous approximation for the pressure field [37]. To circumvent this problem, we adopt a discontinuous polynomial space for pressure discretization. In this paper, unless stated otherwise, we choose stabilized P1/P0 or P2/P1 elements with continuous velocities and discontinuous pressure polynomials. However, this formulation is not restricted at all in choosing higher-order finite-element spaces.

### (a) The mesh and drag force computation

The finite-element mesh adapted in the above formulation is a Delaunay triangulation based on the particle locations. This implies that all interior vertex nodes of the mesh are occupied by particles at all times, while the boundary nodes are inserted only for the convenience of computation and application of boundary conditions (figure 1).

For moving particles the mesh vertices move with the particles, thereby deforming the mesh. Currently, we re-mesh at fixed (short) time intervals in order to maintain the quality of triangles and also use the triangulation for contact detection. This implies that a new triangulation is created from the current particle positions and the solution from the old mesh is transferred to the new mesh using a simple projection scheme. To remain focused, we will not discuss the projection scheme in detail.

In figure 2*a*, the particle overlaps with an element are shown for particles of different sizes. While we do not address polydisperse particles in a fluid flow in this paper, different size particles are shown to highlight the generality of the proposed method. Figure 2*b* shows the drag force contribution from each element. The total drag force and torque acting on the *i*th particle are considered as a sum of contributions from all the overlap elements
3.5
where *e* is the index counting the number of triangular overlaps, *E*_{e}, of the *i*th particle (e.g. *E*_{e}=5 in figure 2*b*).

An important modelling aspect from the numerical point of view is the location of the drag forces, *x*_{e}, computed from equation (2.6). Here, we list a few possibilities for application of the (as shown in figure 3):

(1) at the midpoint of the chord of the respective overlaps;

(2) at the midpoint of the arc in respective overlaps (as in figure 2);

(3) at the intersection of the particle circumference with element edges; and

(4) at the nodal location of the respective particle.

Figure 3 shows several possible sites for the application of the drag force. Unless specified otherwise, for the sake of simplicity, we choose the midpoint of the chord (i.e. figure 3*a*), as it lies close to the fluid–solid interface, where the momentum exchange occurs. We did not experience an impact of the choice of application of force on the numerical results for the various cases tested, however this issue remains a task for future studies.

A force equal in magnitude but opposite in direction is applied to the fluid, i.e. , at exactly the same point in the cell, thereby providing a consistent point-force-based coupling. The total force on the particle owing to the fluid also contains the buoyancy force, which is computed based on the local pressure gradient.

### (b) Local porosity calculation

At this point, the general variational form, i.e. equation (3.1), can be solved using various assumptions for the porosity field. If the particles are fixed and are relatively homogeneously distributed, then one can simplify equation (3.1) by making the assumption that the porosity is a constant for ∀ *Ω*_{e} and there is no temporal variation. Thus, for a locally averaged formulation, one could take a simpler approach and define a porosity for each triangle in the mesh (see figure 2) as
3.6
Although equation (3.6) is computationally efficient and simple to compute, this definition may lead to high fluctuations in the porosity field, thereby adding to the numerical instabilities especially for dynamic meshes. Therefore, in this paper, we use equation (3.6) only for static particles. To remedy this issue, we interpolate the particle number density using an SPH kernel function as given by
3.7
where *h* is the smoothing length. Following Xu *et al*. [21], one can evaluate the porosity and its gradient at an arbitrary point *r* as
3.8
This definition yields a smoother porosity field on a length scale *h* that can be considerably larger than the particle diameter; however, it incurs additional computation at each numerical quadrature point. Furthermore, special attention is required at the boundaries; for example, Shepard correction is required [38], but details are to be studied in the future.

The issue of the validity of equation (2.1) for different length scales is another open issue for future research. Note that an SPH–DEM code, as introduced in [38], was working reliably well for smoothing lengths as small as *h*∼2*d*. Furthermore, the relevance and validity of drag laws such as Ergun's in the case of moving particles and for rather high resolution is an open issue. For a preliminary study on the local field quantities for different resolution see [39], which considered only static particles but provides vast details on the local and non-local flow behaviour across different scales.

### (c) Time integration

After performing spatial integration, a second-order finite-difference scheme is used for time integration of the resulting system of equations. In a general form, this can be written as
3.9
Using the necessary polynomial approximations of test-and-trial functions, the finite-element matrices for each element in the mesh are assembled and the algebraic form of the equations is written as
3.10
where [*M*] represents the mass matrix, [*C*] is the matrix representing the convection term and [*B*] and [*A*] are the matrices due to pressure gradient and diffusion terms. The [*γ*_{i}] matrix is due to pressure penalty terms on the interior boundaries. The terms in {.} denote the corresponding coefficients of the finite-element solution. We discretize in time using a second-order scheme, i.e. equation (3.9), and the *θ*-method (Crank–Nicolson method with *θ*=0.5) for linearizing the convection term as
3.11
where { *f*} represents the sum of the forces and explicit RHS terms. This implies that the drag forces are explicitly calculated. A suitable time-step size for the FEM is chosen according to the Courant–Friedrichs–Lewy condition and the DEM time step is computed based on the natural frequency of particle contacts. In order to allow that at every fluid time step *n* DEM time steps are performed, the integer *n*=Δ*t*_{FEM}/Δ*t*_{DEM} is specified as the input parameter together with Δ*t*_{DEM}.

## 4. Numerical results

Here, numerical results are presented for both verification and—to some extent—validation of the code. The computational framework described in §3*c* is used to simulate several basic test cases for both static and moving particles. In §4*a*, we first present results for static particles before continuing with deforming mesh simulations. The application to a real complex experimental situation is not shown in this study.

### (a) Static particles

This section deals with flow through static porous media for both ordered and disordered cases. The first example is a simplified model of flow through a homogeneous porous medium, which verifies the compatibility between the present model and Darcy's law. In the second example, we compare our mesoscale-resolution simulation with the average velocities obtained from fully resolved ANSYS simulations of flow through both ordered and disordered arrays of static particles [18,19]. The fully resolved simulations were performed using a fine mesh with approximately 10^{5} elements to accurately capture particle geometry and predict the flow around each particle. Our mesoscale approach, by contrast, contains elements of the same order of the number of particles (i.e. only a few hundreds). While the flow is not fully resolved, the comparison reassures us that this scheme efficiently computes average velocities that are in the expected range and captures qualitatively the flow behaviour at the macroscale while keeping details on the mesoscale.

#### (i) Case 0: homogeneous porous medium and Darcy flow

A well-defined multi-phase model should also reproduce the limit behaviour of single phase flow. When combined with a homogenized drag force (as opposed to point forces), the model must reproduce flow predictions from Darcy's law. As a preliminary verification case, we simulate the flow through a simple porous medium using our formulation with a homogeneous body force (drag) in the test domain. We compare the average flow velocity from simulation with analytical results from Darcy's law. Recall that the permeability, *K*, of the medium describes the resistance to the flow and is intrinsically related to the drag coefficient, *β*, via equation (2.5). Substituting equation (2.5) into (2.3) leads to
4.1
Setting *β*=1 kg m^{−3} s^{−1}, *ε*=0.5 and ∇*p*=−1 kg m^{−2} s^{−2}, one obtains 〈** u**〉=0.5 m s

^{−1}. For this special case, equation (2.1) can be simplified to d

**/d**

*u**t*=(−1/

*ρ*)(∇

*p*+

*β*/

*ε*

**), with**

*u**ρ*=1 kg m

^{−3}. Assuming that the fluid is initially at rest (

**(0)=0), the analytical, transient solution of equation (4.1) is 4.2 We also numerically solve equations (2.1) and (2.2), where**

*u*

*f*^{D}=

*β*〈

**〉 acts as the distributed body force with stress-free boundary conditions. The problem set-up sketched in figure 4**

*u**a*,

*b*shows that the flow quickly achieves a steady-state value of 〈

**〉=0.5 m s**

*u*^{−1}, in perfect agreement with the analytical solution above.

#### (ii) Case 1: flow through ordered and random porous media

For this problem in a square domain, the top and the bottom boundaries have no-slip boundary conditions, whereas the left and right boundaries maintain a pressure gradient of 5 kg m^{−2} s^{−2}. In figures 5 and 6, the colours refer to the horizontal velocity in the ordered and random media, respectively. The dark blue regions (see online version) indicate the slow-flow region between the horizontal rows in the 5×5 particle array, whereas the predominant channel for the bulk flow lies between the two adjacent rows of particles. With decreasing porosity, the flow gradually confines itself between the walls and the top and bottom rows of particles as the interior becomes less and less permeable.

For comparison purposes, the average flow velocity is computed for the entire domain and compared with finely resolved FEM simulations. We used the drag law of Yazdchi *et al*. [18], from table 1, in this simulation, which is valid for a wide range of porosities. The average flow predictions for both ordered and random cases agree very well with data from finely resolved FEM simulations (figure 7). We must mention here that the fully resolved simulation is geometrically correct, i.e. particles are represented by holes with no-slip boundary conditions and contain more than 10^{5} degrees of freedom (d.f.). Our simulation, by contrast, relies on a few hundred d.f. only.

### (b) Moving particles

For the case of moving particles, e.g. fluidized beds, the underlying grid deforms as the particles (and mesh-nodes they occupy) move. This is an important feature of our methodology: because the particle positions are known at all times, it reduces the computational overhead associated with finding particles inside the correct cell [12]. For verification, we present two test cases of one- and two-particle sedimentation. To circumvent the solution degeneracy owing to the deforming mesh, we re-mesh at fixed intervals. Re-meshing is essential in this approach, because we wish to preserve the nearest-neighbour property characteristic of the Delaunay triangulation for contact detection at all times. However, this is not too restrictive as the particles do not move much per Δ*t*_{FEM} time step and contact detection with walls is handled separately in our code. Therefore, we do not address the particles escaping the fluid flow region in the present work, which remains a limitation to address in future work.

#### (i) Case 1: single-particle settling

A particle under gravity in a viscous fluid, both initially at rest, will fall until it has reached the settling/terminal velocity, *u*_{s}, calculated using the drag law prescribed in [40]. The parameters are *μ*=1.14 kg m s^{−1}, *ρ*=1.25×10^{3} kg m^{−3}, *ρ*_{p}=7.74×10^{3} kg m^{−3}, *d*=4.8×10^{−3} m with drag force *f*^{D}=4*πμ**u*_{s}/ln (7.4/*Re*), and *Re*=*ρ**u*_{s}*d*/*μ*.

No-slip boundary conditions are used at the top and bottom walls, while friction-less (no shear stress) boundary conditions are used along the left and right walls. The particle is released from *Z*_{0}=0.6*H* m, where *H*=2 m is the height of the box. The mesh is based on the single-particle location (corner points and two additional boundary points on each wall) and consist of only 12 triangular elements, which is rather coarse. As mentioned before, here, we switch to fourth-order polynomials for an increased flow resolution. The settling velocity can be computed when the frictional force, *f*^{D}, combined with the buoyancy force exactly balance the gravitational force (*m*** g**) and is equal to

*u*_{s}≅0.17 m s

^{−1}. Figure 8 shows the deforming mesh as the particle follows its trajectory. Near the particle surface, a halo region with non-zero upwards fluid velocity appears owing to the drag exerted by the falling particle. A trail of this halo is not evident, because viscosity is large, and our approach does not fully resolve the flow. Note that for this particular case no re-meshing was required as the mesh does not entangle throughout the simulation. The single-particle settling was studied in the context of another mesoscale SPH–DEM coupled method in [38] in much more detail, and we refer to that study rather than testing different particle sizes and different fluid properties here.

#### (ii) Case 2: the drafting, kissing, tumbling problem

We illustrate another benchmark case, where two particles are initially separated vertically and start falling under gravity. As in the previous case, both fluid and particles are initially at rest and particles are then released. The simplest Stokes drag law, i.e. *f*^{D}=3*πμd**u*_{s}, with *μ*=10^{−3} kg ms^{−1}, *ρ*=10^{3} kg m^{−3}, *ρ*_{p}=1.01×10^{3} kg m^{−3} and *d*=4×10^{−3} m, is used in the simulation.

Similar to the previous example, no-slip boundary conditions are used at the top and bottom walls, whereas friction-less (no shear stress) boundary conditions are used on the left and right walls. Figure 9 depicts several snapshots of the two-particle settling behaviour. While the bottom particle centre is aligned with the centreline of the box, the top particle's centre location is offset to the centreline by 1% to the right to trigger the instability (a smaller offset would do the same job; if no offset is used, then the flow situation is determined by numerical errors). As the particles fall through the column of this fluid, the top particle is observed to draft behind the first particle and catches up with the first particle (kissing) and then gets past it with a tumbling behaviour. This behaviour is very sensitive to flow resolution around the particles as the draft of one particle affects the other. This behaviour is well captured in using third-order polynomials for fluid resolution in this approach.

A more quantitative study of this test case and the comparison with fully resolved simulations [41,42] is one of the next steps to be done. The present simple cases are rather an illustration that the method is qualitatively working, but do not yet represent a quantitative validation, which is work in progress.

## 5. Summary and conclusion

A mesoscale, two-way fluid–particle interaction framework based on coupling the FEM and a soft particle DEM on an unstructured, moving mesh has been proposed. The key component in our approach is a Delaunay triangulation, which serves both as a contact detection tool and as an FEM mesh generator. The triangulation deforms and changes with the particle motion. This design alleviates any computational overhead purported by existing methods for contact detections, particularly in dense particulate flows. Because particles always occupy nodal positions in our mesh, locating particles inside cells also becomes trivial while at the same time the particles help to keep the mesh from degenerating. Furthermore, duplication of data for storing the mesh and particles as well as their contact detection has been avoided by defining a triangulation based on particle locations.

An FEM-based fluid solver allows for a higher-order interpolation when a better resolution of the flow is desired, whenever the underlying mesh is coarse. On the other hand, dense flows are resolved rather well with low order, because the mesh resolution is refining inversely proportional to the particle density. Different time scales in the DEM and FEM can be coupled through inner iterations of DEM steps. The approach provides the dynamics of the particles and the fluid using a deforming mesh, while reasonably well resolving the fluid flow around the particles. The average velocities are accurately predicted when compared with fully resolved simulations, which indicates that the assumptions made are valid (for the drag laws, the continuum equations and the choices made for mechanical and technical details) even though they might be at their limits. Many open issues remain for more detailed future studies of the algorithm and methodology, as stated in the previous sections.

The next step is to validate the model with more complex experimental test-cases such as fluidized beds, batch sedimentation or particle dispersion [43], which is work in progress. We have already identified one challenge: the empty areas that are generated by particle inhomogeneities require some type of mesh refinement, because they leave too big cells; this could be addressed by adding ‘ghost’ particles that are used for mesh generation, are repelled from each other and from real particles, but affect neither the fluid nor the other particles. In the future, another challenge is to deal with multiple phases and to couple various other physical fields (e.g. temperature, electromagnetic), using the same data structure and to use the same data structure also for coarse-graining or upscaling [39] to achieve an even more versatile multi-scale model.

## Funding statement

This work is sponsored by STW through the MUST project number 10120.

## Acknowledgements

We thank S. Turek for the discussions related to the incompressible solver.

## Footnotes

One contribution of 13 to a Theme Issue ‘Multi-scale systems in fluids and soft matter: approaches, numerics and applications’.

↵1 Because the forces between particles can be dissipative, the choice of an integrator does not have a major impact either on solution quality or on the performance, thus will not be discussed in detail in this paper.

↵2 The Ergun equation is essentially a correction to the Carman–Kozeny [29] drag relation for creeping flows, which also takes into account the inertial drag at higher Reynolds numbers [30].

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.