 Home Engineering  # Applications: Wave and sediment grain interaction by a nonbreaking solitary wave on a steep slope

The second application case reproduces the interaction of a solitary wave with solid particles representing the sediment grains of a beach, modelled by means of the Discrete Element Method (DEM). Particles can collide and interact between them and with the walls, and are coupled with the flow via the individual drag forces and dynamic porosity terms included in the VARANS equations. The base setup is similar to the one depicted in Higuera et al. (2018). The goal of this example is to observe the evolution of the sediment grains in the swash zone and compare the numerical results qualitatively with observations available in literature.

## Introduction to DEM

The purpose of this subsection is to serve as a brief introduction to the physics and implementation of DEM modelling coupled with CFD. Due to practical space limitations, only the basics of DEM will be covered. The reader is referred to Kloss et al. (2012); Smuts (2015); .Ting et al. (2016); Norouzi et al. (2016) for further details concerning the underlying theory.

6.1.1 Overview

The Discrete Element Method, also known as Distinct Element Method, is a numerical method that has been widely used in literature to solve discontinuous mechanics problems such as granular flow or rock mechanics simulations since Cundall and Strack (1979). More recent advances include the so-called extended discrete element method (XDEM) (Peters, 2013), which enhances DEM, allowing solving multi-physics problems such as chemical reactions (e.g., pellets acting as catalysts), heat transfer or electro-magnetic forces (e.g., magnetic particle sorting and separation). Although DEM particles are often spherical in shape, this is not a restriction, as complex shapes can be built from overlapping spheres of different sizes (Ferellec and McDowell, 2010). There are two main approaches for sphere interaction. In the inelastic hard sphere method, the spheres are rigid and particle collisions produce an exchange of momentum through impulsive forces. This method has a limitation, it does not allow sustained contact between particles, therefore, it is suitable for the so-called collisional particle regime only. The other method is the soft sphere approach, which is more common in DEM. In it, particles, although still assumed rigid, are allowed to overlap slightly. This way, the interaction between particles takes a finite time and allows simulating a sustained contact between particles (i.e., frictional forces) as well as regular collisions (Mitarai and Nakanishi, 2003). The soft sphere method will be used in the present work.

In a more flexible approach, DEM has also been coupled with Finite Element Method (FEM) (Munjiza, 2005), introducing tetrahedral meshes to represent arbitrary shapes and enabling calculating internal stresses within the elements (Latham et ah, 2008). The same research group applied the FEM-DEM to model rubble mound coastal structures in Latham et al. (2013); Anastasaki et al. (2016), although their model was not coupled with CFD at that time, and their investigations were limited to seismic-type loading, as they were exciting armour units using a global acceleration field.

DEM coupling with CFD has been recently gaining momentum for research and industrial applications alike, propelled by open source models. For example, Klossr and Goniva (2010) developed CFDEM, a coupling between the CFD solver OpenFOAM® and the DEM solver LIGGGHTS (originally derived from LAMMPS), which has been widely used in literature ever since (e.g., Zhao and Shan (2013); Schmeeckle (2014)). Other references with a specific focus on modelling multiphase flows with the VOF technique include Xiang et al. (2012), .ling et al. (2016) and de Lataillade et al. (2017).

Two approaches exist for coupling CFD and DEM. In the resolved coupling, the particles are larger than the underlying CFD mesh cell size, allowing the detailed flow around the particles to be simulated. This technique requires the immersed boundary or fictitious domain method to be implemented in the CFD solver, so that the cells inside the large particles are blocked, acting like the solid part of a moving obstacle, and the cells intersected by particle boundaries act as moving walls (Guo et al., 2013). Alternatively, in the unresolved coupling, particles are smaller than the cell size and the detailed flow around them is not simulated. Instead, the particles are assumed to be a point with a unique set of flow magnitudes (i.e., VOF, velocity, pressure...)

The benefits of using DEM for sediment transport is that loose granular flow can be represented in a more physical way, rather than as a continuum, provided that particle sizes are large enough. Also, given the physical properties of the material and grain shape, the critical angle of repose and the avalanching processes are captured implicitly. It must also be noted that the main disadvantages of DEM are the large computational costs because of the large number of particles that may be involved in simulations and the challenges posed to parallelization when coupled with mesh decomposition techniques. Modern models are parallelized or even GPU-accelerated (Qi et al., 2015), thus reducing some overheads, but there is still room for improvement.

