MOTION OF SOLID SPHERES UNDER SOLITARY WAVE ATTACK : PHYSICAL AND NUMERICAL MODELING

In this study, physical and numerical test are carried out focusing on the motion of two solid spheres under solitary wave attack. The numerical model CADMAS-SURF/3D-2F-DEM coupling a multiphase flow solver solving Reynolds Averaged Navier-Stokes Equations based on a porous body model and a discrete element method solver for Newton’s equations of motion is validated against the data of physical model experiments carried out in the wave flume of METU Ocean Engineering Research Center. Comparisons of the numerical simulations and physical model experiments show that the numerical model is capable of simulating the motion of solid spheres under solitary wave attack in a reasonably well accuracy.


INTRODUCTION
Numerical simulations using computational fluid dynamics (CFD) tools coupled with discrete element method (DEM) solvers for engineering purposes gained popularity in the last two decades, as a well-developed CFD-DEM solver is capable of simulating the motion of solid particles under flow conditions.Among other studies in the literature, there are several studies in the scope of coastal engineering applications focusing on sediment transport (Sun and Xiao, 2016), breakage of armor units of rubble mound breakwaters (Latham et al. 2013), and the motion of solid object under flow conditions (Arikawa et al., 2011;Canelas et al., 2016).
CFD-DEM solvers could be used to assess motion of rubble stones under wave attack.Most of DEM studies start modeling the solid particles as spheres since the contact detection algorithms for spheres are trivial.Therefore, accurate numerical modeling of motion of solid spheres under wave attack is vital as the first step towards modeling motion of rubble stones.
In this study, the motion of two solid spheres under breaking solitary wave attack are physically and numerically modeled.Physical model experiments are carried out in the wave flume of METU Ocean Engineering Research Center.On the other hand, CADMAS-SURF/3D-2F-DEM is used as the numerical model.The major aim of the study is to extend the understanding of the physical processes involved in the motion of the solid spheres under breaking solitary wave attack and to calibrate/validate the numerical model used in the present study with a simple experimental configuration.
The present paper is structured as follows: After this brief introduction, physical model experiments and analysis technique for tracking the spherical particles are summarized.Next, the numerical model and numerical simulations are described.Results of the numerical simulations in comparison with the experimental data are presented and discussed.Finally, summary and conclusions are given.

Experimental Setup
Physical model experiments are carried out in the wave flume of METU Ocean Engineering Research Center.Length and width of the wave flume are 26 m and 0.8 m, respectively.A slope of 1:20 is placed in the wave flume followed by a horizontal flat area.On the flat area, two spherical billiard balls colored as red and yellow with a diameter of 6 cm are placed setting the center to center distance between spherical particles as 8 cm.The density of the billiard balls is 1.9 g/cm 3 .The water level is set as 33 cm; in other words, the water level is up to half of the billiard balls.The water level inside the wave channel is measured by five wave gauges, and the motion of spheres are tracked using a camera placed at the top of the flat area.A sketch of the experimental setup is given in Figure 1.In the experiments, a solitary wave with a wave height of 6.6 cm at WG1 is used.That wave reaches to 7.8 cm at WG5 and breaks on the billiard balls.The experiment is repeated for five times.

Analysis of Video Recording using Image Processing Techniques
The experiments are recorded with a camera (SONY Cybershot DSC-HX7V) at 25 Hz placed at the top of the experimental setup to track the motion of the spheres.A checkerboard is placed on the flat area to ease the tracking procedure.The paths of the spheres are tracked using color detection algorithms implemented in the MATLAB environment.The algorithm differentiates yellow and red regions in the images based on Hue-Saturation-Value (HSV) at each pixel of the image.The centroids of these regions are calculated to find the center of the spheres in x and y coordinates in a reasonable accuracy.Actual and tracked positions of the spheres at t=11.72 sec and t=13.00 sec are given in Figure 2. In Figure 2, the yellow sphere is tracked by the pink point whereas the red sphere is tracked by the black point.Analysis of the physical model experiment revealed that no collision occurs.The breaking solitary wave pushes the yellow ball (at the right-hand side) faster than the red ball (at the left-hand side).Therefore, the distance between the balls increases with respect to time.

