Simulation of flow past a sphere on a rough bed using smoothed particle hydrodynamics (SPH)

Abstract

This paper presents an investigation of flow past a sphere on a rough surface by means of simulation using the meshless numerical method, smoothed particle hydrodynamics (SPH) in the presence of a free surface. A representative of natural and engineered objects is spherical wall-mounted bodies. These are ideal to study the flow conditions around hydraulic structures, boulders, fish habitat structures or even architectural structures. The aim is to understand their effect on the hydrodynamics around the larger object sphere resting on them and their interaction. The Lagrangian particle-based scheme SPH using the open-source code DualSPHysics is validated against datasets from a laboratory-based flume experiment. The validation case was selected from the family of the flow past a sphere test cases since it reproduces the flow features in a simplified manner that can be observed in the vicinity of natural sediment grains or larger bodies. The validation results include flow velocity profiles in the vicinity of the large sphere and comparison with experiment data. The results stand in overall agreement with the experimental velocity and force measurements that demonstrates the applicability of SPH in aquatic environments. An SPH investigation on a rough bed in combination with objects/body near the bed in the presence of a free-surface flow has not been shown before in the literature. This is a novel application with insight into the fluid mechanics made possible by using DualSPHysics solver.

1 Introduction

A commonly used representative of natural and engineered objects is spherical wall-mounted bodies. These are ideal to study the flow conditions around boulders in fluvial beds, fish habitat structures or even architectural structures [1]. Even if the geometry of a sphere appears simple, the flow field around it is a complex, three-dimensional, unsteady and still not fully understood. In addition, experimental studies conducted with a wall-mounted sphere are rare [1,2,3]. Many physical and numerical studies [2, 4,5,6] have examined the flow around such a body. Typical flow features that can be observed around bluff bodies are the flow separation in the vicinity of the sphere and the development of a wake in the rear of such an object body [4,5,6]. Investigations in the laboratory are usually the first choice to quantitatively derive parameters with various measuring devices. To assess the capability of numerical simulations, we performed a set of numerical experiments to investigate the effects of a single sphere on a rough bed, i.e. bluff body on the prevailing hydrodynamics. The acquired data are used to validate a 3-D numerical model against a dataset of Papanicolaou et al. [2].

Numerical simulation techniques that have been used in the past to investigate the flow past a sphere include large-eddy simulations (LESs) and Reynolds-averaged Navier–Stokes equations (RANS) [7, 8]. The general idea of LES is to reduce the computational cost by modelling the smallest length scales, which are the most computationally expensive to resolve, via low-pass filtering of the Navier–Stokes equations. In the case of RANS, the instantaneous quantity is decomposed into its time-averaged and fluctuating quantities. This approach can be used with approximations based on knowledge of the properties of flow to approximate time-averaged solutions of the Navier–Stokes equations [7, 9]. Liu et al. [8] used a LES to analyse the flow over submerged boulders that were placed on a rough bed. However, these Eulerian techniques require a mesh, and so the use of these methods can be inefficient due to inherent difficulties of mesh generation and the relatively high computational cost in particular, when it comes to the meshing of cases in the presence of a free surface and where the flow needs to be resolved downstream of an obstacle.

An alternative to Eulerian mesh-based schemes are meshless Lagrangian schemes where the computation points are particles which move according to the governing dynamics. Smoothed particle hydrodynamics (SPH) is a relatively young computational fluid dynamics simulation technique. Unlike most Eulerian numerical methods (e.g. finite differences (FDM) or finite volumes (FVM)), SPH does not need an interconnected mesh or grid to discretize the domain in order to calculate the spatial derivatives of the governing equations [10,11,12,13]. Instead, the physical interaction between the individual particles is solved through a three-dimensional convolution-based discretization [11, 14]. Since SPH is a meshless method, it is ideally suited to simulate problems dominated by complex boundary dynamics and free-surface flows which are the ultimate focus of this research.

The SPH solver used in this study, namely DualSPHysics [15, 16], is open-source solver that exploits the computation acceleration provided by a graphics processing unit (GPU) for the calculation of the hydrodynamic equations. This parallel architecture is shown to increase the computation speed by orders of magnitude [17,18,19]. Therefore, relative complex simulations can be performed on a single-user workstation in a reasonable time frame without demanding access to HPC clusters or other high-performance facilities. SPH has been rigorously tested and validated for a large range of scenarios, and it is capable of simulating accurately hydrodynamics and rigid body interactions [10, 12, 13, 15, 20,21,22,23].

It needs to be emphasized that although flow over rough bed has been investigated before in SPH [24,25,26], to the best of our knowledge none of these studies include a rough bed in combination with objects/body near the bed in the presence of a free-surface flow. This is new and novel application with insight into the fluid mechanics using DualSPHysics. Therefore, our initial aim is to test the validity of the SPH model as the basis for follow-on studies for flows past bluff objects on rough beds with a free surface. In this study, our numerical experiment will investigate and validate the flow past a stationary sphere in resting on a rough surface that can be considered as flow past a boulder [5]. This ultimately allows for a comparison of the numerical model measurements and the flume data of Papanicolaou et al. [2]. The primarily goal of this research effort is to examine the suitability of SPH for such flows and its application to near bed granular flows for fluvial applications which involve free-surface flows. The numerical results are compared with the results obtained by a laboratory-based flume experiment of Papanicolaou et al. [2] and hence serve as the basis for similar approaches such as the flow past a non-spherical object to be used with SPH in the future.