In this work, the DEM model developed in OpenFOAM®, consisting of solid spheres, will be used to simulate individual sediment (sand) particles. Sediment transport studies are common in the DEM literature. Although most of them consider a single phase flow (Schmeeckle, 2014), they cannot represent a moving shoreline, but they can include various types of particle shapes (Sun and Xiao, 2017). Moreover, the native DEM model in OpenFOAM® has previously been validated in a backward-facing step setup for a single phase flow (Greifzu et al., 2016).

6.1.2 Particle physics

In DEM every particle is defined by its physical properties (e.g., shape, inertia, material...) and its position, orientation and linear and angular momentum (i.e., velocity and rotation velocity). Under certain conditions, some of the previous information can be neglected, as for example the orientation of the particles when considering the soft sphere approach.

In this particular case in which the particles are considered rigid spheres, the shape will not change in time and the contact points with their neighbours and the walls will be straightforward to calculate. In spite of being considered rigid, the particles are allowed to overlap, thus through contact models, the forces acting on each particle can be calculated. Particles follow Newton’s 2nd law, therefore, acceleration (linear and angular) can be easily calculated given the forces and the properties of the particle, as shown in the following equations: in which m is the mass and I is the moment of inertia of the particle, v is velocity, oj is the angular velocity. Y F represents the sum of all forces acting on the particle, which include the gravitational force, buoyant force, drag force, contact forces (normal and tangential), added mass force and/or Magnus (lift due to rotation) and Saffman (lift due to shear) forces. Y M represents the sum of all moments (torques) acting on the particle.

Given the previous equations, the velocity and position of the particle in the new time step can be obtained directly by numerical integration.

6.1.3 Numerical implementation

In this subsection we will first introduce the particulars regarding the numerical implementation of the unresolved DEM model in OpenFOAM® and the coupling of DEM with the two-phase flow solver.

The DEM physics has been integrated in a solver derived from olaFlow, in which the particle evolution loop has been included within the main solving loop before advecting the VOF function. Since this is an unresolved method, flow properties of a particle are taken from the cell in which the centre of the particle is located.

The particle effects have been added to the VOF, conservation of mass and conservation of momentum equations as explained in Section 3.3. In this case the contribution to the different cells, in the form of porosity and source terms, is proportional to the volume of the particle inside each cell over the volume of the cell.

Regarding the forces acting on the particles, gravity and buoyancy are straightforward to consider, as they only depend on the mass and volume of the particle and the density of the fluid in which it is submerged. There are two types of forces that are modelled: drag and contact forces, which are the most significant. The rest of the forces, i.e., added mass, Magnus and Saffman, are not taken into consideration due to the unresolved approach used, although they are expected to be of second order importance in this case.

The drag coefficient {Cp) is modelled according to a complex expression that fits the experimental data for a single smooth sphere by Schlichting (1979). The equation, presented in Morrison (2013), is as follows: where Re is the Reynolds number. This model can be used from creeping flow (Re ~ 10-1) up to a turbulent value of Re ~ 10°.

Contact forces account for continuous contacts (e.g., particle resting or rolling on a wall) and for instantaneous contacts (e.g., collisions). These forces are modelled with the spring-slider-dashpot model, which takes into account the normal force, the friction force and the damping of energy, respectively. The inputs for the model are the Young’s modulus and shear modulus, the coefficients of friction (tangential sliding) and restitution and the cohesion energy density, to model cohesive materials. For further reference the reader is referred to .Tiang et al. (2009).

Solving the equations for particle motion (20-21) has a lower complexity as compared with solving Navier-Stokes equations. However, the time step of the DEM model needs to be much smaller than that of the CFD to limit the amount of overlapping between particles. Otherwise, large overlapping would yield excessively large forces. Moreover, short time steps ensure that the energy is transferred at a maximum of the Rayleigh wave speed (Otsubo et al., 2017). Therefore, for every CFD time step, which is on the order of 10~!, the DEM solver can calculate up to thousands of sub-time steps.

 Related topics