Description of the Numerical Model
The numerical model CADMAS-SURF/3D-2F-DEM is developed by Arikawa et al. (2011) coupling a flow solver called CADMAS-SURF/3D-2F with a Discrete Element Method solver called DEM.
The flow solver CADMAS-SURF/3D-2F solves multiphase Reynolds Averaged Navier-Stokes Equations (RANS) based on a porous body model.Furthermore, the effect of the air compression is taken into account with a specific term that is introduced in the continuity and momentum equations.The continuity equation is given in Equation 1 and the momentum equations are given in Equations 2-4.
In Equations 1-4, x, y and z are horizontal, transverse and vertical axes; u, v, and w are the velocities in x, y and z-direction, respectively; t is time, ρ is the density of water; ρG is the density of air; G is the derivative of the density of air; p is pressure; F is the ratio of fluid volume in a cell to cell volume; νe is the effective viscosity (summation of the molecular viscosity, ν and eddy viscosity, νt); γv is the porosity; γx, γy and γz are permeability in x, y and z directions, respectively; Dx, Dy and Dz are coefficients for energy dissipation in x, y and z directions, respectively; Sρ, Su, Sv and Sw are source terms for wave generation and g is the gravitational acceleration.λv, λx, λy and λz are expressed by using the inertia coefficient CM and permeability coefficients as given in Equation 5. ( Rx, Ry and Rz are resistance terms for a porous structure where dx, dy and dz are the mesh sizes in x, y and z directions, respectively, presented in Equation 6. To consider the effect of the air phase, the air phase densities are temporally and spatially distributed.The derivative of the air phase ( G ) in Equations 1-4 is defined in Equation 7 as an advection term.
In CADMAS-SURF/3D-2F, the free-surface is tracked using Volume of Fluid (VOF) method in the flow solver by solving the VOF-advection equation given in Equation 8where SF is a source term for wave generation.
In this flow solver, turbulence is addressed using k- turbulence closure.The equations presented above are solved using the finite difference method on staggered grids.Wave generating boundary conditions are available in the numerical model.Furthermore, energy dissipation zone (sponge layer) can be defined in the numerical model to construct an almost non-reflecting boundary.
The discrete element method solver DEM solves Newton's equations of motion given in Equations 9 and 10.In Equations 9 and 10, m is the mass of the particle; u pi is the velocity of the particle; F pi is the forces acting on the particle; t is time; I is moment of inertia of the particle;  is the angular velocity of the particle, and T pi is the torque acting on the particle.Forces acting on the particles are defined as the forces due to the fluid motion, collision forces between the particles and the reaction forces from the walls (boundaries) of the computational domain.
In this numerical model, the forces due to the fluid motion are calculated using the wave pressure.DEM is originally developed for solid bodies represented by several spherical particles, and these particles are connected with bonds.The model evaluates the forces acting on the solid body to assess damage by "breaking" the bonds in between particles.In the present study, the motion of the two solid spheres is studied.Therefore, there are no bonds in between the particles that would be broken by forces applied to the solid body.DEM model is slightly modified to account for the pressure forces at the surrounding cells of the particle.The modification is simply interpolating pressure forces at the surrounding cells on to the particle itself.A sketch of the modification in the x-direction is presented in Figure 3, and the pressure force updating equation in the direction of motion is given in Equation 11 for x-direction.The extension to other dimensions is straightforward.
In Equation 11, fx is the pressure force in the x-direction, f1 is the interpolated force at the negative x-direction, and f2 is the interpolated force at the positive x-direction.
Forces arising from the collision of the particles and reaction of the walls are modeled with an elastic spring and viscous dashpot system given by Tsuji et al. (1993).Coefficients of the spring-dashpot system could be calibrated and validated with the use of experimental data.
The flow solver and the DEM solver are coupled by transfer of information between these solvers.DEM particles in the computational domain are represented by porous regions.Information on the flow properties around these porous regions are sent to DEM solver, DEM solver evaluates information to calculate the forces acting on the particles and the new positions of the particles, and finally, DEM sends the new positions of the particles to the flow solver.

Numerical Setting for the Simulation
Side view of the computational domain is given in Figure 4.The length (x-dir), height (z-dir) and width (y-dir) of the computational domain are set as 18 m, 0.5 m, and 0.4 m, respectively.In all directions, grid size is selected as 2 cm resulting in 1025000 cells.A 2.5 m long energy dissipation zone (sponge region) is placed at the end of the channel to reduce the effect of reflecting waves.Inlet boundary condition is set as a theoretical water surface and velocity profile of a solitary wave with a wave height of 6.6 cm.Simulation duration is selected as 15 sec whereas the sampling interval is 0.05 sec.In this simulation, k- turbulence model is used.

Results and Discussions
The numerical model is validated against the experimental data by changing spring and viscous dashpot coefficients which can be referred to as DEM parameters.The parameters used in the validation simulation are presented in Table 1.In this case, as there is no collision between the particles, these DEM parameters are only used to simulate reaction forces by the flat bottom where the balls move on.Snapshots from the simulation are presented in Figure 5.As it can be observed from the figures, the solitary wave breaks on the particles, and the water particle velocities at the front of the wave are in the order of 0.5 m/s.Similar to the physical model experiments, no collision occurs in the numerical simulation.The yellow ball (on the side in the numerical simulation) is moves faster due to the wave action.The distance between the particles increases with respect to time.Water surface elevation is measured using five wave gauges along the wave channel presented in Figure 1.Comparison of water surface elevations from the physical model experiments and the numerical model are given in Figure 6.In Figure 6, it is seen the peak water surface elevations along the wave channel are captured in reasonably well agreement with the experimental results.However, there are clear differences especially in the second half of the solitary waves.These differences are due to the inlet boundary condition.As the theoretical water surface elevation and the water particle velocities are forced as the boundary conditions, deviations from the experimental data are observed.When the actual water surface elevation and water particle velocity data were used, these deviations would be minimized (Guler et al., 2018).
Next comparison between the experimental data and the results of the numerical model is on the xcoordinates of the red and yellow spheres presented in Figure 7. Y-coordinates of the spheres are not discussed in this case as the movement in that axis is negligible in both experimental and numerical results.In Figure 7, it is observed that the x-coordinates of the particles are captured in a good accuracy for both spheres.However, the results can still be improved by using the measured water surface elevation and water particle velocity time-series data from the experiments as the inlet boundary condition.The results from the numerical modeling study are compared with the data from the physical model experiments.It is seen that the water surface elevations and the x-coordinates of the spheres are captured in a reasonably well agreement.However, improvements in the results could be achieved using the measured water surface elevation and water particle velocity data as the inlet boundary condition.

Figure 1 .
Figure 1.Experimental Setup (Figure is not to scale!)

Figure 2 .
Figure 2. Actual and Tracked Positions of the Spheres

Figure 3 .
Figure 3. Sketch of the Modification in X-Direction

Figure 4 .
Figure 4. Actual and Tracked Positions of the Spheres Figure 5. Snapshots from the Simulation

Figure 6 .
Figure 6.Comparison of Water Surface Elevations of the Numerical Model (dashed orange lines) with the Experimental Results (solid blue lines) with respect to Time