Although the flow of the experiment of Papanicolaou et al. [2] is reported to be turbulent, in this work we use the artificial viscosity formulation [27] which involves a Laplacian operator and a free parameter where we tune our artificial viscosity using an equivalent Reynolds number based on the artificial viscosity formulation instead of a RANS or LES approach. As discussed by Meringolo et al. [28] in relation to the artificial viscosity term in SPH, when $Re_{dp} = O\left( 1 \right)$, where $dp$ is the particle spacing, the artificial viscosity term resolves the main vorticity scales for that specific ratio of sphere diameter to particle spacing where the Mach number $Ma \approx 0.1$ or lower and $Re_{dp} \approx 1$. Thus, the artificial term has been used to approximate the experimental Reynolds number using the $Re_{dp}$ which for the finer resolution is $Re_{dp} = 1$ herein. This modelling simplification may not reproduce flow characteristics near the walls and downstream to high accuracy, however, has significant advantages in computational cost and ease of use.

This paper is structured as follows: Sect. 2 presents the SPH methodology and numerical set-up; Sect. 3 presents the numerical results, followed by a discussion in Sect. 4; and finally, Sect. 5 presents the conclusions obtained from this work.

2 SPH methodology and numerical set-up

2.1 The SPH formalism

The numerical technique SPH was used to model the present 3-D validation case. According to Gingold and Monaghan [29], the SPH formulation is in principle a kernel or integral approximation of a function $f$ that can represent a physical or numerical variable defined over a domain Ω at a point x,

$$ f\left( x \right) \approx \mathop \int \limits_{{\Omega }}^{{}} f\left( {{\mathbf{x}}^{^{\prime}} } \right)W\left( {{\mathbf{x}} - {\mathbf{x}}^{^{\prime}} ,h} \right){\text{d}}x^{^{\prime}}, $$

(1)

where h represents the smoothing length that defines the size of the support domain of the kernel and W is the weighting or kernel function. The kernel function W is chosen to be a smooth, isotropic and even function with compact support (with a finite radius of influence around x). Similar to Fourtakas and Rogers [30], the fifth-order Wendland kernel with compact support of 2 h was used [31]

$$ W\left( {R,h} \right) = a_{d} \left( {1 - \frac{R}{2}} \right)^{4} \left( {2R + 1} \right), $$

(2)

where R is the distance $\left| {{\mathbf{x}} - {\mathbf{x}}^{^{\prime}} } \right|$. The constant $a_{d}$ represents the normalization constant that is $3/4h$, $7/4h\pi^{2}$ and $21/16h\pi^{3}$, in the 1-D, 2-D and 3-D, respectively.

Within a defined domain of interest, Eq. (1) can be approximated by using an SPH summation of the following form:

$$ f\left( {\mathbf{x}} \right) = \mathop \sum \limits_{j}^{N} f\left( {{\mathbf{x}}_{j} } \right)W\left( {{\mathbf{x}} - {\mathbf{x}}_{j} ,h} \right)V_{j}, $$

(3)

where $f_{j} = f\left( {x_{j} } \right)$and V is the volume of the Lagrangian particle that is expressed as the ratio of mass $m$ to density $\rho$. The number of particles within the compact support is denoted by $N$. In the present study, the subscript $i$ represents the interpolating particle and $j$ the neighbouring particles. The symbol $\ldots$ stands for the SPH interpolation and will be dropped for brevity in the rest of the manuscript. The resulting form of the particle approximation in a discrete form reads

$$ f\left( {x_{i} } \right) = \mathop \sum \limits_{j}^{N} \frac{{m_{j} }}{{\rho_{j} }}f_{j} W_{ij} , $$

(4)

where $W_{ij} = W\left( {x_{i} - x_{j} ,h} \right)$. The reader is referred to Gingold and Monaghan [29], Violeau [32] or Violeau and Rogers [13], on more details regarding the SPH formalism. The flow characteristics within the domain are described by the Lagrangian form of the Navier–Stokes equations discretized using the SPH scheme. Using superscripts α and β to denote coordinate directions employing Einstein's summation, the continuity and momentum equations in Lagrangian form can be written as

$$ \frac{{d{\uprho }}}{dt} + { }\rho \frac{{\partial u^{{\upalpha }} }}{{\partial x_{{\upalpha }} }}{ } { } = { }0 $$

(5)

and

$$ \frac{{du^{\alpha } }}{dt} = \frac{1}{\rho }\frac{{\partial \sigma^{\alpha \beta } }}{{\partial x^{\beta } }} + g^{\alpha }, $$

(6)

where $u$ denotes the velocity vector, $g$ the gravitational acceleration and $\sigma$ the total stress tensor. In a fluidic approach, as in the present case the total stress tensor can be written as the isotropic pressure $p$ and the viscous stresses $\tau$, this results in the following form:

$$ \sigma^{\alpha \beta } = - p\delta^{\alpha \beta } + \tau^{\alpha \beta }. $$

(7)

The continuity and momentum equations can be written in the following SPH form [32]

$$ \frac{{d\rho_{i} }}{dt} = \rho_{i} \mathop \sum \limits_{j}^{N} \frac{{m_{j} }}{{\rho_{j} }}\left( {u_{i}^{\alpha } - u_{j}^{\alpha } } \right)\frac{{\partial W_{ij} }}{{\partial x^{\alpha } }} + D_{i} $$

(8)

and

$$ \frac{{du_{i}^{\alpha } }}{dt} = \mathop \sum \limits_{j}^{N} m_{j} \left( {\frac{{p_{i} + p_{j} }}{{\rho_{i} \rho_{j} }} + \Pi_{ij}^{\alpha } } \right)\frac{{\partial W_{ij} }}{{\partial x^{\alpha } }} + g_{i}^{\alpha } . $$

(9)

The last term on the right-hand side of Eq. (8) is the density diffusion term as described by Molteni and Colagrossi [33], which reads,

$$ D_{a} = 2\delta hc_{0} \mathop \sum \limits_{j}^{N} \frac{{m_{j} }}{{\rho_{j} }}\left( {\rho_{j} - \rho_{i} } \right)\frac{{x_{ij}^{a} }}{{\left| {x_{ij} } \right|^{2} }}\frac{{\partial W_{ij} }}{{\partial x^{\alpha } }}, $$

(10)

which is a second-order density filter necessary to remove spurious density and consequently pressure oscillations. The term $\delta$ is a free parameter usually set to $\delta = 0.1$. Herein, the artificial viscosity formulation of Monaghan [11] has been used to stabilize the solution and provide an artificial viscous force denoted as $\Pi_{ij}$ in Eq. (9),

$$ \Pi_{ij}^{\alpha } = \left\{ {\begin{array}{*{20}l} {\frac{{ - \alpha_{\pi } \overline{c}_{ij} hu^{a} x^{a} }}{{\overline{\rho }_{ij} \left| {x^{a} } \right|}} \quad u^{a} x^{a} < 0 } \\ \,\,\,\,\,\,{ 0 \quad{\text{otherwise}}} \\ \end{array} } \right., $$

(11)

where $\alpha_{\pi }$ is the free parameter and the overbar denotes the average values of the i and j particles such that $\overline{c}_{ij} = \frac{1}{2}\left( {c_{i} + c_{j} } \right)$ and $\overline{\rho }_{ij} = \frac{1}{2}\left( {\rho_{i} + \rho_{j} } \right)$ are the average speed of sound and density, respectively. Herein, we do not use a dedicated turbulence modelling approach for the turbulent stresses but instead model the flow by tuning the $\alpha_{\pi }$ free parameter to the flow characteristics through an artificial equivalent Reynolds number.

Furthermore, in this work the weakly compressible SPH (WCSPH) approach was used to link pressure and density by using the Tait’s equation of state, resulting in a weakly compressible fluid:

$$ p = B\left( {\left( {\frac{\rho }{{\rho_{0} }}} \right)^{\gamma } - 1} \right), $$

(12)

where $\rho_{0}$ is the reference density and $B$ is related to the compressibility of the fluid which is proportional to the speed of sound:

$$ B = \frac{{C_{s0}^{2} \rho_{0} }}{\gamma }, $$

(13)

where $\gamma$ is the polytrophic index and was set to 7 in the present study. The variable $C_{s0}$ is referred to the speed of sound that is computed as:

$$ C_{s0} \ge 10u_{max}, $$

(14)

where $u_{max}$ is the maximum velocity magnitude expected in the domain. The resultant Mach number from Eq. (13) is 0.1 which allows for 1% compressibility in density. Further information about the WCSPH approach can be found in [27]. Apart from the boundary conditions which are described in more detail in the following, the above equations are integrated in time via an explicit second-order predictor–corrector integrator scheme. The corresponding values are corrected using the forces at half time steps, followed by the evaluation of the values at the end time step (Gomez-Gesteira et al. 2012). Moreover, the scheme is bounded by the CFL condition, the maximum force term and the numerical speed of sound as shown in Monaghan [34]. The CFL condition is written as:

$$ \Delta t = C_{0} \mathop {\min }\limits_{i} \left( { \sqrt {\frac{h}{{\left| {F_{i} } \right|}},} \frac{h}{{C_{s0} }}} \right), $$

(15)

where $F_{i}$ is the force per unit mass of the particle $i$ and $C_{0}$ is the Courant number that is set to 0.2 in this study.

In this work, we use the so-called particle shifting algorithm to improve the particle anisotropy which leads to numerical error in the SPH approximation as described by [35] in incompressible SPH and lately in weakly compressible SPH [35]. The shifting algorithm shifts particles to new positions in order to improve the particle distribution based on a Fickian approach, which reads,

$$ \delta x_{i}^{a} = - D\nabla C, $$

(16)

where $\delta x_{i}^{a}$ is the sifting distance of particle $i$ and $D$ is the diffusion coefficient, which reads [36],

$$ D = - Ah\left| v \right|_{i} dt, $$

(17)

with a free parameter $A$ set to 2 in this study.

2.2 Boundary Conditions

The wall boundary condition applied in our model is the dynamic boundary conditions (DBCs) [37] where particles representing the wall are organized in a staggered arrangement and satisfy the same equations as the fluid particles, but their position and velocity are prescribed. The advantages of the DBC include the straightforward computational implementation and the treatment of arbitrary complex geometries. In order to achieve a steady flow, open boundary conditions available in DualSPHysics [38] that served as an inlet and outlet were assigned in our model. The algorithm is based on the use of buffer layers to discretize the open boundaries, and the assignment of flow variables in these regions can be made in two ways: either imposing a flow quantity a priori or extrapolating quantities from the domain interior. The latter is carried out by computing a corrected SPH interpolation at a reference ghost node located inside the fluid domain and then a linearized projection of the interpolated quantity to the corresponding buffer particle [38]. As the implementation of boundary conditions into SPH is not the focus of the present research, for more information the reader is referred to more recent work [38, 39].

2.3 Model configuration

The numerical model domain was set up according to the test section of the original flume experiment by Papanicolaou et al. [2] with dimensions of (L_x, L_y, L_z) = (1.14, 0.3, 0.247) m. A schematic view of the 3-D model to aid visual representation is shown in Fig. 1.

Fig. 1

Model set-up (not to scale), indicating the position and location of the measured flow velocity profiles

A bottom wall was placed at the bottom of the domain with dimensions of L_x = 0.54 m, L_y = 0.3 m and L_z = 0.003 m. Two solid blocks of dimensions L_x’ = 0.3 m, L_y’ = 0.3 m and L_z’ = 0.019 m were positioned at the bottom in the upstream and downstream bed regions of the domain. The rough bed is represented by a total of 522 spheres with a diameter of D_r = 0.019 m that served as roughness elements shown in Fig. 1 to perturbate the flow and were placed between both blocks. Notice that for comparability of the numerical model with the experiments of Papanicolaou et al. [2], the roughness elements were placed in the same hexagonal arrangement as given in Papanicolaou et al. [2]. This was achieved by using a lattice where the coordinates indicating the centre of each sphere were arranged in hexagonal close packing. Each sphere has radius r and is arranged in a hexagonal pattern such that alternate rows of spheres in the x direction have their centres shifted by a distance r in the x direction. This procedure was implemented in MATLAB and repeated until the entire bottom was filled with spheres.

At the top of the roughness elements, a sphere with a larger diameter of D_s = 0.055 m was positioned at the centre of the domain at x = 0.57 m, y = 0.152 m and z = 0.0275 m. Note that both the smaller spheres and the large sphere follow exactly the dimensions as those given in Papanicolaou et al. [2]. At the front and back of the domain, the two solid walls in the dimensions of L_x = 1.14 m, L_y = 0.003 m and L_z = 0.247 m were positioned, whereas the top of the domain was a free surface. Notice that the porpoise using solid front and back walls was to replicate the settings of the flume experiment of Papanicolaou et al. [2].

All solid boundaries, including the spheres on the rough bed, are represented in the DualSPHysics solver using immobile boundary particles (nfixed). After generating the walls, the fluid particles (nfluid) were generated. In order to generate the hexagonal arrangement of the rough bed, we used the ‘sphere’ command in the program GenCase, by which the position of the sphere as well as the corresponding radius can be pre-defined. After defining all elements required to set up the case, the program GenCase was run that transfers the elements into particles placed on a regularly spaced lattice. Further parameters of potential interest that were used to run the presented simulations are summarized in Table 1. Notice that particle shifting to maintain free-surface flows was enabled, according to Mokos et al. (2016).

Table 1 Summary of the simulation parameters

In order to drive the flow, an inlet boundary condition [38] was applied at the left side of the domain, whereas an outlet boundary condition was applied at the right side of the domain (Fig. 1). According to Tafuni et al. [38], particles are created at the inflow, where a particle has just crossed the inlet and becomes a fluid particle. A new inflow particle is then created distance 2 h + ε from the inlet. This zone is referred as a buffer. Conversely, particles that populate the buffers adjacent to the outlet of the computational domain are called outflow particles, where the buffer is 2 h + ε located from the outlet. Here , 2 h is the kernel radius adopted in DualSPHysics, and ε > 0 is an arbitrarily small constant that ensures full kernel support for the particles in the near proximity of an inlet or outlet. The reader is referred to Tafuni et al. [38] for more details about the implementation of open boundary conditions in DualSPHysics.

To reproduce the flow of the flume experiment the same pattern of the flow velocity values with respect to height at the inlet of the domain, a Dirichlet parabolic velocity profile was fitted to the dataset by Papanicolaou et al. [2]. The discretization of geometric shapes such as the walls the spheres and the fluid, into particles, as well as the definition of the boundary conditions was realized with the program GenCase 4.0 that is part of the DualSPHysics software environment (version 4.3).

A spatial convergence study using three different particle spacings (dp) was set up in order to investigate the effect of particle spacing. The cases include three resolutions: dp/D_s = 0.073, 0.055 and 0.036 (Table 1) giving 14, 18 and 27 particles across the sphere, respectively. Notice that the Reynolds number was computed according to the equation ${\text{Re}}_{eq} = 5{*}\left( {\frac{{{\text{Ma}}}}{{\upalpha }}} \right)\left( {\frac{{D_{s} }}{dp}} \right),$ where Ma is the Mach number, α the artificial viscosity and dp the particle spacing.

In total, the numerical model simulated t = 5.0 s of physical time resulting in 120.4 h of computational time for the highest resolution case (dp / D_s = 0.036). The models were computed on an Intel Xeon E5 server with a NVidia Tesla K20c graphics card.

2.4 Post-processing of simulation data

In order to reproduce sufficiently the experimental set-up of Papanicolaou et al. [2], only the data obtained within the range between the upstream and downstream blocks and above the roughness element spheres of the numerical model (Fig. 1) were considered for further analysis.

The discharge as well as the average flow velocity of fluid entering the domain and leaving the domain was quantified using the DualSPHysics post-processing tool Flowtool4. For the measurement of the discharge, two measurement volumes with a defined volume of 0.64 m³ each at the upstream and downstream parts of the domain were defined. These served as reference volumes by which the discharge and average velocity entering and leaving the domain were monitored over the entire simulation time.

To compare the numerical simulation results with the measurements obtained from the flume experiment, 11 profiles were extracted from the model. All profiles were positioned at the same locations as the experimental set-up as shown in Table 2 given in Papanicolaou et al. [2].

Table 2 Numerical profile locations according to Papanicolaou et al. [2]

The velocity data computed by the SPH model were processed with the program MeasureTool4 that is part of the DualSPHysics post-processing suite (see Crespo et al. [15]). Within the program, a regular raster can be defined to which the SPH data are mapped using a Wendland kernel. Therefore, data in three, two and one dimensions can be post-processed by interpolation onto a Eulerian grid. In our case, one-dimensional profiles were defined (Table 1). The spacing between the points defining the profiles was 0.001 m. As a result, the velocity data of the flume experiment and the data of the numerical runs were plotted over height ranging from z = 0.0 m to 0.2 m to allow the validation of the numerical model results.

To compare the velocity data of the numerical model to the data of Papanicolaou et al. [2], the coefficient of determination $R^{2} = 1 - SS_{reg} /SS_{tot}$ was computed. Therefore, the sum of squares of the regression ($SS_{reg}$), where $SS_{reg} = \sum \left( {y_{i} - f_{i} } \right)^{2}$ of the model data ($y_{i}$) and the experimental data ($f_{i}$) and the sum of squares total ($SS_{tot}$), where $SS_{tot} = \sum \left( {y_{i} - y} \right)^{2}$ of the model ($y_{i}$) and the mean values of the model ($y$) were calculated.

In order to extract the drag and lift forces affecting the sphere and hence to allow for comparison with the data given in Papanicolaou et al. [2], the pressure drag was computed for the sphere and the roughness spheres at the bottom of the domain. To aid visualization, an isosurface was generated for both the particles of the representing the sphere and the roughness elements using the DualSPHysics post-processing tool Isosurface 4.0.

3 Results

3.1 Time-averaged flow field

Figure 2 shows the time-averaged flow velocity magnitude for dp / D_s = 0.055 for two slices at the centre of the sphere with a cross section of the model domain at y = 0.15 m and z = 0.035 m resulting in (a) a y-normal side view and (b) a z-normal top view, respectively. Here, time averaging refers to an averaging of instantaneous velocity over t = 2.5 to 5.0 s. The equivalent Reynolds number is Re = 733 with a freestream velocity of 0.9 m/s. A typical stagnation point in front of the sphere is developed as shown in the velocity field (Fig. 2a, b). The velocity profiles are discussed in detail in the following.

Fig. 2

Time-averaged velocity magnitudes: a top, y normal slice, b bottom, z normal slice cut through the model domain

The time-averaged velocity distribution of the discharge entering and leaving the domain was quantified in a time interval between t = 2.5 s and 5 s. The inflow into the domain was 0.0544 m³/s (i.e. 54 l/s), and the measurement at the outflow also revealed the same value, indicating a constant transport over time. The corresponding average velocity value in the downstream entering and leaving the domain was 0.7559 m/s and is referred here as U₀.

In the region close to the sphere at x = 0.4 to 0.45 m, the flow velocity magnitude is minimum at the front of the sphere approaching zero. In contrast, the flow velocity values at the top of the sphere were highest (x = 0.42 to 0.45 m and z = 0.06 to 0.08 m). The values in this region ranged between 0.91 (1.2 U/U₀) and 0.95 m/s (1.6 U/U₀). These values indicate generally higher flow speeds of 0.05 m/s compared to those measured in the free stream. From the lee side of the sphere towards the end of the domain, a different pattern becomes evident at x = 0.45 to 0.55 m. In the region immediately behind the sphere, the flow speed values appeared lower between z = 0.0 m and z = 0.075 m (Fig. 2a). The flow speed values in this low-velocity region averaged 0.18 m/s (0.24 U/U₀). Further downstream, the values appeared to approach to the flow velocity values that were measured in the freestream region.

A similar pattern of the flow velocity value distribution can be observed in the case where a slice of the time-averaged velocity field in the z-normal direction through the model domain is shown in Fig. 2b. In the upstream region of the sphere, a relatively homogeneous flow field can be observed. This is characterized by velocity values in the range of 0.45 m/s (0.6 U/U₀) at x = 0.3 to 0.4 m. In the region close to the front of the sphere, a region with relatively low flow velocity values is also evident. The flow velocity values in this area decrease towards 0.04 m/s. In the regions close to the lower and upper sides of the sphere (x = 0.425 to 0.45 m), two regions showing high velocity values become visible. Values within the two regions were 0.65 m/s (0.86 U/U₀). This indicates an increase of 0.2 m/s (0.26 U/U₀) as compared to the values obtained at the sides of the domain. In the region immediately behind the sphere from x = 0.45 m towards the end of the domain, a characteristic patch of low velocity values in the range of 0.25 m/s (0.33 U/U₀) becomes evident. Similar to Fig. 2a, further downstream direction, the velocity values approach the freestream values. In summary, regions that are characterized by low flow velocity values are located close to the front and at the rear of the sphere, whereas regions indicating high velocity values are located close to the top and the sides of the sphere.

3.2 Fluid velocity profiles

In total, 11 profiles corresponding to P1–P11 in Fig. 1 were extracted from the centreline of the model domain along the downstream direction in a y-normal transect. The exact positions of the profiles are indicated in Fig. 3 for the three resolutions of dp/D_s and given in Table 2. The grey lines in Fig. 3 indicate all instantaneous velocity profiles of t = 4.95 s that were measured. The solid black lines in Fig. 3 indicate time-averaged velocity profiles averaged over the duration of 2.5 to 5 s. This procedure was applied to all tested resolutions in Table 1. The solid grey line refers to the average of the instantaneous velocity profiles (u_x) extracted from the numerical model. The solid dots indicate the data points obtained with the laboratory-based flume experiment measured by Papanicolaou et al. [2].

Fig. 3

Fluid profiles extracted at 11 positions from the centre of the model domain as a transect. Solid black lines indicate the time-averaged velocity values, and the grey-coloured solid lines indicate the instantaneous velocity values obtained from the numerical model. The black dotted lines indicate the data point of the laboratory-based flume experiment of Papanicolaou et al. [2]. The R² refers to the coefficient of determination derived from data of the model and the laboratory-based flume experiment

With all resolutions, in the case of profile 1 the velocity values show a parabolic distribution that is characterized by values in the range of 0.92 m/s (at z = 0.2 m). In the region towards the bottom of the domain, the flow velocity values tend to decrease successively in a parabolic manner and approach zero at the bottom. The distribution of the flow velocity values shown in profiles 2–4 indicates a similar pattern. In the case of profile 5, the shape of the profile appeared different compared to the previous profiles. Here, a decrease in the flow velocity values at the front of the sphere can be observed (z = 0.06 and z = 0.02 m). Herein, the velocity values lowered by 0.4 m/s as compared to the flow velocity values measured further upstream.

In the case of profile 6, measured directly above the top of the sphere a different trend can be observed (z = 0.11 m and z = 0.06 m). Generally, higher velocity values compared to the upstream reach peaking at maximum of 1.0 m/s (z = 0.065 m) were measured. In the direct vicinity of the sphere, the flow velocity values appear to level off to 0.0 m/s in a parabolic-shaped manner.

In the case of the profiles that were measured in the rear of the sphere (i.e. profiles 7–11), a similar trend to the parabolic-shaped profiles that were measured upstream of the sphere becomes evident (i.e. between z = 0.2 m and 0.06 m). However, below z = 0.06 m a different pattern can be observed compared to the distribution of the flow speed values in the upstream reach of the sphere. Generally, lower velocity values behind the sphere are observed compared to the upstream reach of the sphere. In the case of profile 7, a relatively sharp decline of the flow velocity values towards 0.0 m/s is apparent (i.e. below z = 0.06 m). The vertical extent of this pattern appears to level off further into the downstream. For example, in the case of profile 9 the flow velocities reduce from 0.92 m/s to 0.0 m/s within the range of z = 0.05 m to 0.0 m. However, in the case of profile 11, the shape of the flow velocity distribution appeared similar to the profiles extracted upstream of the sphere again. This indicates that the average flow velocity values approached the upstream distribution.

To investigate the effect of particle size, a convergence study was performed with three different particle resolutions dp / D_s. The velocity profiles appear similar in shape with regard to the far field in the upstream reach of the sphere. For a clear comparison of the flow speed values with the tested resolutions, Fig. 4 shows the velocity values of profile 10 as an example case. As the resolution increases (i.e. decreasing dp), the profiles show convergence. The shape of all profiles appears similar above 0.03 m. Slight differences in the shape become visible in the range between 0.0 and 0.03 m with regard to the tested resolutions. The medium- and high-resolution cases show a closer agreement approaching the bottom wall. This indicates that in the near-wall region, a high resolution (with more particles) is needed to resolve the flow than in the freestream. In Fig. 3, minor differences appear at the bottom of the domain in the rear of the sphere observable in profiles 9–11. In the high-resolution case, the instantaneous velocity values (grey lines in Fig. 3) appear to level off to higher and lower values of 0.04 m/s with regard to the low- and medium-resolution cases, whereas the time-averaged profile represents the shape of the low- and medium-resolution cases. These results indicate that a resolution of dp / D_s = 0.055 is sufficient to resolve the time-averaged flow.

Fig. 4

Fluid profiles extracted from position 10. The results show the time-averaged velocity values the with regard to the low-, medium- and high-resolution cases obtained with the numerical model. The black dotted lines indicate the data point of the laboratory-based flume experiment of Papanicolaou et al. [2]

3.3 Drag and lift forces affecting the sphere

The respective forces that were computed by integration of the pressure drag and lift values affecting the sphere are shown in Fig. 5, as a time series and were plotted over 2.5 s. With respect to the drag force, the values fluctuated between C_d =  − 0.1535 (F_drag =  − 0.029 N) and C_d = 1.7436 (F_drag = 0.3314 N) and averaged into a time-averaged drag force of C_d = 0.73 (F_drag = 0.1393 N). In terms of the lift forces, a similar pattern of relatively constant pressure forces over time can be observed but was generally stronger pronounced. The corresponding values showed a minimum of C_l =  − 7.5567 (F_lift = −1.4363 N) and a maximum of − 5.7254 (F_lift =  − 1.0882 N), whereby the time-averaged value was C_l =  − 6.5687 (F_lift =  − 1.2495 N).

Fig. 5

Time series of the drag forces (top panel) and lift forces (bottom panel) affecting the sphere

3.4 Streamlines and pressure distributions

The drag and lift components acting on the spheres were computed by multiplying pressure acting over a surface area by the x and z components of the surface normal and integrating this over the entire surface. The flow velocity distribution, lift and drag pressure distributions at the surface of the sphere and the roughness elements are shown with the streamlines in Fig. 6.

Fig. 6

The streamlines and the flow field cut through the domain are shown as a y—normal oriented slice: a pressure drag values affecting the sphere as well as the roughness elements, b the pressure lift distribution and streamlines

The streamlines show a pattern of parallel oriented lines before the flow encountered the region close to the sphere facing in the downstream direction at approximately x = 0.3–0.4 m. As expected, the streamlines are deflected towards the cross-stream (y direction) and vertical direction surrounding the sphere. This can be observed in the range between x = 0.4 m and x = 0.425 m. Towards the rear of the sphere, the streamlines were arranged in chaotic pattern (x = 0.4 m to x = 0.5 m). Further downstream direction, the streamlines re-join the freestream direction. This also matches with the distribution of the flow velocity values shown in the profiles that indicate the same trend.

The distribution of the lift and drag pressures affecting the sphere as well as the roughness elements at the bottom of the domain is shown in Fig. 5a, b. The results of the pressure drag indicated a pattern of negative pressure drag values imposing the front of the sphere that is facing the flow of − 2400 Pa peak drag pressure. In the rear of the sphere that is oriented with the flow, highest pressure drag values were observed of 2700 Pa peak lift pressure.

A similar pattern for the large sphere is observed in Fig. 5b for the individual roughness elements, where in front of the sphere relatively high negative pressure drag values are present, and in the rear of the spheres, there are relatively high positive pressure drag values. The value range was generally lower as obtained from the single sphere and ranged between − 1288 Pa and 1400 Pa.

In the case of the pressure lift distribution, a maximum of negative pressure lift values is located at the top of the sphere (Fig. 5b). The values in this region range between − 2000 and 2200 Pa. In the region towards the bottom of the sphere towards the roughness elements, values were at a maximum which is indicated by values of 2700 Pa. A transition between both maxima is evident in the at midpoint of the sphere, where the values ranged between − 1288 Pa and 1400 Pa. In terms of the pressure lift value distribution of the roughness elements, values were only measured in the upper half of the spheres where the highest negative values were measured at top of the spheres. Nonetheless, values were generally higher as in the case of the sphere peaking at − 2400 Pa at the top.

4 Discussion

The distribution of the velocity values that were derived from the numerical model shown in Fig. 3 stands in good agreement with the flow velocity measurements of Papanicolaou et al. [2]. In particular, the coefficient of determination in the order of R² = 0.9 (Fig. 3) indicated the effectiveness the models of all profiles, whereby the parabolic shape of the velocity values that was assigned at the inlet is represented well with respect to height. Minor derivations of the numerically obtained data points to the experimentally obtained dataset are evident in the range close to the bed in the first 0.02 m, in the cases of profiles 3, 4, 9 and 11. This might be due to a low sampling rate of 0.01 s chosen to derive a time-averaged value. Another cause for this issue could be due to the fitting of the experimental data of profile 1 that was then assigned to the inlet of the domain. In particular, a slight difference in the fit and the original dataset appears within the first 0.02 m above the bottom of the model domain. Therefore, it is expected that a higher sampling rate with respect to a longer simulation time as well as a more idealized fit would smooth out these minor differences. Another issue that is often the cause of such minor derivations is the dynamic boundary conditions that are implemented within the present version of DualSPHysics. This issue becomes evident when comparing the profiles that were processed directly above the sphere (Fig. 3, profile 6) of the low-resolution case and the high-resolution case. The profile of the low-resolution case indicates a slightly lower coefficient of confidence of 0.8948, whereas the high-resolution case indicates a higher coefficient of confidence of 0.9321, indicating a higher resolution of the flow. Nonetheless, following Crespo et al. [15], as the particle resolution increases, the error due to the dynamic boundary conditions reduces becomes neglectable. This is clearly evident in Fig. 4 where the profiles of the medium- and high-resolution case show a similar pattern. Another explanation for this discrepancy could be related to the absence of a turbulence model, which will be the topic of a future study. Despite these minor differences, a satisfactory match can be observed in the cases of the profiles that were extracted in the rear of the sphere demonstrating a high fidelity of the data produced by the numerical model being capable of reproducing the measurements of a physical experiment.

A direct comparison of the C_d values of the flume experiment and our results appears difficult as the C_d values given in Papanicolaou et al. [2] were experimentally estimated and based on the assumption of a C_d given in Coleman [40], suggesting a value C_d = 0.48 and a characteristic velocity of 0.40 m/s, which results in F_D = 0.076 N. However, the C_D values established from our simulation were derived from a direct integration of the pressure values acting normal to the sphere. This procedure hence may explain the slight difference to the C_d value derived from the SPH model.

A similar numerical attempt to demonstrate the fidelity of a numerical model was done by Lui et al. (2017). Their model is based on the finite volume method (FVM) and solves the flow equations using a Eulerian, i.e. grid-based approach. In accordance with the data derived from our model, their results also stand in a good agreement with the flow velocity values derived with the flume experiment of Papanicolaou et al. [2]. Also, a good match of the pattern of the streamlines behind the sphere and our numerical results demonstrates the fidelity of the results of the Lagrangian modelling approach employed here. A significant difference in their modelling approach is that the FVM model uses periodic boundary conditions to simulate a recycling flow, whereas our model uses absorbing inlet and outlet boundary conditions with a mapped fluid profile at the inlet to replicate the flume experiment of Papanicolaou et al. [2]. Furthermore, the runtime of 120 h for the finest resolution case, dp / D_s = 0.036, performed on a desktop GPU, is significantly shorter than the 10³ h required by Lui et al. (2017).

With respect to the hydrodynamics, the SPH model reproduces well the typical flow features that are often described in the literature in terms of the flow past a bluff bodies [4, 5, 41, 42]. This includes a laminar flow that is characterized by parallel streamlines in the region upstream of the sphere. A transition into an unsteady flow that is characterized by a flow separation zone, which is indicated by deflected streamlines, into the region in which the streamlines appeared highly distorted and chaotic, which is also referred as the wake or shedding zone, is also well reproduced by the numerical model.

Despite the promising agreement with the experimental dataset, there are some limitations of the present study. It is clearly acknowledged that the present model does not have the accuracy of high-fidelity numerical large-eddy simulation (LES) models with resolved fluid flow. However, the presented model results are sufficient to reproduce the flow around a sphere, a phenomenon that is widely observed in nature. Turbulence was not investigated in the present study, and this will serve as the topic of a future study.

5 Conclusions

In this study, a numerical simulation using SPH was performed to demonstrate the feasibility of a Lagrangian numerical scheme capable of simulating the flow past a sphere on a rough bed. This test case appeared ideal for validation of the numerical results due to the fact that such an experiment represents a simplified scenario by which the most common flow features can be observed that one would also expect when it comes to the investigation of the hydrodynamics in the vicinity of a boulder that sits on a rough bed. The flow over a rough bed has been investigated before in SPH [24,25,26], none of these studies include a rough bed in combination with objects/body near the bed in the presence of a free-surface flow. This is a new and novel application with insight into the fluid mechanics only made possible by using the recent versions of DualSPHysics. The numerical model results show promising agreement with the flow speed data of Papanicolaou et al. [2]. Moreover, a qualitative comparison with the data of Lui et al. (2017) further indicated the fidelity of the SPH model by indicating similar distribution of the flow speed values as well as a good match with the deflection of the streamlines within the rear of the sphere. Moreover, typical flow features that have been reported in the literature such as the flow separation leading to the generation of a wake are reproduced well by the numerical model. Another goal of this study was to demonstrate the applicability of a SPH model by the open-source software environment DualSPHysics as a suitable tool to perform such type of simulations. Moreover, it was shown that due to the GPU-based approach to solve the flow equations implemented in the DualSPHysics environment and hence appeared as an ideal software environment for single users without access to high-performance computation faculties but are in demand to run high-resolution models on a small budget or time frame. It is planned to use the parameters of the present study as a solid basis for future simulations such as the simulation of the entrainment of moving boulders.